scalar_sweeper.hpp

Go to the documentation of this file.

#ifndef _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
#define _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_

#include <complex>
#include <vector>
using namespace std;

#include <pfasst/logging.hpp>
#include <pfasst/encap/imex_sweeper.hpp>
#include <pfasst/encap/vector.hpp>


namespace pfasst
{
  namespace examples
  {
    namespace scalar
    {
      template<typename time = pfasst::time_precision>
      class ScalarSweeper
        : public encap::IMEXSweeper<time>
      {
        private:
          typedef encap::Encapsulation<time> encap_type;

          typedef encap::VectorEncapsulation<complex<double>> complex_vector_type;

          complex<double> lambda, u0;

          const complex<double> i_complex = complex<double>(0, 1);

          double error;

          size_t n_f_expl_eval, n_f_impl_eval, n_impl_solve;

        public:
          ScalarSweeper(const complex<double>& lambda, const complex<double>& u0)
            :   lambda(lambda)
              , u0(u0)
              , error(0.0)
              , n_f_expl_eval(0)
              , n_f_impl_eval(0)
              , n_impl_solve(0)
          {}

          virtual ~ScalarSweeper()
          {
            ML_LOG(INFO, "Final error:                   " << this->error);
            ML_LOG(INFO, "Number of explicit evaluations:" << this->n_f_expl_eval);
            ML_LOG(INFO, "Number of implicit evaluations:" << this->n_f_impl_eval);
            ML_LOG(INFO, "Number of implicit solves:     " << this->n_impl_solve);
          }

          void echo_error(time t)
          {
            auto& qend = encap::as_vector<complex<double>, time>(this->get_end_state());

            complex_vector_type qex(qend.size());

            this->exact(qex, t);
            double max_err = abs(qend[0] - qex[0]) / abs(qex[0]);
            ML_LOG(INFO, "err:" << max_err);
            this->error = max_err;
          }

          double get_errors()
          {
            return this->error;
          }

          void post_predict() override
          {
            time t  = this->get_controller()->get_time();
            time dt = this->get_controller()->get_step_size();
            this->echo_error(t + dt);
          }

          void post_sweep() override
          {
            time t  = this->get_controller()->get_time();
            time dt = this->get_controller()->get_step_size();
            this->echo_error(t + dt);
          }

          void exact(complex_vector_type& q, time t)
          {
            q[0] = this->u0 * exp(this->lambda * t);
          }

          void exact(shared_ptr<encap_type> q_encap, time t)
          {
            auto& q = encap::as_vector<complex<double>, time>(q_encap);
            this->exact(q, t);
          }

          void f_expl_eval(shared_ptr<encap_type> f_encap,
                           shared_ptr<encap_type> q_encap, time t) override
          {
            UNUSED(t);
            auto& f = encap::as_vector<complex<double>, time>(f_encap);
            auto& q = encap::as_vector<complex<double>, time>(q_encap);

            // f_expl = multiply with imaginary part of lambda
            f[0] = this->i_complex * imag(this->lambda) * q[0];

            this->n_f_expl_eval++;
          }

          void f_impl_eval(shared_ptr<encap_type> f_encap,
                           shared_ptr<encap_type> q_encap, time t) override
          {
            UNUSED(t);
            auto& f = encap::as_vector<complex<double>, time>(f_encap);
            auto& q = encap::as_vector<complex<double>, time>(q_encap);

            // f_impl = multiply with real part of lambda
            f[0] = real(this->lambda) * q[0];

            this->n_f_impl_eval++;
          }

          void impl_solve(shared_ptr<encap_type> f_encap,
                          shared_ptr<encap_type> q_encap, time t, time dt,
                          shared_ptr<encap_type> rhs_encap) override
          {
            UNUSED(t);
            auto& f = encap::as_vector<complex<double>, time>(f_encap);
            auto& q = encap::as_vector<complex<double>, time>(q_encap);
            auto& rhs = encap::as_vector<complex<double>, time>(rhs_encap);

            // invert f_impl = multiply with inverse of real part of lambda
            double inv = 1.0 / (1.0 - double(dt) * real(this->lambda));
            q[0] = inv * rhs[0];

            // compute f_impl_eval of q[0], i.e. multiply with real(lambda)
            f[0] = real(this->lambda) * q[0];

            this->n_impl_solve++;
          }
      };
    }  // ::pfasst::examples::scalar
  }  // ::pfasst::examples
}  // ::pfasst

#endif  // _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_