Program Listing for File navier_stokes.h¶
↰ Return to documentation for file (include/immersx/physics/navier_stokes.h)
// ---------------------------------------------------------------------
//
// Copyright (C) 2026 by Luca Heltai
//
// This file is part of the ImmersX application, based on
// the deal.II library.
//
// ---------------------------------------------------------------------
#ifndef immersx_navier_stokes_h
#define immersx_navier_stokes_h
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/parsed_convergence_table.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/tria_base.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_description.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <immersx/algebra/linear_algebra.h>
#include <immersx/core/time_parameters.h>
#include <list>
#include <memory>
#include <string>
#include <utility>
#include <variant>
#include <vector>
namespace ImmersX
{
namespace LA
{
using namespace ImmersXLA;
} // namespace LA
template <int dim, int spacedim = dim>
class NavierStokesParameters : public dealii::ParameterAcceptor
{
public:
explicit NavierStokesParameters(
const std::string &subsection = "/Navier-Stokes/");
std::string output_directory = ".";
std::string output_name = "navier_stokes";
TimeIntervalParameters time_parameters;
FixedStepParameters fixed_step_parameters;
IDAParameters ida_parameters;
unsigned int velocity_degree = 2;
unsigned int pressure_degree = 1;
unsigned int initial_refinement = 1;
std::list<dealii::types::boundary_id> dirichlet_ids{0};
std::string name_of_grid = "hyper_cube";
std::string arguments_for_grid = "-1: 1: false";
std::string triangulation_type = "distributed";
double density = 1.0;
double viscosity = 1.0;
bool include_convective_term = true;
std::string analytical_solution_expression =
dim == 2 ? "0; 0; 0" : "0; 0; 0; 0";
mutable dealii::ParsedConvergenceTable convergence_table;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
rhs;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
bc;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
initial_condition;
mutable dealii::ParameterAcceptorProxy<dealii::ReductionControl>
solver_control;
unsigned int inner_solver_max_steps = 200;
double inner_solver_tolerance = 1.e-10;
bool log_solver_iterations = false;
void
set_time(const double time) const;
};
template <int dim, int spacedim = dim>
class NavierStokesSolver : public dealii::EnableObserverPointer
{
static_assert(dim == spacedim,
"NavierStokesSolver currently supports full-dimensional "
"meshes only.");
static_assert(dim == 2 || dim == 3,
"NavierStokesSolver supports only 2D and 3D meshes.");
public:
using VectorType = LA::MPI::Vector;
using BlockVectorType = LA::MPI::BlockVector;
explicit NavierStokesSolver(
const NavierStokesParameters<dim, spacedim> &par);
void
make_grid();
void
setup_fe();
void
setup_system();
void
assemble_system();
void
solve();
void
output_results() const;
void
advance_one_timestep();
void
run();
void
accept_state(const VectorType &velocity,
const VectorType &pressure,
double time,
unsigned int step_number);
dealii::types::global_dof_index
n_dofs() const;
unsigned int
n_time_steps() const;
unsigned int
timestep_number() const;
double
solution_l2_norm() const;
bool
solution_is_finite() const;
double
system_residual_l2_norm() const;
double
divergence_l2_norm() const;
double
current_time() const;
double
time_step() const;
const dealii::parallel::TriangulationBase<dim, spacedim> &
triangulation() const;
const dealii::DoFHandler<dim, spacedim> &
dof_handler() const;
const dealii::FiniteElement<dim, spacedim> &
finite_element() const;
const dealii::Mapping<dim, spacedim> &
mapping() const;
const dealii::AffineConstraints<double> &
constraints() const;
const LA::MPI::BlockSparseMatrix &
system_matrix() const;
const LA::MPI::BlockSparseMatrix &
mass_matrix() const;
const LA::MPI::BlockSparseMatrix &
continuous_operator() const;
const LA::MPI::SparseMatrix &
velocity_mass_matrix() const;
const LA::MPI::SparseMatrix &
pressure_metric_matrix() const;
void
velocity_forcing_at_time(double time, LA::MPI::Vector &destination) const;
const BlockVectorType &
system_rhs() const;
const BlockVectorType &
solution() const;
const BlockVectorType &
previous_solution() const;
const BlockVectorType &
locally_relevant_solution() const;
const dealii::IndexSet &
locally_owned_dofs() const;
const dealii::IndexSet &
locally_relevant_dofs() const;
const std::vector<dealii::IndexSet> &
locally_owned_dofs_by_block() const;
const std::vector<dealii::IndexSet> &
locally_relevant_dofs_by_block() const;
const dealii::FEValuesExtractors::Vector &
velocity_extractor() const;
const dealii::FEValuesExtractors::Scalar &
pressure_extractor() const;
const dealii::ComponentMask &
velocity_component_mask() const;
double
density() const;
double
viscosity() const;
bool
include_convective_term() const;
dealii::types::global_dof_index
velocity_block_size() const;
void
set_solution(const BlockVectorType &new_solution);
private:
using DistributedTriangulation =
dealii::parallel::distributed::Triangulation<dim, spacedim>;
using FullyDistributedTriangulation =
dealii::parallel::fullydistributed::Triangulation<dim, spacedim>;
using TriangulationVariant =
std::variant<DistributedTriangulation, FullyDistributedTriangulation>;
static TriangulationVariant
make_triangulation_storage(MPI_Comm mpi_communicator);
bool
uses_fully_distributed_triangulation() const;
void
update_constraints();
void
interpolate_initial_condition();
void
update_locally_relevant_solution();
void
initialize_time_control();
const NavierStokesParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
dealii::ConditionalOStream pcout;
mutable dealii::TimerOutput computing_timer;
TriangulationVariant triangulation_storage;
dealii::parallel::TriangulationBase<dim, spacedim> *tria;
std::unique_ptr<dealii::FiniteElement<dim, spacedim>> fe;
std::unique_ptr<dealii::Quadrature<dim>> quadrature;
dealii::MappingQ1<dim, spacedim> mapping_storage;
dealii::DoFHandler<dim, spacedim> dh;
dealii::IndexSet owned_dofs;
dealii::IndexSet relevant_dofs;
std::vector<dealii::IndexSet> owned_dofs_by_block;
std::vector<dealii::IndexSet> relevant_dofs_by_block;
std::vector<dealii::types::global_dof_index> dofs_per_block;
dealii::AffineConstraints<double> constraints_storage;
LA::MPI::BlockSparseMatrix system_matrix_storage;
LA::MPI::BlockSparseMatrix mass_matrix_storage;
LA::MPI::BlockSparseMatrix continuous_operator_storage;
BlockVectorType solution_storage;
BlockVectorType previous_solution_storage;
BlockVectorType system_rhs_storage;
BlockVectorType locally_relevant_solution_storage;
LA::MPI::PreconditionAMG velocity_preconditioner;
LA::MPI::PreconditionAMG pressure_preconditioner;
dealii::FEValuesExtractors::Vector velocity;
dealii::FEValuesExtractors::Scalar pressure;
dealii::ComponentMask velocity_mask;
double current_time_storage = 0.0;
double time_step_storage = 0.0;
unsigned int n_time_steps_storage = 0;
unsigned int timestep_number_storage = 0;
mutable std::vector<std::pair<double, std::string>> times_and_names;
mutable unsigned int output_cycle = 0;
};
} // namespace ImmersX
#endif // immersx_navier_stokes_h