bdf_time_integration.hpp

Go to the documentation of this file.

#pragma once
 
 
 
#include <deal.II/lac/la_parallel_vector.h>
 
#include <meltpooldg/linear_algebra/newton_raphson_solver.hpp>
#include <meltpooldg/linear_algebra/preconditioner.hpp>
#include <meltpooldg/linear_algebra/preconditioner_trilinos_wrapper.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 <functional>
#include <memory>
#include <vector>
 
namespace MeltPoolDG::TimeIntegration
{
  inline static constexpr std::array<TimeIntegratorSchemes, 6> bdf_supported_schemes{
    {TimeIntegratorSchemes::bdf_1,
     TimeIntegratorSchemes::bdf_2,
     TimeIntegratorSchemes::bdf_3,
     TimeIntegratorSchemes::bdf_4,
     TimeIntegratorSchemes::bdf_5,
     TimeIntegratorSchemes::bdf_6}};
 
  template <int dim, typename number>
  class BDFIntegrator final : public TimeIntegratorBase<number>
  {
  public:
    using VectorType = dealii::LinearAlgebra::distributed::Vector<number>;
 
    using JacobianType =
      std::function<void(number time_step, VectorType &src, const VectorType &dst)>;
 
    using ResidualType = std::function<void(number            time,
                                            number            time_step,
                                            const VectorType &src,
                                            VectorType       &dst,
                                            const VectorType &old_solution)>;
 
    using DistributeConstraintsType = std::function<void(VectorType &dst)>;
 
    struct SolverFunctions
    {
      JacobianType compute_jacobian;
 
      ResidualType compute_residual;
 
      DistributeConstraintsType distribute_constraints;
    };
 
    explicit BDFIntegrator(
      const TimeIntegratorData<number>         &time_integrator_data,
      const SolverFunctions                    &solver_functions,
      Preconditioner<dim, VectorType, number> &&preconditioner_in =
        Preconditioner<dim, VectorType, number>(IdentityPreconditioner<dim, VectorType, number>()));
 
    unsigned int
    required_solution_history_size() const override;
 
    void
    reinit(const VectorType &vector_template) override;
 
    void
    reinit(const SolutionHistory<VectorType> &solution_history) override;
 
    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 &)>
                        &stage_pre_processing,
                      const std::function<void(number, number, VectorType &, const VectorType &)>
                        &stage_post_processing) override;
 
  private:
    ResidualType compute_residual;
 
    JacobianType compute_jacobian;
 
    DistributeConstraintsType distribute_constraints;
 
    std::unique_ptr<NewtonRaphsonSolver<number, VectorType>> solver;
 
    Preconditioner<dim, VectorType, number> preconditioner;
 
    struct
    {
      number rhs;
      std::vector<number> old_solutions;
    } bdf_weights;
 
    VectorType summed_old_solution;
 
    bool preconditioner_update_flag = false;
 
    unsigned n_steps_performed = 0;
 
    void
    compute_weighted_old_solution_sum(SolutionHistory<VectorType> &solution_history);
 
    void
    set_up_bdf_parameters(const unsigned int bdf_scheme);
  };
} // namespace MeltPoolDG::TimeIntegration