This program was contributed by Jie Cheng <chengjiehust@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
Pictures from this code gallery program
Annotated version of Readme.md
Time-dependent Navier-Stokes
General description of the problem
We solve the time-dependent Navier-Stokes equations with implicit-explicit (IMEX) scheme. Here is the equations we want to solve:
\begin{eqnarray*}
{\mathbf{u}}_{,t} - \nu {\nabla}^2\mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \nabla p = \mathbf{f}
\end{eqnarray*}
The idea is as follows: we use backward Euler time for time discretization. The diffusion term is treated implicitly and the convection term is treated explicitly. Let \((u, p)\) denote the velocity and pressure, respectively and \((v, q)\) denote the corresponding test functions, we end up with the following linear system:
\begin{eqnarray*}
m(u^{n+1}, v) + \Delta{t}\cdot a((u^{n+1}, p^{n+1}), (v, q))=m(u^n, v)-\Delta{t}c(u^n;u^n, v)
\end{eqnarray*}
where \(a((u, p), (v, q))\) is the bilinear form of the diffusion term plus the pressure gradient and its transpose (the divergence constraints):
\begin{eqnarray*}
a((u, p), (v, q)) = \int_\Omega \nu\nabla{u}\nabla{v}-p\nabla\cdot v-q\nabla\cdot ud\Omega
\end{eqnarray*}
\(m(u, v)\) is the mass matrix:
\begin{eqnarray*}
m(u, v) = \int_{\Omega} u \cdot v d\Omega
\end{eqnarray*}
and \(c(u;u, v)\) is the convection term:
\begin{eqnarray*}
c(u;u, v) = \int_{\Omega} (u \cdot \nabla u) \cdot v d\Omega
\end{eqnarray*}
Substracting \(m(u^n, v) + \Delta{t}a((u^n, p^n), (v, q))\) from both sides of the equation, we have the incremental form:
\begin{eqnarray*}
m(\Delta{u}, v) + \Delta{t}\cdot a((\Delta{u}, \Delta{p}), (v, q)) = \Delta{t}(-a(u^n, p^n), (q, v)) - \Delta{t}c(u^n;u^n, v)
\end{eqnarray*}
The system we want to solve can be written in matrix form:
\begin{eqnarray*}
\left(
\begin{array}{cc}
A & B^{T} \\
B & 0 \\
\end{array}
\right)
\left(
\begin{array}{c}
U \\
P \\
\end{array}
\right)
=
\left(
\begin{array}{c}
F \\
0 \\
\end{array}
\right)
\end{eqnarray*}
Grad-Div stablization
Similar to step-57, we add \(\gamma B^T M_p^{-1} B\) to the upper left block of the system. This is a term that is consistent, i.e., the corresponding operators applied to the exact solution would be zero. (This is because \(\gamma B^T M_p^{-1} B\) applied to the velocity vector corresponds to the operator \(\gamma\text{grad}\;\text{div}\) applied to the velocity field – which is of course zero because of the incompressibility constraint \(\text{div}\;\mathbf{u}=0\). On the other hand, the term is not zero when applied to a finite element approximation of the exact velocity.) With this, the system becomes:
\begin{eqnarray*}
\left(
\begin{array}{cc}
\tilde{A} & B^{T} \\
B & 0 \\
\end{array}
\right)
\left(
\begin{array}{c}
U \\
P \\
\end{array}
\right)
=
\left(
\begin{array}{c}
F \\
0 \\
\end{array}
\right)
\end{eqnarray*}
where \(\tilde{A} = A + \gamma B^T M_p^{-1} B\).
A detailed explanation of the Grad-Div stablization can be found in [1].
Block preconditioner
The block preconditioner is pretty much the same as in step-22, except for two additional terms, namely the inertial term (mass matrix) and the Grad-Div term.
The block preconditioner can be written as:
\begin{eqnarray*}
P^{-1}
=
\left(
\begin{array}{cc}
{\tilde{A}}^{-1} & 0 \\
{\tilde{S}}^{-1}B{\tilde{A}}^{-1} & -{\tilde{S}}^{-1} \\
\end{array}
\right)
\end{eqnarray*}
where \({\tilde{S}}\) is the Schur complement of \({\tilde{A}}\), which can be decomposed into the Schur complements of the mass matrix, diffusion matrix, and the Grad-Div term:
\begin{eqnarray*}
{\tilde{S}}^{-1}
\approx
{S_{mass}}^{-1} + {S_{diff}}^{-1} + {S_{Grad-Div}}^{-1}
\approx
[B^T (diag M)^{-1} B]^{-1} + \Delta{t}(\nu + \gamma)M_p^{-1}
\end{eqnarray*}
For more information about preconditioning incompressible Navier-Stokes equations, please refer to [1] and [2].
Test case
We test the code with a classical benchmark case, flow past a cylinder, in both 2D and 3D. The geometry setup of the case can be found on this webpage. The video shows the 2D flow when \(Re = 100\), where mesh refinement is periodically performed. To test the parallel scaling, a 3D case with 1009804 degrees of freedom was ran for 10 time steps on different number of (Xeon E5-2560) processors, results are shown in the graph.
Acknowledgements
Thanks go to Wolfgang Bangerth, Timo Heister and Martin Kronbichler for their helpful discussions on my numerical formulation and implementation.
References
[1] Timo Heister. A massively parallel finite element framework with application to incompressible flows. Doctoral dissertation, University of Gottingen, 2011.
[2] M. Kronbichler, A. Diagne and H. Holmgren. A fast massively parallel two-phase flow solver for microfluidic chip simulation, International Journal of High Performance Computing Applications, 2016.
Annotated version of time_dependent_navier_stokes.cc
#include <fstream>
#include <iostream>
#include <sstream>
namespace fluid
{
Create the triangulation
The code to create triangulation is copied from Martin Kronbichler's code with very few modifications.
Helper function
bool compute_in_2d = true)
{
left,
std::vector<unsigned int>({3U, 4U}),
false);
right,
std::vector<unsigned int>({18U, 4U}),
false);
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
const ::Triangulation< dim, spacedim > & tria
Create middle part first as a hyper shell.
++cell)
for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
{
bool is_inner_rim = true;
for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_face; ++v)
{
Point<2> &vertex = cell->face(f)->vertex(v);
{
is_inner_rim = false;
break;
}
}
if (is_inner_rim)
cell->face(f)->set_manifold_id(1);
}
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
cell_iterator begin(const unsigned int level=0) const
void refine_global(const unsigned int times=1)
cell_iterator end() const
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
void reset_all_manifolds()
void hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Then move the vertices to the points where we want them to be to create a slightly asymmetric cube with a hole:
++cell)
for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
{
if (
std::abs(vertex[0] - 0.7) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.6) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.6) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.5) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.5) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.4) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.4) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.3) < 1e-10 &&
else if (
std::abs(vertex[0] - 0.56379) < 1e-4 &&
else if (
std::abs(vertex[0] - 0.56379) < 1e-4 &&
else if (
std::abs(vertex[0] - 0.43621) < 1e-4 &&
else if (
std::abs(vertex[0] - 0.43621) < 1e-4 &&
}
Refine once to create the same level of refinement as in the neighboring domains:
Must copy the triangulation because we cannot merge triangulations with refinement:
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
Left domain is requred in 3d only.
if (compute_in_2d)
{
}
else
{
}
}
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false)
2D flow around cylinder triangulation
{
create_triangulation_2d(
tria);
Set the left boundary (inflow) to 0, the right boundary (outflow) to 1, upper to 2, lower to 3 and the cylindrical surface to 4.
++cell)
{
for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary())
{
if (
std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(1);
}
else if (
std::abs(cell->face(f)->center()[0] - 0.3) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(0);
}
else if (
std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(3);
}
else if (
std::abs(cell->face(f)->center()[1]) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(2);
}
else
{
cell->face(f)->set_all_boundary_ids(4);
}
}
}
}
}
3D flow around cylinder triangulation
{
create_triangulation_2d(tria_2d, false);
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
Set the ids of the boundaries in x direction to 0 and 1; y direction to 2 and 3; z direction to 4 and 5; the cylindrical surface 6.
++cell)
{
for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary())
{
if (
std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(1);
}
else if (
std::abs(cell->face(f)->center()[0]) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(0);
}
else if (
std::abs(cell->face(f)->center()[1] - 0.41) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(3);
}
else if (
std::abs(cell->face(f)->center()[1]) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(2);
}
else if (
std::abs(cell->face(f)->center()[2] - 0.41) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(5);
}
else if (
std::abs(cell->face(f)->center()[2]) < 1e-12)
{
cell->face(f)->set_all_boundary_ids(4);
}
else
{
cell->face(f)->set_all_boundary_ids(6);
}
}
}
}
}
Time stepping
This class is pretty much self-explanatory.
class Time
{
public:
Time(const double time_end,
const double delta_t,
const double output_interval,
const double refinement_interval)
: timestep(0),
time_current(0.0),
time_end(time_end),
delta_t(delta_t),
output_interval(output_interval),
refinement_interval(refinement_interval)
{
}
double current() const { return time_current; }
double end()
const {
return time_end; }
double get_delta_t() const { return delta_t; }
unsigned int get_timestep() const { return timestep; }
bool time_to_output() const;
bool time_to_refine() const;
void increment();
private:
unsigned int timestep;
double time_current;
const double time_end;
const double delta_t;
const double output_interval;
const double refinement_interval;
};
bool Time::time_to_output() const
{
unsigned int delta = static_cast<unsigned int>(output_interval / delta_t);
return (timestep >= delta && timestep % delta == 0);
}
bool Time::time_to_refine() const
{
unsigned int delta = static_cast<unsigned int>(refinement_interval / delta_t);
return (timestep >= delta && timestep % delta == 0);
}
void Time::increment()
{
time_current += delta_t;
++timestep;
}
VectorType::value_type * end(VectorType &V)
Boundary values
Dirichlet boundary conditions for the velocity inlet and walls.
template <int dim>
class BoundaryValues :
public Function<dim>
{
public:
BoundaryValues() :
Function<dim>(dim + 1) {}
const unsigned int component) const;
};
template <int dim>
double BoundaryValues<dim>::value(
const Point<dim> &p,
const unsigned int component) const
{
Assert(component < this->n_components,
ExcIndexRange(component, 0, this->n_components));
double left_boundary = (dim == 2 ? 0.3 : 0.0);
if (component == 0 &&
std::abs(p[0] - left_boundary) < 1
e-10)
{
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
#define Assert(cond, exc)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
For a parabolic velocity profile, \(U_\mathrm{avg} = 2/3
U_\mathrm{max}\) in 2D, and \(U_\mathrm{avg} = 4/9 U_\mathrm{max}\) in 3D. If \(\nu = 0.001\), \(D = 0.1\), then \(Re = 100 U_\mathrm{avg}\).
double Uavg = 1.0;
double Umax = (dim == 2 ? 3 * Uavg / 2 : 9 * Uavg / 4);
double value = 4 * Umax * p[1] * (0.41 - p[1]) / (0.41 * 0.41);
if (dim == 3)
{
value *= 4 * p[2] * (0.41 - p[2]) / (0.41 * 0.41);
}
return value;
}
return 0;
}
template <int dim>
void BoundaryValues<dim>::vector_value(
const Point<dim> &p,
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryValues<dim>::value(p, c);
}
Block preconditioner
The block Schur preconditioner can be written as the product of three matrices: \(
P^{-1} = \begin{pmatrix} \tilde{A}^{-1} & 0\\ 0 & I\end{pmatrix}
\begin{pmatrix} I & -B^T\\ 0 & I\end{pmatrix}
\begin{pmatrix} I & 0\\ 0 & \tilde{S}^{-1}\end{pmatrix}
\) \(\tilde{A}\) is symmetric since the convection term is eliminated from the LHS. \(\tilde{S}^{-1}\) is the inverse of the Schur complement of \(\tilde{A}\), which consists of a reaction term, a diffusion term, a Grad-Div term and a convection term. In practice, the convection contribution is ignored, namely \(\tilde{S}^{-1} = -(\nu + \gamma)M_p^{-1} -
\frac{1}{\Delta{t}}{[B(diag(M_u))^{-1}B^T]}^{-1}\) where \(M_p\) is the pressure mass, and \({[B(diag(M_u))^{-1}B^T]}\) is an approximation to the Schur complement of (velocity) mass matrix \(BM_u^{-1}B^T\).
Same as the tutorials, we define a vmult operation for the block preconditioner instead of write it as a matrix. It can be seen from the above definition, the result of the vmult operation of the block preconditioner can be obtained from the results of the vmult operations of \(M_u^{-1}\), \(M_p^{-1}\), \(\tilde{A}^{-1}\), which can be transformed into solving three symmetric linear systems.
{
public:
BlockSchurPreconditioner(
double gamma,
double viscosity,
double dt,
const std::vector<IndexSet> &owned_partitioning,
private:
const double gamma;
const double viscosity;
const double dt;
system_matrix;
As discussed, \({[B(diag(M_u))^{-1}B^T]}\) and its inverse need to be computed. We can either explicitly compute it out as a matrix, or define it as a class with a vmult operation. The second approach saves some computation to construct the matrix, but leads to slow convergence in CG solver because it is impossible to apply a preconditioner. We go with the first route.
BlockSchurPreconditioner::BlockSchurPreconditioner
Input parameters and system matrix, mass matrix as well as the mass schur matrix are needed in the preconditioner. In addition, we pass the partitioning information into this class because we need to create some temporary block vectors inside.
BlockSchurPreconditioner::BlockSchurPreconditioner(
double gamma,
double viscosity,
double dt,
const std::vector<IndexSet> &owned_partitioning,
: timer(timer),
gamma(gamma),
viscosity(viscosity),
dt(dt),
system_matrix(&system),
mass_matrix(&mass),
mass_schur(&schur)
{
The schur complemete of mass matrix is actually being computed here.
tmp1.
reinit(owned_partitioning, mass_matrix->get_mpi_communicator());
tmp2.
reinit(owned_partitioning, mass_matrix->get_mpi_communicator());
tmp1 = 1;
tmp2 = 0;
void reinit(const unsigned int n_blocks, const MPI_Comm &communicator, const size_type block_size, const size_type locally_owned_size, const bool omit_zeroing_entries=false)
Jacobi preconditioner of matrix A is by definition \({diag(A)}^{-1}\), this is exactly what we want to compute.
system_matrix->block(1, 0).mmult(
mass_schur->block(1, 1), system_matrix->block(0, 1), tmp2.
block(0));
}
BlockType & block(const unsigned int i)
BlockSchurPreconditioner::vmult
The vmult operation strictly follows the definition of BlockSchurPreconditioner introduced above. Conceptually it computes \(u = P^{-1}v\).
void BlockSchurPreconditioner::vmult(
{
Temporary vectors
This block computes \(u_1 = \tilde{S}^{-1} v_1\), where CG solvers are used for \(M_p^{-1}\) and \(S_m^{-1}\).
{
mass_schur->get_mpi_communicator());
real_type l2_norm() const
\(-(\nu + \gamma)M_p^{-1}v_1\)
Mp_preconditioner.
initialize(mass_matrix->block(1, 1));
cg_mp.solve(
mass_matrix->block(1, 1), tmp, src.
block(1), Mp_preconditioner);
tmp *= -(viscosity + gamma);
}
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
\(-\frac{1}{dt}S_m^{-1}v_1\)
{
mass_schur->get_mpi_communicator());
PreconditionBlockJacobi works find on Sm if we do not refine the mesh. Because after refine_mesh is called, zero entries will be created on the diagonal (not sure why), which prevents PreconditionBlockJacobi from being used.
Sm_preconditioner.
initialize(mass_schur->block(1, 1));
cg_sm.solve(
mass_schur->block(1, 1), dst.
block(1), src.
block(1), Sm_preconditioner);
}
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Adding up these two, we get \(\tilde{S}^{-1}v_1\).
Compute \(v_0 - B^T\tilde{S}^{-1}v_1\) based on \(u_1\).
system_matrix->block(0, 1).vmult(utmp, dst.
block(1));
utmp *= -1.0;
Finally, compute the product of \(\tilde{A}^{-1}\) and utmp using another CG solver.
{
mass_schur->get_mpi_communicator());
We do not use any preconditioner for this block, which is of course slow, only because the performance of the only two preconditioners available PreconditionBlockJacobi and PreconditionBoomerAMG are even worse than none.
A_preconditioner.
initialize(system_matrix->block(0, 0));
cg_a.solve(
system_matrix->block(0, 0), dst.
block(0), utmp, A_preconditioner);
}
}
The incompressible Navier-Stokes solver
Parallel incompressible Navier Stokes equation solver using implicit-explicit time scheme. This program is built upon dealii tutorials step-57, step-40, step-22, and step-20. The system equation is written in the incremental form, and we treat the convection term explicitly. Therefore the system equation is linear and symmetric, which does not need to be solved with Newton's iteration. The system is further stablized and preconditioned with Grad-Div method, where GMRES solver is used as the outer solver.
template <int dim>
class InsIMEX
{
public:
void run();
private:
void setup_dofs();
void make_constraints();
void initialize_system();
void assemble(bool use_nonzero_constraints, bool assemble_system);
std::pair<unsigned int, double> solve(bool use_nonzero_constraints,
bool assemble_system);
void refine_mesh(const unsigned int, const unsigned int);
void output_results(const unsigned int) const;
double viscosity;
double gamma;
const unsigned int degree;
std::vector<types::global_dof_index> dofs_per_block;
QGauss<dim - 1> face_quad_formula;
void print_summary() const
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
System matrix to be solved
Mass matrix is a block matrix which includes both velocity mass matrix and pressure mass matrix.
The schur complement of mass matrix is not a block matrix. However, because we want to reuse the partition we created for the system matrix, it is defined as a block matrix where only one block is actually used.
The latest known solution.
The increment at a certain time step.
System RHS
The IndexSets of owned velocity and pressure respectively.
std::vector<IndexSet> owned_partitioning;
The IndexSets of relevant velocity and pressure respectively.
std::vector<IndexSet> relevant_partitioning;
The IndexSet of all relevant dofs.
The BlockSchurPreconditioner for the entire system.
std::shared_ptr<BlockSchurPreconditioner> preconditioner;
Time time;
};
InsIMEX::InsIMEX
template <int dim>
: viscosity(0.001),
gamma(0.1),
degree(1),
fe(
FE_Q<dim>(degree + 1), dim,
FE_Q<dim>(degree), 1),
volume_quad_formula(degree + 2),
face_quad_formula(degree + 2),
mpi_communicator(MPI_COMM_WORLD),
pcout(
std::cout,
Utilities::MPI::this_mpi_process(mpi_communicator) == 0),
time(1e0, 1e-3, 1e-2, 1e-2),
timer(
{
}
InsIMEX::setup_dofs
template <int dim>
void InsIMEX<dim>::setup_dofs()
{
The first step is to associate DoFs with a given mesh.
dof_handler.distribute_dofs(fe);
We renumber the components to have all velocity DoFs come before the pressure DoFs to be able to split the solution vector in two blocks which are separately accessed in the block preconditioner.
std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Partitioning.
unsigned int dof_u = dofs_per_block[0];
unsigned int dof_p = dofs_per_block[1];
owned_partitioning.resize(2);
owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, dof_u);
owned_partitioning[1] =
dof_handler.locally_owned_dofs().get_view(dof_u, dof_u + dof_p);
relevant_partitioning.resize(2);
relevant_partitioning[0] = locally_relevant_dofs.get_view(0, dof_u);
relevant_partitioning[1] =
locally_relevant_dofs.get_view(dof_u, dof_u + dof_p);
pcout << " Number of active fluid cells: "
<< " Number of degrees of freedom: " << dof_handler.n_dofs() << " ("
<< dof_u << '+' << dof_p << ')' << std::endl;
}
InsIMEX::make_constraints
template <int dim>
void InsIMEX<dim>::make_constraints()
{
Because the equation is written in incremental form, two constraints are needed: nonzero constraint and zero constraint.
nonzero_constraints.clear();
zero_constraints.clear();
nonzero_constraints.reinit(locally_relevant_dofs);
zero_constraints.reinit(locally_relevant_dofs);
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Apply Dirichlet boundary conditions on all boundaries except for the outlet.
std::vector<unsigned int> dirichlet_bc_ids;
if (dim == 2)
dirichlet_bc_ids = std::vector<unsigned int>{0, 2, 3, 4};
else
dirichlet_bc_ids = std::vector<unsigned int>{0, 2, 3, 4, 5, 6};
for (auto id : dirichlet_bc_ids)
{
id,
BoundaryValues<dim>(),
nonzero_constraints,
fe.component_mask(velocities));
dof_handler,
id,
zero_constraints,
fe.component_mask(velocities));
}
nonzero_constraints.close();
zero_constraints.close();
}
InsIMEX::initialize_system
template <int dim>
void InsIMEX<dim>::initialize_system()
{
preconditioner.reset();
system_matrix.clear();
mass_matrix.clear();
mass_schur.clear();
sparsity_pattern.copy_from(dsp);
dsp,
dof_handler.locally_owned_dofs(),
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(owned_partitioning, dsp, mpi_communicator);
mass_matrix.reinit(owned_partitioning, dsp, mpi_communicator);
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Only the \((1, 1)\) block in the mass schur matrix is used. Compute the sparsity pattern for mass schur in advance. The only nonzero block has the same sparsity pattern as \(BB^T\).
schur_dsp.block(1, 1).compute_mmult_pattern(sparsity_pattern.block(1, 0),
sparsity_pattern.block(0, 1));
mass_schur.reinit(owned_partitioning, schur_dsp, mpi_communicator);
present_solution is ghosted because it is used in the output and mesh refinement functions.
present_solution.reinit(
owned_partitioning, relevant_partitioning, mpi_communicator);
solution_increment is non-ghosted because the linear solver needs a completely distributed vector.
solution_increment.reinit(owned_partitioning, mpi_communicator);
system_rhs is non-ghosted because it is only used in the linear solver and residual evaluation.
system_rhs.reinit(owned_partitioning, mpi_communicator);
}
InsIMEX::assemble
Assemble the system matrix, mass matrix, and the RHS. It can be used to assemble the entire system or only the RHS. An additional option is added to determine whether nonzero constraints or zero constraints should be used. Note that we only need to assemble the LHS for twice: once with the nonzero constraint and once for zero constraint. But we must assemble the RHS at every time step.
template <int dim>
void InsIMEX<dim>::assemble(bool use_nonzero_constraints,
bool assemble_system)
{
if (assemble_system)
{
system_matrix = 0;
mass_matrix = 0;
}
system_rhs = 0;
volume_quad_formula,
face_quad_formula,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = volume_quad_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<Tensor<1, dim>> current_velocity_values(n_q_points);
std::vector<Tensor<2, dim>> current_velocity_gradients(n_q_points);
std::vector<double> current_velocity_divergences(n_q_points);
std::vector<double> current_pressure_values(n_q_points);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
for (auto cell = dof_handler.begin_active(); cell != dof_handler.end();
++cell)
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
if (assemble_system)
{
local_matrix = 0;
local_mass_matrix = 0;
}
local_rhs = 0;
fe_values[velocities].get_function_values(present_solution,
current_velocity_values);
fe_values[velocities].get_function_gradients(
present_solution, current_velocity_gradients);
fe_values[velocities].get_function_divergences(
present_solution, current_velocity_divergences);
fe_values[pressure].get_function_values(present_solution,
current_pressure_values);
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Assemble the system matrix and mass matrix simultaneouly. The mass matrix only uses the \((0, 0)\) and \((1, 1)\) blocks.
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
div_phi_u[k] = fe_values[velocities].divergence(k, q);
grad_phi_u[k] = fe_values[velocities].gradient(k, q);
phi_u[k] = fe_values[velocities].value(k, q);
phi_p[k] = fe_values[pressure].value(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
if (assemble_system)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
local_matrix(i, j) +=
(viscosity *
div_phi_u[i] * phi_p[j] -
phi_p[i] * div_phi_u[j] +
gamma * div_phi_u[j] * div_phi_u[i] +
phi_u[i] * phi_u[j] / time.get_delta_t()) *
fe_values.JxW(q);
local_mass_matrix(i, j) +=
(phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j]) *
fe_values.JxW(q);
}
}
local_rhs(i) -=
grad_phi_u[i]) -
current_velocity_divergences[q] * phi_p[i] -
current_pressure_values[q] * div_phi_u[i] +
gamma * current_velocity_divergences[q] * div_phi_u[i] +
current_velocity_gradients[q] *
current_velocity_values[q] * phi_u[i]) *
fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
use_nonzero_constraints ? nonzero_constraints : zero_constraints;
if (assemble_system)
{
local_rhs,
local_dof_indices,
system_matrix,
system_rhs);
local_mass_matrix, local_dof_indices, mass_matrix);
}
else
{
local_rhs, local_dof_indices, system_rhs);
}
}
}
if (assemble_system)
{
}
}
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
long double gamma(const unsigned int n)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
InsIMEX::solve
Solve the linear system using FGMRES solver with block preconditioner. After solving the linear system, the same AffineConstraints object as used in assembly must be used again, to set the constrained value. The second argument is used to determine whether the block preconditioner should be reset or not.
template <int dim>
std::pair<unsigned int, double>
InsIMEX<dim>::solve(bool use_nonzero_constraints, bool assemble_system)
{
if (assemble_system)
{
preconditioner.reset(new BlockSchurPreconditioner(timer,
gamma,
viscosity,
time.get_delta_t(),
owned_partitioning,
system_matrix,
mass_matrix,
mass_schur));
}
system_matrix.m(), 1e-8 * system_rhs.l2_norm(), true);
Because PETScWrappers::SolverGMRES only accepts preconditioner derived from PETScWrappers::PreconditionBase, we use dealii SolverFGMRES.
The solution vector must be non-ghosted
gmres.solve(system_matrix, solution_increment, system_rhs, *preconditioner);
use_nonzero_constraints ? nonzero_constraints : zero_constraints;
return {solver_control.last_step(), solver_control.last_value()};
}
void distribute(VectorType &vec) const
InsIMEX::run
template <int dim>
void InsIMEX<dim>::run()
{
pcout << "Running with PETSc on "
<< " MPI rank(s)..." << std::endl;
setup_dofs();
make_constraints();
initialize_system();
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Time loop.
bool refined = false;
while (time.end() - time.current() > 1e-12)
{
if (time.get_timestep() == 0)
{
output_results(0);
}
time.increment();
std::cout.precision(6);
std::cout.width(12);
pcout << std::string(96, '*') << std::endl
<< "Time step = " << time.get_timestep()
<< ", at t = " << std::scientific << time.current() << std::endl;
Resetting
Only use nonzero constraints at the very first time step
bool apply_nonzero_constraints = (time.get_timestep() == 1);
We have to assemble the LHS for the initial two time steps: once using nonzero_constraints, once using zero_constraints, as well as the steps imediately after mesh refinement.
bool assemble_system = (time.get_timestep() < 3 || refined);
refined = false;
assemble(apply_nonzero_constraints, assemble_system);
auto state = solve(apply_nonzero_constraints, assemble_system);
Note we have to use a non-ghosted vector to do the addition.
tmp.
reinit(owned_partitioning, mpi_communicator);
tmp = present_solution;
tmp += solution_increment;
present_solution = tmp;
pcout << std::scientific << std::left << " GMRES_ITR = " << std::setw(3)
<< state.first << " GMRES_RES = " << state.second << std::endl;
Output
if (time.time_to_output())
{
output_results(time.get_timestep());
}
if (time.time_to_refine())
{
refine_mesh(0, 4);
refined = true;
}
}
}
InsIMEX::output_result
template <int dim>
void InsIMEX<dim>::output_results(const unsigned int output_index) const
{
pcout << "Writing results..." << std::endl;
std::vector<std::string> solution_names(dim, "velocity");
solution_names.push_back("pressure");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
data_component_interpretation.push_back(
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
@ component_is_part_of_vector
vector to be output must be ghosted
solution_names,
data_component_interpretation);
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
Partition
for (unsigned int i = 0; i < subdomain.size(); ++i)
{
}
std::string basename =
std::string filename =
basename +
".vtu";
std::ofstream output(filename);
static std::vector<std::pair<double, std::string>> times_and_names;
{
for (unsigned int i = 0;
++i)
{
times_and_names.push_back(
{time.current(),
}
std::ofstream pvd_output("navierstokes.pvd");
}
}
virtual void build_patches(const unsigned int n_subdivisions=0)
void write_vtu(std::ostream &out) const
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
InsIMEX::refine_mesh
template <int dim>
void InsIMEX<dim>::refine_mesh(const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
pcout << "Refining mesh..." << std::endl;
face_quad_formula,
{},
present_solution,
estimated_error_per_cell,
fe.component_mask(velocity));
{
++cell)
{
cell->clear_refine_flag();
}
}
++cell)
{
cell->clear_coarsen_flag();
}
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, typename InputVector::value_type > * > &neumann_bc, const InputVector &solution, Vector< float > &error, const ComponentMask &component_mask=ComponentMask(), const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void refine_and_coarsen_fixed_fraction(parallel::distributed::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::NormType::L1_norm)
Prepare to transfer
trans(dof_handler);
trans.prepare_for_coarsening_and_refinement(present_solution);
Refine the mesh
Reinitialize the system
setup_dofs();
make_constraints();
initialize_system();
Transfer solution Need a non-ghosted vector for interpolation
tmp = 0;
trans.interpolate(tmp);
present_solution = tmp;
}
}
main function
int main(int argc, char *argv[])
{
try
{
using namespace fluid;
create_triangulation(
tria);
flow.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
1.9.5