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 <compile_time_options.h>
 
#include "electrostatic_configuration_library.h"
#include "laplace_operator.h"
 
#include <hyperbolic_module.h>
#include <initial_values.h>
#include <mpi_ensemble.h>
#include <offline_data.h>
 
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/matrix_free/matrix_free.h>
 
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 Description, int dim, typename Number>
    class ParabolicModule final : public dealii::ParameterAcceptor
    {
    public:
 
      using HyperbolicSystem = typename Description::HyperbolicSystem;
 
      using View = typename HyperbolicSystem::template View<dim, Number>;
 
      using ParabolicSystem = typename Description::ParabolicSystem;
 
      using StateVector = typename View::StateVector;
 
      using ScalarHostVector = Vectors::ScalarHostVector<Number>;
 
      using BlockHostVector = Vectors::BlockHostVector<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<Description, dim, Number> &initial_values,
          const std::string &subsection = "/ParabolicModule");
 
      void prepare();
 
 
 
      void reinit_state_vector(StateVector &state_vector) const;
 
      void prepare_state_vector(StateVector &state_vector, Number t) const;
 
      template <int stages>
      void
      backward_euler_step(const StateVector &old_state_vector,
                          const Number old_t,
                          std::array<std::reference_wrapper<const StateVector>,
                                     stages> stage_state_vectors,
                          const std::array<Number, stages> stage_weights,
                          StateVector &new_state_vector,
                          Number tau) const;
 
      void crank_nicolson_step(const StateVector &old_state_vector,
                               const Number old_t,
                               StateVector &new_state_vector,
                               Number tau) const;
 
      void set_id_violation_strategy(const IDViolationStrategy &strategy) const
      {
        id_violation_strategy_ = strategy;
      }
 
 
 
      void print_solver_statistics(std::ostream &output) const;
 
      ACCESSOR_READ_ONLY(n_restarts)
 
      
      ACCESSOR_READ_ONLY(n_corrections)
 
      
      ACCESSOR_READ_ONLY(n_warnings)
 
    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<Description, dim, Number>>
          initial_values_;
 
      ElectrostaticConfigurationLibrary::
          electrostatic_configuration_list_type<dim, Number>
              electrostatic_configuration_list_;
 
      std::shared_ptr<ElectrostaticConfigurationLibrary::
                          ElectrostaticConfiguration<dim, Number>>
          selected_electrostatic_configuration_;
 
      mutable IDViolationStrategy id_violation_strategy_;
 
      mutable unsigned int cycle_;
      mutable double n_iterations_gauss_;
      mutable double n_iterations_step_;
 
      mutable unsigned int n_restarts_;
      mutable unsigned int n_corrections_;
      mutable unsigned int n_warnings_;
 
      dealii::MatrixFree<dim, Number> matrix_free_;
      dealii::AffineConstraints<Number> affine_constraints_potential_;
 
      LaplaceOperator<dim, Number> laplace_operator_;
      dealii::DiagonalMatrix<ScalarHostVector> diagonal_preconditioner_;
      MGSmoother<dim, Number> multigrid_preconditioner_;
      UpdateOperator<dim, Number> update_operator_;
 
      mutable ScalarHostVector potential_rhs_;
      mutable ScalarHostVector density_;
      mutable ScalarHostVector background_density_;
      mutable BlockHostVector magnetic_field_;
      mutable BlockHostVector velocity_rhs_;
 
      mutable bool potential_initialized_;
      mutable Number t_background_density_;
      mutable Number t_magnetic_field_;
 
 
      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 &state_vector) const;
 
      void enforce_magnetic_drift_velocity(StateVector &state_vector) const;
 
      void step(const StateVector &old_state_vector,
                const Number old_t,
                StateVector &new_state_vector,
                Number tau,
                const bool crank_nicolson_extrapolation) const;
 
    };
  } // namespace EulerPoisson
} /* namespace ryujin */