parabolic_solver.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2023 - 2024 by the ryujin authors
//
 
#pragma once
 
#include "hyperbolic_module.h"
 
#include <compile_time_options.h>
 
#include <convenience_macros.h>
#include <initial_values.h>
#include <mpi_ensemble.h>
#include <offline_data.h>
#include <simd.h>
#include <sparse_matrix_simd.h>
 
#include "parabolic_solver_gmg_operators.h"
 
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/timer.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;
 
    template <typename Description, int dim, typename Number = double>
    class ParabolicSolver final : public dealii::ParameterAcceptor
    {
    public:
 
      using HyperbolicSystem = typename Description::HyperbolicSystem;
 
      using View =
          typename Description::template HyperbolicSystemView<dim, Number>;
 
      using ParabolicSystem = typename Description::ParabolicSystem;
 
      using ScalarNumber = typename View::ScalarNumber;
 
      static constexpr auto problem_dimension = View::problem_dimension;
 
      using state_type = typename View::state_type;
 
      using StateVector = typename View::StateVector;
 
      using ScalarVector = Vectors::ScalarVector<Number>;
 
      using BlockVector = Vectors::BlockVector<Number>;
 
 
 
      ParabolicSolver(
          const MPIEnsemble &mpi_ensemble,
          std::map<std::string, dealii::Timer> &computing_timer,
          const HyperbolicSystem &hyperbolic_system,
          const ParabolicSystem &parabolic_system,
          const OfflineData<dim, Number> &offline_data,
          const InitialValues<Description, dim, Number> &initial_values,
          const std::string &subsection = "ParabolicSolver");
 
      void prepare();
 
 
 
      void backward_euler_step(const StateVector &old_state_vector,
                               const Number old_t,
                               StateVector &new_state_vector,
                               Number tau,
                               const IDViolationStrategy id_violation_strategy,
                               const bool reinitialize_gmg) const;
 
      void crank_nicolson_step(const StateVector &old_state_vector,
                               const Number old_t,
                               StateVector &new_state_vector,
                               Number tau,
                               const IDViolationStrategy id_violation_strategy,
                               const bool reinitialize_gmg) const;
 
      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_;
      ACCESSOR_READ_ONLY(use_gmg_velocity)
 
      bool use_gmg_internal_energy_;
      ACCESSOR_READ_ONLY(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_;
 
 
 
      void step(const StateVector &old_state_vector,
                const Number old_t,
                StateVector &new_state_vector,
                Number tau,
                const IDViolationStrategy id_violation_strategy,
                const bool reinitialize_gmg,
                const bool crank_nicolson_extrapolation) const;
 
 
 
      // FIXME: refactor
      static constexpr unsigned int order_fe = 1;
      static constexpr unsigned int order_quad = 2;
 
      const MPIEnsemble &mpi_ensemble_;
      std::map<std::string, dealii::Timer> &computing_timer_;
 
      dealii::SmartPointer<const HyperbolicSystem> hyperbolic_system_;
      dealii::SmartPointer<const ParabolicSystem> parabolic_system_;
      dealii::SmartPointer<const ryujin::OfflineData<dim, Number>>
          offline_data_;
      dealii::SmartPointer<
          const ryujin::InitialValues<Description, dim, Number>>
          initial_values_;
 
      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 BlockVector velocity_;
      mutable BlockVector velocity_rhs_;
      mutable ScalarVector internal_energy_;
      mutable ScalarVector internal_energy_rhs_;
      mutable ScalarVector 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_;
 
    };
 
  } // namespace NavierStokes
} /* namespace ryujin */