vdp_sweeper.hpp

Go to the documentation of this file.

#ifndef _EXAMPLES__VANDERPOL__VDP_SWEEPER_HPP_
#define _EXAMPLES__VANDERPOL__VDP_SWEEPER_HPP_

#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

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


namespace pfasst
{
  namespace examples
  {
    namespace vdp
    {
      template<typename time = pfasst::time_precision>
      class VdpSweeper
        : public encap::ImplicitSweeper<time>
      {
        private:
          typedef encap::Encapsulation<time> encap_type;

          typedef encap::VectorEncapsulation<double> real_vector_type;

          double nu;

          double x0, y0;

          size_t newton_maxit;

          double newton_tol;

          size_t n_f_impl_eval, n_impl_solve, n_newton_iter;

          fstream output_file;

          double error;

        public:
          VdpSweeper(double nu, double x0, double y0)
            :   nu(nu)
              , x0(x0)
              , y0(y0)
              , newton_maxit(50)
              , newton_tol(1e-12)
              , n_f_impl_eval(0)
              , n_impl_solve(0)
              , n_newton_iter(0)
          {
            this->output_file.open("./vanderpol.txt", ios_base::out);

            // The file should also contain the initial value, so perform a first write here
            this->output_file << x0 << "    " << y0 << endl;
          }

          virtual ~VdpSweeper()
          {
            ML_LOG(INFO, "Number of implicit evaluations:" << this->n_f_impl_eval);
            ML_LOG(INFO, "Number of implicit solves:     " << this->n_impl_solve);
            this->output_file.close();
          }

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

            real_vector_type qex(qend.size());

            this->exact(qex, t);

            if (this->nu==0)
            {
              double max_err = max(abs(qend[0] - qex[0])/abs(qex[0]) , abs(qend[1]-qex[1])/abs(qex[1]) );
              ML_LOG(INFO, "error:" << max_err);
              this->error = max_err;
            }
            this->output_file << qend[0] << "    " << qend[1] << endl;
          }

          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);
          }

          // I don't know why Doxygen is eating up so many backslashes here ...
          void exact(real_vector_type& q, time t)
          {
            if (this->nu==0)
            {
              q[0] = this->y0*sin(t) + this->x0*cos(t);
              q[1] = -this->x0*sin(t) + this->y0*cos(t);
            }
            else
            {
              // For the nonlinear case, there is no analytic solution for the vdP oscillator, so
              // simply return the initial value
              q[0] = this->x0;
              q[1] = this->y0;
            }
          }

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

          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<double, time>(f_encap);
            auto& q = encap::as_vector<double, time>(q_encap);

            // x' = y
            f[0] = q[1];
            // y' = nu(1-x^2)*y - x
            f[1] = this->nu*(1.0-q[0]*q[0])*q[1] - q[0];

            this->n_f_impl_eval++;
          }

          // I don't know why Doxygen is eating up so many backslashes here ...
          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<double, time>(f_encap);
            auto& q = encap::as_vector<double, time>(q_encap);
            auto& rhs = encap::as_vector<double, time>(rhs_encap);

            // I don't know why Doxygen is eating up so many backslashes here ...
            // initialize residual
            double residual = this->newton_tol + 1.0;
            size_t iter     = 0.0;

            // Initial value for q is just rhs: For small dt, P is approximately the identity
            q[0] = rhs[0];
            q[1] = rhs[1];

            // NEWTON ITERATION: q_new = q + inv(J(q))*(-f(q))
            do {
              // Store -P(q) in f0, f1
              double f0 = -( q[0] - dt*q[1] - rhs[0] );
              double f1 = -( q[1] - dt*( this->nu*(1-q[0]*q[0])*q[1]-q[0]) - rhs[1] );

              double a = dt * q[0] * q[0] - dt + 1.0;
              double b = -2.0 * dt * this->nu * q[0] * q[1] - dt;
              double c = 2.0 * this->nu * q[0] * q[1] * dt * dt + dt * dt + dt * q[0] * q[0] - dt + 1.0;

              // Compute inv(J(q))*(-f(q)) and store it in f
              f[0] = (1.0/c)*( a*f0 + dt*f1 );
              f[1] = (1.0/c)*( b*f0 + f1 );

              // Perform Newton update q_new = q + inv(J(q))*(-f(q))
              q[0] += f[0];
              q[1] += f[1];

              // Compute relative residual; note that f = q_new - q is the update
              residual = fmax( abs(f[0]), abs(f[1]) ) / fmax( abs(q[0]), abs(q[1]) );

              iter++;
              this->n_newton_iter++;

            } while ( (iter<this->newton_maxit) && (residual>this->newton_tol) );

            // Say something, only if the residual tolerance has not been reach in the maximum
            // number of iterations
            if (residual >=this->newton_tol)
            {
              cout << "Newton failed to converge: res = " << scientific << residual << " -- n_iter = " << iter << " of maxit = " << this->newton_maxit << endl;
            }

            // Set f to f(q)
            f[0] = q[1];
            f[1] = this->nu*(1.0-q[0]*q[0])*q[1] - q[0];

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

#endif  // _EXAMPLES__VANDERPOL__VDP_SWEEPER_HPP_