one_step_theta.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/config.h>
 
#include <deal.II/lac/la_parallel_vector.h>
 
#include <meltpooldg/linear_algebra/linear_solver.hpp>
#include <meltpooldg/linear_algebra/linear_solver_data.hpp>
#include <meltpooldg/time_integration/solution_history.hpp>
#include <meltpooldg/time_integration/time_integrator_base.hpp>
#include <meltpooldg/time_integration/time_integrator_data.hpp>
 
#include <meltpooldg/utilities/matrix_type_wrapper.h>
 
#include <functional>
#include <ostream>
 
namespace MeltPoolDG::TimeIntegration
{
  inline static constexpr std::array<TimeIntegratorSchemes, 3> one_step_theta_supported_schemes{
    {TimeIntegratorSchemes::explicit_euler,
     TimeIntegratorSchemes::implicit_euler,
     TimeIntegratorSchemes::crank_nicolson}};
 
  template <typename number, typename PDEOperator>
  class OneStepTheta final : public TimeIntegratorBase<number>
  {
  public:
    using VectorType      = dealii::LinearAlgebra::distributed::Vector<number>;
    using BlockVectorType = dealii::LinearAlgebra::distributed::BlockVector<number>;
 
    OneStepTheta(const PDEOperator                &pde_operator,
                 const TimeIntegratorData<number> &time_integrator_data_in,
                 const LinearSolverData<number>   &linear_solver_data_in)
      : TimeIntegratorBase<number>(time_integrator_data_in)
      , pde_operator(pde_operator)
      , linear_solver_data(linear_solver_data_in)
    {
      switch (this->time_integrator_data.integrator_type)
        {
            case TimeIntegratorSchemes::explicit_euler: {
              Theta_ = 0.0;
              break;
            }
 
            case TimeIntegratorSchemes::implicit_euler: {
              Theta_ = 1.0;
              break;
            }
 
            case TimeIntegratorSchemes::crank_nicolson: {
              Theta_ = 0.5;
              break;
            }
            default: {
              AssertThrow(
                false,
                dealii::ExcMessage(
                  "The chosen time integration scheme is not supported by the OneStepTheta class."));
            }
        }
    };
 
    unsigned int
    required_solution_history_size() const override
    {
      return 1;
    }
 
    void
    perform_time_step(
      const number                                                                 current_time,
      const number                                                                 time_step,
      SolutionHistory<VectorType>                                                 &solution_history,
      const std::function<void(number, number, VectorType &, const VectorType &)> &pre_processing,
      const std::function<void(number, number, VectorType &, const VectorType &)> &post_processing)
      override
    {
      Assert(solution_history.size() >= required_solution_history_size(),
             dealii::ExcMessage(
               "The size of the solution history object does not fit the requirements of the "
               "chosen time integration scheme."));
      if (pre_processing)
        pre_processing(current_time,
                       time_step,
                       right_hand_side_,
                       solution_history.get_recent_old_solution());
      dt_       = time_step;
      old_time_ = current_time - time_step;
      // Create the right hand side
      right_hand_side_.reinit(solution_history.get_recent_old_solution());
      pde_operator.apply_operator(old_time_,
                                  right_hand_side_,
                                  solution_history.get_recent_old_solution(),
                                  std::function<void(unsigned int, unsigned int)>());
      right_hand_side_ *= (1.0 - Theta_);
      buffer = 0.;
      pde_operator.apply_dirichlet_boundary_operator(old_time_,
                                                     buffer,
                                                     solution_history.get_recent_old_solution());
      right_hand_side_.add(1.0 - Theta_, buffer);
      buffer = 0.;
      pde_operator.apply_dirichlet_boundary_operator(old_time_ + dt_,
                                                     buffer,
                                                     solution_history.get_recent_old_solution());
      right_hand_side_.add(Theta_, buffer);
      right_hand_side_ *= dt_;
      right_hand_side_.add(1.0, solution_history.get_recent_old_solution());
 
      if (this->monitoring_vector != nullptr)
        *this->monitoring_vector = right_hand_side_;
 
      const MatrixTypeObject<VectorType> linear_solver_matrix(
        [&](VectorType &dst, const VectorType &src) -> void {
          dst = 0.; // TODO: Initial guess
          pde_operator.apply_operator(old_time_ + dt_,
                                      dst,
                                      src,
                                      std::function<void(unsigned int, unsigned int)>());
          dst *= -1.0;
          dst *= Theta_ * dt_;
          dst.add(1.0, src);
        });
 
 
      /* TODO: add preconditioner*/
      LinearSolver::solve<VectorType>(linear_solver_matrix,
                                      solution_history.get_current_solution(),
                                      right_hand_side_,
                                      linear_solver_data);
 
      if (post_processing)
        post_processing(current_time + time_step,
                        time_step,
                        solution_history.get_current_solution(),
                        solution_history.get_current_solution());
    }
 
 
    void
    reinit(const VectorType &vector_template) override
    {
      buffer.reinit(vector_template);
    }
 
    void
    reinit(const SolutionHistory<VectorType> &solution_history) override
    {
      reinit(solution_history.get_current_solution());
    }
 
 
  private:
    number Theta_;
 
    mutable number dt_;
    mutable number old_time_;
 
    mutable VectorType old_solution_;
    mutable VectorType right_hand_side_;
    mutable VectorType buffer;
 
    const PDEOperator              &pde_operator;
    const LinearSolverData<number> &linear_solver_data;
  };
} // namespace MeltPoolDG::TimeIntegration