runge_kutta.h

// This file is part of Hermes2D.
//
// Hermes2D is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// Hermes2D is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D.  If not, see <http://www.gnu.org/licenses/>.

#ifndef __H2D_RUNGE_KUTTA_H
#define __H2D_RUNGE_KUTTA_H

#include "global.h"
#include "weakform/weakform.h"
#include "function/filter.h"
#include "exceptions.h"
#include "mixins2d.h"
namespace Hermes
{
  namespace Hermes2D
  {
    // TODO LIST:
    //
    // (1) With explicit and diagonally implicit methods, the matrix is treated
    //     in the same way as with fully implicit ones. To make this more
    //     efficient, with explicit and diagonally implicit methods one should
    //     first only solve for the upper left block, then eliminate all blocks
    //     under it, then solve for block at position 22, eliminate all blocks
    //     under it, etc. Currently this is not done and everything is left to
    //     the matrix solver.
    //
    // (2) In example 03-timedep-adapt-space-and-time with implicit Euler
    //     method, Newton's method takes much longer than in 01-timedep-adapt-space-only
    //     (that also uses implicit Euler method). This means that the initial guess for
    //     the K_vector should be improved (currently it is zero).
    //
    // (3) At the end of rk_time_step_newton(), the previous time level solution is
    //     projected onto the space of the new_ time-level solution so that
    //     it can be added to the stages. This projection is slow so we should
    //     find a way to do this differently. In any case, the projection
    //     is not necessary when no adaptivity in space takes place and the
    //     two spaces are the same (but it is done anyway).
    //
    // (4) We do not take advantage of the fact that all blocks in the
    //     Jacobian matrix have the same structure. Thus it is enough to
    //     assemble the matrix M (one block) and copy the sparsity structure
    //     into all remaining nonzero blocks (and diagonal blocks). Right
    //     now, the sparsity structure is created expensively in each block
    //     again.
    //
    // (5) If space does not change, the sparsity does not change. Right now
    //     we discard everything at the end of every time step, we should not
    //     do it.
    //
    // (6) If the problem does not depend explicitly on time, then all the blocks
    //     in the Jacobian matrix of the stationary residual are the same up
    //     to a multiplicative constant. Thus they do not have to be aassembled
    //     from scratch.
    //
    //
    // (7) In practice, Butcher's tables are being transformed to the
    //     Jordan canonical form (I think) for better performance. This
    //     can be found, I think, in newer Butcher's papers or presentation
    //     (he has them online), and possibly in his book.
    template<typename Scalar>
    class HERMES_API RungeKutta :
      public Hermes::Mixins::Loggable,
      public Hermes::Mixins::TimeMeasurable,
      public Hermes::Mixins::IntegrableWithGlobalOrder,
      public Hermes::Mixins::SettableComputationTime,
      public Hermes::Hermes2D::Mixins::SettableSpaces<Scalar>,
      public Hermes::Algebra::Mixins::MatrixRhsOutput < Scalar >
    {
    public:
      RungeKutta(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces, ButcherTable* bt);

      RungeKutta(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space, ButcherTable* bt);

      void set_start_from_zero_K_vector();
      void set_residual_as_solutions();
      void set_block_diagonal_jacobian();

      ~RungeKutta();

      // Perform one explicit or implicit time step using the Runge-Kutta method
      // corresponding to a given Butcher's table. If err_vec != nullptr then it will be
      // filled with an error vector calculated using the second B-row of the Butcher's
      // table (the second B-row B2 must be nonzero in that case). The negative default
      // values for newton_tol and newton_max_iter are for linear problems.
      // Many improvements are needed, a todo list is presented at the beginning of
      // the corresponding .cpp file.
      // freeze_jacobian... if true then the Jacobian is not recalculated in each
      //                    iteration of the Newton's method.
      // block_diagonal_jacobian... if true then the tensor product block Jacobian is
      //                            reduced to just the diagonal blocks.
      void rk_time_step_newton(std::vector<MeshFunctionSharedPtr<Scalar> > slns_time_prev, std::vector<MeshFunctionSharedPtr<Scalar> > slns_time_new, std::vector<MeshFunctionSharedPtr<Scalar> > error_fns);
      void rk_time_step_newton(MeshFunctionSharedPtr<Scalar> slns_time_prev, MeshFunctionSharedPtr<Scalar> slns_time_new, MeshFunctionSharedPtr<Scalar> error_fn);

      // This is a wrapper for the previous function if error_fn is not provided
      // (adaptive time stepping is not wanted).
      void rk_time_step_newton(std::vector<MeshFunctionSharedPtr<Scalar> > slns_time_prev, std::vector<MeshFunctionSharedPtr<Scalar> > slns_time_new);
      void rk_time_step_newton(MeshFunctionSharedPtr<Scalar> sln_time_prev, MeshFunctionSharedPtr<Scalar> sln_time_new);

      void set_freeze_jacobian();
      void set_tolerance(double newton_tol);
      void set_max_allowed_iterations(int newton_max_iter);
      void set_newton_damping_coeff(double newton_damping_coeff);
      void set_newton_max_allowed_residual_norm(double newton_max_allowed_residual_norm);

      virtual void set_spaces(std::vector<SpaceSharedPtr<Scalar> > spaces);
      virtual void set_space(SpaceSharedPtr<Scalar> space);
      virtual std::vector<SpaceSharedPtr<Scalar> > get_spaces();

      void set_filters_to_reinit(std::vector<Filter<Scalar>*> filters_to_reinit);

    protected:
      void init();

      void multiply_as_diagonal_block_matrix(SparseMatrix<Scalar>* matrix_left, int num_stages,
        Scalar* stage_coeff_vec, Scalar* vector_left);

      void create_stage_wf(unsigned int size, bool block_diagonal_jacobian);

      void update_stage_wf(std::vector<MeshFunctionSharedPtr<Scalar> > slns_time_prev);

      // Prepare u_ext_vec.
      void prepare_u_ext_vec();

      Hermes::Algebra::SparseMatrix<Scalar>* matrix_left;

      Hermes::Algebra::SparseMatrix<Scalar>* matrix_right;
      Hermes::Algebra::Vector<Scalar>* vector_right;

      Hermes::Solvers::LinearMatrixSolver<Scalar>* solver;

      const WeakFormSharedPtr<Scalar> wf;

      std::vector<SpaceSharedPtr<Scalar> > spaces;
      std::vector<unsigned int> spaces_seqs;

      ButcherTable* bt;

      unsigned int num_stages;

      // For the main part equation (written on the right),
      WeakFormSharedPtr<Scalar> stage_wf_right;
      DiscreteProblem<Scalar>* stage_dp_right;

      WeakFormSharedPtr<Scalar> stage_wf_left;
      DiscreteProblem<Scalar>* stage_dp_left;

      bool start_from_zero_K_vector;
      bool block_diagonal_jacobian;
      bool residual_as_vector;

      unsigned int iteration;

      bool freeze_jacobian;
      double newton_tol;
      int newton_max_iter;
      double newton_damping_coeff;
      double newton_max_allowed_residual_norm;

      std::vector<MeshFunctionSharedPtr<Scalar> > residuals_vector;

      Scalar* K_vector;

      Scalar* u_ext_vec;

      Scalar* vector_left;

      std::vector<Filter<Scalar>*> filters_to_reinit;
    };
  }
}
#endif