parabolic_module.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2023 - 2026 by the ryujin authors
//
 
#pragma once
 
#include <ryujin/base/compile_time_options.h>
#include <ryujin/base/convenience_macros.h>
#include <ryujin/base/mpi_ensemble.h>
#include <ryujin/base/observer_pointer.h>
#include <ryujin/base/simd.h>
#include <ryujin/discretization/offline_data.h>
#include <ryujin/interface/parabolic_module.h>
#include <ryujin/navier_stokes/parabolic_module_gmg_operators.h>
#include <ryujin/typed/initial_values.h>
 
#include <deal.II/base/mg_level_object.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/lac/sparse_matrix.templates.h>
#include <deal.II/lac/vector.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
 
#include <limits>
 
namespace ryujin
{
  namespace Euler
  {
    struct HyperbolicDescription;
  } // namespace Euler
 
  namespace NavierStokes
  {
    template <int, typename>
    class DiagonalMatrix;
 
    template <int dim, typename Number = double>
    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 HyperbolicDescription = Euler::HyperbolicDescription;
 
      using HyperbolicSystem = Euler::HyperbolicSystem;
 
      using View = Euler::HyperbolicSystemView<dim, Number>;
 
      using ParabolicSystem = NavierStokes::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_ensemble,
                      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
      {
        // do nothing
      }
 
      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:
 
 
      bool use_gmg_velocity_;
      bool use_gmg_internal_energy_;
 
      Number tolerance_;
      bool tolerance_linfty_norm_;
 
      unsigned int gmg_max_iter_vel_;
      unsigned int gmg_max_iter_en_;
      double gmg_smoother_range_vel_;
      double gmg_smoother_range_en_;
      double gmg_smoother_max_eig_vel_;
      double gmg_smoother_max_eig_en_;
      unsigned int gmg_smoother_degree_;
      unsigned int gmg_smoother_n_cg_iter_;
      unsigned int gmg_min_level_;
 
 
 
      // 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_;
 
      mutable unsigned int cycle_;
      mutable double n_iterations_velocity_;
      mutable double n_iterations_internal_energy_;
 
      mutable dealii::MatrixFree<dim, Number> matrix_free_;
 
      mutable BlockHostVector velocity_;
      mutable BlockHostVector velocity_rhs_;
      mutable ScalarHostVector<Number> internal_energy_;
      mutable ScalarHostVector<Number> internal_energy_rhs_;
      mutable ScalarHostVector<Number> density_;
 
      mutable dealii::MGLevelObject<dealii::MatrixFree<dim, float>>
          level_matrix_free_;
      mutable dealii::MGConstrainedDoFs mg_constrained_dofs_;
      mutable dealii::MGLevelObject<
          dealii::LinearAlgebra::distributed::Vector<float>>
          level_density_;
      mutable MGTransferVelocity<dim, float> mg_transfer_velocity_;
      mutable dealii::MGLevelObject<VelocityMatrix<dim, float, Number>>
          level_velocity_matrices_;
      mutable MGTransferEnergy<dim, float> mg_transfer_energy_;
      mutable dealii::MGLevelObject<EnergyMatrix<dim, float, Number>>
          level_energy_matrices_;
 
      mutable dealii::mg::SmootherRelaxation<
          dealii::PreconditionChebyshev<
              VelocityMatrix<dim, float, Number>,
              dealii::LinearAlgebra::distributed::BlockVector<float>,
              DiagonalMatrix<dim, float>>,
          dealii::LinearAlgebra::distributed::BlockVector<float>>
          mg_smoother_velocity_;
 
      mutable dealii::mg::SmootherRelaxation<
          dealii::PreconditionChebyshev<
              EnergyMatrix<dim, float, Number>,
              dealii::LinearAlgebra::distributed::Vector<float>>,
          dealii::LinearAlgebra::distributed::Vector<float>>
          mg_smoother_energy_;
 
 
 
      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 step(const StateVector<Number> &old_state_vector,
                StateVector<Number> &new_state_vector,
                Number t,
                Number tau,
                const bool crank_nicolson_extrapolation) const;
 
    };
 
  } // namespace NavierStokes
} // namespace ryujin