<commit_before>/*
 * Advection/diffusion sweeper.
 */

#ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
#define _ADVECTION_DIFFUSION_SWEEPER_HPP_

#include <complex>
#include <vector>
#include <cassert>
#include <ostream>

#include <pfasst/encap/imex_sweeper.hpp>

#include "fft.hpp"

#define PI     3.1415926535897932385
#define TWO_PI 6.2831853071795864769

using namespace std;

template<typename time = pfasst::time_precision>
class AdvectionDiffusionSweeper
  : public pfasst::encap::IMEXSweeper<time>
{
    typedef pfasst::encap::Encapsulation<time> Encapsulation;
    typedef pfasst::encap::VectorEncapsulation<double> DVectorT;

    FFT fft;

    vector<complex<double>> ddx, lap;

    double v  = 1.0;
    time   t0 = 1.0;
    double nu = 0.02;
    size_t nf1evals = 0;

  public:
    AdvectionDiffusionSweeper(size_t nvars)
    {
      ddx.resize(nvars);
      lap.resize(nvars);
      for (size_t i = 0; i < nvars; i++) {
        double kx = TWO_PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
        ddx[i] = complex<double>(0.0, 1.0) * kx;
        lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
      }
    }

    ~AdvectionDiffusionSweeper()
    {
      cout << "number of f1 evals: " << nf1evals << endl;
    }

    void exact(shared_ptr<Encapsulation> q, time t)
    {
      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
      assert(q_cast);
      this->exact(q_cast, t);
    }

    void exact(shared_ptr<DVectorT> q, time t)
    {
      size_t n = q->size();
      double a = 1.0 / sqrt(4 * PI * nu * (t + t0));

      for (size_t i = 0; i < n; i++) {
        q->data()[i] = 0.0;
      }

      for (int ii = -2; ii < 3; ii++) {
        for (size_t i = 0; i < n; i++) {
          double x = double(i) / n - 0.5 + ii - t * v;
          q->data()[i] += a * exp(-x * x / (4 * nu * (t + t0)));
        }
      }
    }

    void echo_error(time t, bool predict = false)
    {
      shared_ptr<DVectorT> qend = dynamic_pointer_cast<DVectorT>(this->get_state(this->get_nodes().size() - 1));
      assert(qend);
      shared_ptr<DVectorT> qex = make_shared<DVectorT>(qend->size());

      exact(qex, t);

      double max = 0.0;
      for (size_t i = 0; i < qend->size(); i++) {
        double d = abs(qend->data()[i] - qex->data()[i]);
        if (d > max) { max = d; }
      }

      auto n = this->get_controller()->get_step();
      auto k = this->get_controller()->get_iteration();
      cout << "err: " << n << " " << k << " " << scientific << max
           << " (" << qend->size() << ", " << predict << ")"
           << endl;
    }

    void predict(bool initial)
    {
      pfasst::encap::IMEXSweeper<time>::predict(initial);
      time t  = this->get_controller()->get_time();
      time dt = this->get_controller()->get_time_step();
      echo_error(t + dt, true);
    }

    void sweep()
    {
      pfasst::encap::IMEXSweeper<time>::sweep();
      time t  = this->get_controller()->get_time();
      time dt = this->get_controller()->get_time_step();
      echo_error(t + dt);
    }

    void f1eval(shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t)
    {
      shared_ptr<DVectorT> f_cast = dynamic_pointer_cast<DVectorT>(f);
      assert(f_cast);
      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
      assert(q_cast);

      this->f1eval(f_cast, q_cast, t);
    }
<commit_msg>advection_diffusion: Save errors in a map: step, iter -> error.<commit_after>/*
 * Advection/diffusion sweeper.
 */

#ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
#define _ADVECTION_DIFFUSION_SWEEPER_HPP_

#include <cassert>
#include <complex>
#include <map>
#include <ostream>
#include <vector>

#include <pfasst/encap/imex_sweeper.hpp>

#include "fft.hpp"

#define PI     3.1415926535897932385
#define TWO_PI 6.2831853071795864769

using namespace std;

template<typename time = pfasst::time_precision>
class AdvectionDiffusionSweeper
  : public pfasst::encap::IMEXSweeper<time>
{
    typedef pfasst::encap::Encapsulation<time> Encapsulation;
    typedef pfasst::encap::VectorEncapsulation<double> DVectorT;

    FFT fft;

    vector<complex<double>> ddx, lap;
    map<pair<size_t, size_t>, double> errors;

    double v  = 1.0;
    time   t0 = 1.0;
    double nu = 0.02;
    size_t nf1evals = 0;

  public:
    AdvectionDiffusionSweeper(size_t nvars)
    {
      ddx.resize(nvars);
      lap.resize(nvars);
      for (size_t i = 0; i < nvars; i++) {
        double kx = TWO_PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
        ddx[i] = complex<double>(0.0, 1.0) * kx;
        lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
      }
    }

    ~AdvectionDiffusionSweeper()
    {
      cout << "number of f1 evals: " << nf1evals << endl;
    }

    void exact(shared_ptr<Encapsulation> q, time t)
    {
      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
      assert(q_cast);
      this->exact(q_cast, t);
    }

    void exact(shared_ptr<DVectorT> q, time t)
    {
      size_t n = q->size();
      double a = 1.0 / sqrt(4 * PI * nu * (t + t0));

      for (size_t i = 0; i < n; i++) {
        q->data()[i] = 0.0;
      }

      for (int ii = -2; ii < 3; ii++) {
        for (size_t i = 0; i < n; i++) {
          double x = double(i) / n - 0.5 + ii - t * v;
          q->data()[i] += a * exp(-x * x / (4 * nu * (t + t0)));
        }
      }
    }

    void echo_error(time t, bool predict = false)
    {
      shared_ptr<DVectorT> qend = dynamic_pointer_cast<DVectorT>(this->get_state(this->get_nodes().size() - 1));
      assert(qend);
      shared_ptr<DVectorT> qex = make_shared<DVectorT>(qend->size());

      exact(qex, t);

      double max = 0.0;
      for (size_t i = 0; i < qend->size(); i++) {
        double d = abs(qend->data()[i] - qex->data()[i]);
        if (d > max) { max = d; }
      }

      auto n = this->get_controller()->get_step();
      auto k = this->get_controller()->get_iteration();
      cout << "err: " << n << " " << k << " " << scientific << max
           << " (" << qend->size() << ", " << predict << ")"
           << endl;

      errors.insert(pair<pair<size_t, size_t>, double>
		    (pair<size_t, size_t>(n, k), max));
    }

    void predict(bool initial)
    {
      pfasst::encap::IMEXSweeper<time>::predict(initial);
      time t  = this->get_controller()->get_time();
      time dt = this->get_controller()->get_time_step();
      echo_error(t + dt, true);
    }

    void sweep()
    {
      pfasst::encap::IMEXSweeper<time>::sweep();
      time t  = this->get_controller()->get_time();
      time dt = this->get_controller()->get_time_step();
      echo_error(t + dt);
    }

    void f1eval(shared_ptr<Encapsulation> f, shared_ptr<Encapsulation> q, time t)
    {
      shared_ptr<DVectorT> f_cast = dynamic_pointer_cast<DVectorT>(f);
      assert(f_cast);
      shared_ptr<DVectorT> q_cast = dynamic_pointer_cast<DVectorT>(q);
      assert(q_cast);

      this->f1eval(f_cast, q_cast, t);
    }

{'fraction_non_alphanumeric': 0.16640207822196565, 'fraction_numerical': 0.023957280992928272, 'mean_word_length': 2.6, 'pattern_counts': {'":': 0, '<': 143, '<?xml version=': 0, '>': 116, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3631166', 'n_tokens_mistral': 2784, 'n_tokens_neox': 2648, 'n_words': 533}