parabolic_module.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2026 by the ryujin authors
//
 
#pragma once
 
#include <ryujin/base/compile_time_options.h>
#include <ryujin/base/mpi_ensemble.h>
#include <ryujin/discretization/offline_data.h>
#include <ryujin/euler_poisson/electrostatic_configuration_library.h>
#include <ryujin/euler_poisson/laplace_operator.h>
#include <ryujin/euler_poisson/parabolic_system.h>
#include <ryujin/interface/parabolic_module.h>
#include <ryujin/typed/initial_values.h>
 
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/matrix_free/matrix_free.h>
 
#include <limits>
#include <memory>
#include <vector>
 
namespace ryujin
{
  namespace EulerPoisson
  {
    enum GaussLawRestartStrategy : dealii::types::boundary_id {
      no_restart,
 
      full_restart,
 
      correction,
 
      static_no_restart,
 
      static_full_restart,
 
    };
 
  } // namespace EulerPoisson
} // namespace ryujin
 
#ifndef DOXYGEN
DECLARE_ENUM(
    ryujin::EulerPoisson::GaussLawRestartStrategy,
    LIST({ryujin::EulerPoisson::GaussLawRestartStrategy::no_restart,
          "no restart"},
         {ryujin::EulerPoisson::GaussLawRestartStrategy::full_restart,
          "full restart"},
         {ryujin::EulerPoisson::GaussLawRestartStrategy::correction,
          "correction"},
         {ryujin::EulerPoisson::GaussLawRestartStrategy::static_no_restart,
          "static no restart"},
         {ryujin::EulerPoisson::GaussLawRestartStrategy::static_full_restart,
          "static full restart"}));
#endif
 
namespace ryujin
{
  namespace EulerPoisson
  {
    template <typename HyperbolicDescription, int dim, typename Number>
    class ParabolicModule final : public Interface::ParabolicModule<dim, Number>
    {
    public:
 
      using StageVectors =
          typename Interface::ParabolicModule<dim, Number>::StageVectors;
 
      using StageWeights =
          typename Interface::ParabolicModule<dim, Number>::StageWeights;
 
      using HyperbolicSystem = typename HyperbolicDescription::HyperbolicSystem;
 
      using View = typename HyperbolicSystem::template View<dim, Number>;
 
      using ParabolicSystem = EulerPoisson::ParabolicSystem;
 
      using BlockHostVector =
          dealii::LinearAlgebra::distributed::BlockVector<Number>;
 
      using ScalarNumber = typename View::ScalarNumber;
 
      static constexpr auto problem_dimension = View::problem_dimension;
 
      using state_type = typename View::state_type;
 
 
 
      ParabolicModule(const MPIEnsemble &mpi_ensemle,
                      const OfflineData<dim, Number> &offline_data,
                      const HyperbolicSystem &hyperbolic_system,
                      const ParabolicSystem &parabolic_system,
                      const InitialValues<HyperbolicDescription, dim, Number>
                          &initial_values,
                      const std::string &subsection = "/ParabolicModule");
 
      void prepare() override;
 
 
 
      void
      reinit_state_vector(StateVector<Number> &state_vector) const override;
 
      void prepare_state_vector(StateVector<Number> &state_vector,
                                Number t) const override;
 
      Number backward_euler_step(
          const StateVector<Number> &old_state_vector,
          StageVectors stage_state_vectors,
          StageWeights stage_weights,
          StateVector<Number> &new_state_vector,
          Number t,
          Number tau,
          Number tau_max = std::numeric_limits<Number>::max()) const override;
 
      Number crank_nicolson_step(
          const StateVector<Number> &old_state_vector,
          StateVector<Number> &new_state_vector,
          Number t,
          Number tau,
          Number tau_max = std::numeric_limits<Number>::max()) const override;
 
 
 
      void print_solver_statistics(std::ostream &output) const override;
 
    private:
 
 
      GaussLawRestartStrategy gauss_law_restart_strategy_;
 
      unsigned int gmg_max_iter_;
      double gmg_smoother_range_;
      double gmg_smoother_max_eig_;
      unsigned int gmg_smoother_degree_;
      unsigned int gmg_smoother_n_cg_iter_;
      unsigned int gmg_min_level_;
 
      Number tolerance_;
      bool tolerance_linfty_norm_;
 
 
 
      // FIXME: refactor
      static constexpr unsigned int order_fe = 1;
      static constexpr unsigned int order_quad = 2;
 
      const MPIEnsemble &mpi_ensemble_;
 
      dealii::ObserverPointer<const HyperbolicSystem> hyperbolic_system_;
      dealii::ObserverPointer<const ParabolicSystem> parabolic_system_;
      dealii::ObserverPointer<const ryujin::OfflineData<dim, Number>>
          offline_data_;
      dealii::ObserverPointer<
          const ryujin::InitialValues<HyperbolicDescription, dim, Number>>
          initial_values_;
 
      std::vector<std::unique_ptr<ElectrostaticConfiguration<dim, Number>>>
          electrostatic_configuration_list_;
 
      ElectrostaticConfiguration<dim, Number>
          *selected_electrostatic_configuration_ = nullptr;
 
      mutable unsigned int cycle_;
      mutable double n_iterations_gauss_;
      mutable double n_iterations_step_;
 
      dealii::MatrixFree<dim, Number> matrix_free_;
      dealii::AffineConstraints<Number> affine_constraints_potential_;
 
      LaplaceOperator<dim, Number> laplace_operator_;
      dealii::DiagonalMatrix<ScalarHostVector<Number>> diagonal_preconditioner_;
      MGSmoother<dim, Number> multigrid_preconditioner_;
      UpdateOperator<dim, Number> update_operator_;
 
      mutable ScalarHostVector<Number> potential_rhs_;
      mutable ScalarHostVector<Number> density_;
      mutable ScalarHostVector<Number> background_density_;
      mutable BlockHostVector magnetic_field_;
      mutable BlockHostVector velocity_rhs_;
 
      mutable bool potential_initialized_;
      mutable Number t_background_density_;
      mutable Number t_magnetic_field_;
 
 
      template <int stages>
      void backward_euler_step(const StateVector<Number> &old_state_vector,
                               StageVectors stage_state_vectors,
                               StageWeights stage_weights,
                               StateVector<Number> &new_state_vector,
                               Number t,
                               Number tau) const;
 
      void create_constraints();
 
      void update_background_density(const Number t) const;
 
      void update_magnetic_field(const Number t) const;
 
      void compute_potential(const Number t,
                             StateVector<Number> &state_vector) const;
 
      void
      enforce_magnetic_drift_velocity(StateVector<Number> &state_vector) const;
 
      void step(const StateVector<Number> &old_state_vector,
                StateVector<Number> &new_state_vector,
                Number t,
                Number tau,
                const bool crank_nicolson_extrapolation) const;
 
    };
  } // namespace EulerPoisson
} // namespace ryujin