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 <compile_time_options.h>
 
#include "hyperbolic_module.h"
 
#include "observer_pointer.h"
#include <convenience_macros.h>
#include <initial_values.h>
#include <mpi_ensemble.h>
#include <offline_data.h>
#include <simd.h>
 
#include "parabolic_module_gmg_operators.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>
 
namespace ryujin
{
  namespace NavierStokes
  {
    template <int, typename>
    class DiagonalMatrix;
 
    struct Description;
 
    template <int dim, typename Number = double>
    class ParabolicModule final : public dealii::ParameterAcceptor
    {
    public:
 
      using HyperbolicSystem = Euler::HyperbolicSystem;
 
      using View = Euler::HyperbolicSystemView<dim, Number>;
 
      using ParabolicSystem = NavierStokes::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_ensemble,
          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
      {
        // do nothing
      }
 
      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:
 
 
      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<Description, dim, Number>>
          initial_values_;
 
      mutable IDViolationStrategy id_violation_strategy_;
 
      mutable unsigned int cycle_;
 
      mutable unsigned int n_restarts_;
      mutable unsigned int n_corrections_;
      mutable unsigned int n_warnings_;
      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 internal_energy_;
      mutable ScalarHostVector internal_energy_rhs_;
      mutable ScalarHostVector 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_;
 
 
 
      void step(const StateVector &old_state_vector,
                const Number old_t,
                StateVector &new_state_vector,
                Number tau,
                const bool crank_nicolson_extrapolation) const;
 
    };
 
  } // namespace NavierStokes
} /* namespace ryujin */