Reference documentation for deal.II version 9.5.0
Loading...
Searching...
No Matches
The step-45 tutorial program

This tutorial depends on step-6.

Table of contents

This program was contributed by Daniel Arndt and Matthias Maier.

# Introduction

In this example we present how to use periodic boundary conditions in deal.II. Periodic boundary conditions are algebraic constraints that typically occur in computations on representative regions of a larger domain that repeat in one or more directions.

An example is the simulation of the electronic structure of photonic crystals, because they have a lattice-like structure and, thus, it often suffices to do the actual computation on only one box of the lattice. To be able to proceed this way one has to assume that the model can be periodically extended to the other boxes; this requires the solution to have a periodic structure.

# Procedure

deal.II provides a number of high level entry points to impose periodic boundary conditions. The general approach to apply periodic boundary conditions consists of three steps (see also the Glossary entry on periodic boundary conditions):

1. Create a mesh
2. Identify those pairs of faces on different parts of the boundary across which the solution should be symmetric, using GridTools::collect_periodic_faces()
3. Add the periodicity information to the mesh using parallel::distributed::Triangulation::add_periodicity()
4. Add periodicity constraints using DoFTools::make_periodicity_constraints()

The second and third step are necessary for parallel meshes using the parallel::distributed::Triangulation class to ensure that cells on opposite sides of the domain but connected by periodic faces are part of the ghost layer if one of them is stored on the local processor. If the Triangulation is not a parallel::distributed::Triangulation, these steps are not necessary.

The first step consists of collecting matching periodic faces and storing them in a std::vector of GridTools::PeriodicFacePair. This is done with the function GridTools::collect_periodic_faces() that can be invoked for example like this:

b_id1,
b_id2,
direction,
matched_pairs,
offset = <default value>,
matrix = <default value>,
first_vector_components = <default value>);
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const unsigned int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator > > &matched_pairs, const Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())

This call loops over all faces of the container dof_handler on the periodic boundaries with boundary indicator b_id1 and b_id2, respectively. (You can assign these boundary indicators by hand after creating the coarse mesh, see Boundary indicator. Alternatively, you can also let many of the functions in namespace GridGenerator do this for if you specify the "colorize" flag; in that case, these functions will assign different boundary indicators to different parts of the boundary, with the details typically spelled out in the documentation of these functions.)

Concretely, if $$\text{vertices}_{1/2}$$ are the vertices of two faces $$\text{face}_{1/2}$$, then the function call above will match pairs of faces (and dofs) such that the difference between $$\text{vertices}_2$$ and $$matrix\cdot \text{vertices}_1+\text{offset}$$ vanishes in every component apart from direction and stores the resulting pairs with associated data in matched_pairs. (See GridTools::orthogonal_equality() for detailed information about the matching process.)

Consider, for example, the colored unit square $$\Omega=[0,1]^2$$ with boundary indicator 0 on the left, 1 on the right, 2 on the bottom and 3 on the top faces. (See the documentation of GridGenerator::hyper_cube() for this convention on how boundary indicators are assigned.) Then,

/*b_id1*/ 0,
/*b_id2*/ 1,
/*direction*/ 0,
matched_pairs);

would yield periodicity constraints such that $$u(0,y)=u(1,y)$$ for all $$y\in[0,1]$$.

If we instead consider the parallelogram given by the convex hull of $$(0,0)$$, $$(1,1)$$, $$(1,2)$$, $$(0,1)$$ we can achieve the constraints $$u(0,y)=u(1,y+1)$$ by specifying an offset:

/*b_id1*/ 0,
/*b_id2*/ 1,
/*direction*/ 0,
matched_pairs,
Tensor<1, 2>(0.,1.));

or

/*b_id1*/ 0,
/*b_id2*/ 1,
/*arbitrary direction*/ 0,
matched_pairs,
Tensor<1, 2>(1.,1.));

Here, again, the assignment of boundary indicators 0 and 1 stems from what GridGenerator::parallelogram() documents.

The resulting matched_pairs can be used in DoFTools::make_periodicity_constraints for populating an AffineConstraints object with periodicity constraints:

DoFTools::make_periodicity_constraints(matched_pairs, constraints);
void make_periodicity_constraints(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, AffineConstraints< number > &constraints, const ComponentMask &component_mask=ComponentMask(), const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >(), const number periodicity_factor=1.)

Apart from this high level interface there are also variants of DoFTools::make_periodicity_constraints available that combine those two steps (see the variants of DofTools::make_periodicity_constraints).

There is also a low level interface to DoFTools::make_periodicity_constraints if more flexibility is needed. The low level variant allows to directly specify two faces that shall be constrained:

using namespace DoFTools;
make_periodicity_constraints(face_1,
face_2,
affine_constraints,
component_mask = <default value>;
face_orientation = <default value>,
face_flip = <default value>,
face_rotation = <default value>,
matrix = <default value>);

Here, we need to specify the orientation of the two faces using face_orientation, face_flip and face_orientation. For a closer description have a look at the documentation of DoFTools::make_periodicity_constraints. The remaining parameters are the same as for the high level interface apart from the self-explaining component_mask and affine_constraints.

# A practical example

In the following, we show how to use the above functions in a more involved example. The task is to enforce rotated periodicity constraints for the velocity component of a Stokes flow.

On a quarter-circle defined by $$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}$$ we are going to solve the Stokes problem

\begin{eqnarray*} -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\ -\textrm{div}\; \textbf{u}&=&0,\\ \textbf{u}|_{\Gamma_1}&=&{\bf 0}, \end{eqnarray*}

where the boundary $$\Gamma_1$$ is defined as $$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}$$. For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.

\begin{align*} u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\ u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1]. \end{align*}

The mesh will be generated by GridGenerator::quarter_hyper_shell(), which also documents how it assigns boundary indicators to its various boundaries if its colorize argument is set to true.

# The commented program

This example program is a slight modification of step-22 running in parallel using Trilinos to demonstrate the usage of periodic boundary conditions in deal.II. We thus omit to discuss the majority of the source code and only comment on the parts that deal with periodicity constraints. For the rest have a look at step-22 and the full source code at the bottom.

In order to implement periodic boundary conditions only two functions have to be modified:

• StokesProblem<dim>::setup_dofs(): To populate an AffineConstraints object with periodicity constraints
• StokesProblem<dim>::create_mesh(): To supply a distributed triangulation with periodicity information.

The rest of the program is identical to step-22, so let us skip this part and only show these two functions in the following. (The full program can be found in the "Plain program" section below, though.)

### Setting up periodicity constraints on distributed triangulations

template <int dim>
void StokesProblem<dim>::create_mesh()
{
const double inner_radius = .5;
const double outer_radius = 1.;

triangulation, center, inner_radius, outer_radius, 0, true);

Definition point.h:112
Point< 3 > center
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Before we can prescribe periodicity constraints, we need to ensure that cells on opposite sides of the domain but connected by periodic faces are part of the ghost layer if one of them is stored on the local processor. At this point we need to think about how we want to prescribe periodicity. The vertices $$\text{vertices}_2$$ of a face on the left boundary should be matched to the vertices $$\text{vertices}_1$$ of a face on the lower boundary given by $$\text{vertices}_2=R\cdot \text{vertices}_1+b$$ where the rotation matrix $$R$$ and the offset $$b$$ are given by

\begin{align*} R=\begin{pmatrix} 0&1\\-1&0 \end{pmatrix}, \quad b=\begin{pmatrix}0&0\end{pmatrix}. \end{align*}

The data structure we are saving the resulting information into is here based on the Triangulation.

periodicity_vector;

FullMatrix<double> rotation_matrix(dim);
rotation_matrix[0][1] = 1.;
rotation_matrix[1][0] = -1.;

2,
3,
1,
periodicity_vector,
rotation_matrix);

typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:270

Now telling the triangulation about the desired periodicity is particularly easy by just calling parallel::distributed::Triangulation::add_periodicity.

triangulation.add_periodicity(periodicity_vector);

triangulation.refine_global(4 - dim);
}

template <int dim>
void StokesProblem<dim>::setup_dofs()
{
dof_handler.distribute_dofs(fe);

std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise(dof_handler, block_component);

const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];

{
owned_partitioning.clear();
IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u + n_p));

relevant_partitioning.clear();
const IndexSet locally_relevant_dofs =
relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
relevant_partitioning.push_back(
locally_relevant_dofs.get_view(n_u, n_u + n_p));

constraints.clear();
constraints.reinit(locally_relevant_dofs);

const FEValuesExtractors::Vector velocities(0);

DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dof_handler,
0,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
dof_handler,
1,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));

void clear()
Definition index_set.h:1637
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())

After we provided the mesh with the necessary information for the periodicity constraints, we are now able to actual create them. For describing the matching we are using the same approach as before, i.e., the $$\text{vertices}_2$$ of a face on the left boundary should be matched to the vertices $$\text{vertices}_1$$ of a face on the lower boundary given by $$\text{vertices}_2=R\cdot \text{vertices}_1+b$$ where the rotation matrix $$R$$ and the offset $$b$$ are given by

\begin{align*} R=\begin{pmatrix} 0&1\\-1&0 \end{pmatrix}, \quad b=\begin{pmatrix}0&0\end{pmatrix}. \end{align*}

These two objects not only describe how faces should be matched but also in which sense the solution should be transformed from $$\text{face}_2$$ to $$\text{face}_1$$.

FullMatrix<double> rotation_matrix(dim);
rotation_matrix[0][1] = 1.;
rotation_matrix[1][0] = -1.;

Tensor<1, dim> offset;

For setting up the constraints, we first store the periodicity information in an auxiliary object of type std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> . The periodic boundaries have the boundary indicators 2 (x=0) and 3 (y=0). All the other parameters we have set up before. In this case the direction does not matter. Due to $$\text{vertices}_2=R\cdot \text{vertices}_1+b$$ this is exactly what we want.

std::vector<
periodicity_vector;

const unsigned int direction = 1;

2,
3,
direction,
periodicity_vector,
offset,
rotation_matrix);

Next, we need to provide information on which vector valued components of the solution should be rotated. Since we choose here to just constraint the velocity and this starts at the first component of the solution vector, we simply insert a 0:

std::vector<unsigned int> first_vector_components;
first_vector_components.push_back(0);

After setting up all the information in periodicity_vector all we have to do is to tell make_periodicity_constraints to create the desired constraints.

DoFTools::make_periodicity_constraints<dim, dim>(periodicity_vector,
constraints,
fe.component_mask(
velocities),
first_vector_components);
}

constraints.close();

{
TrilinosWrappers::BlockSparsityPattern bsp(owned_partitioning,
owned_partitioning,
relevant_partitioning,
mpi_communicator);

Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if (!((c == dim) && (d == dim)))
coupling[c][d] = DoFTools::always;
else
coupling[c][d] = DoFTools::none;

coupling,
bsp,
constraints,
false,
mpi_communicator));

bsp.compress();

system_matrix.reinit(bsp);
}

{
owned_partitioning,
owned_partitioning,
relevant_partitioning,
mpi_communicator);

Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if ((c == dim) && (d == dim))
preconditioner_coupling[c][d] = DoFTools::always;
else
preconditioner_coupling[c][d] = DoFTools::none;

preconditioner_coupling,
preconditioner_bsp,
constraints,
false,
mpi_communicator));

preconditioner_bsp.compress();

preconditioner_matrix.reinit(preconditioner_bsp);
}

system_rhs.reinit(owned_partitioning, mpi_communicator);
solution.reinit(owned_partitioning,
relevant_partitioning,
mpi_communicator);
}

void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161

The rest of the program is then again identical to step-22. We will omit it here now, but as before, you can find these parts in the "Plain program" section below.

# Results

The created output is not very surprising. We simply see that the solution is periodic with respect to the left and lower boundary:

Without the periodicity constraints we would have ended up with the following solution:

# The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2008 - 2023 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Daniel Arndt, Matthias Maier, 2015
*
* Based on @ref step_22 "step-22" by Wolfgang Bangerth and Martin Kronbichler
*/
namespace Step45
{
using namespace dealii;
template <int dim>
class StokesProblem
{
public:
StokesProblem(const unsigned int degree);
void run();
private:
void create_mesh();
void setup_dofs();
void assemble_system();
void solve();
void output_results(const unsigned int refinement_cycle) const;
void refine_mesh();
const unsigned int degree;
MPI_Comm mpi_communicator;
DoFHandler<dim> dof_handler;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
TrilinosWrappers::BlockSparseMatrix preconditioner_matrix;
MappingQ<dim> mapping;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues()
: Function<dim>(dim + 1)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> & /*p*/,
const unsigned int component) const
{
(void)component;
Assert(component < this->n_components,
ExcIndexRange(component, 0, this->n_components));
return 0;
}
template <int dim>
void BoundaryValues<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryValues<dim>::value(p, c);
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide()
: Function<dim>(dim + 1)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
};
template <int dim>
double RightHandSide<dim>::value(const Point<dim> & p,
const unsigned int component) const
{
const Point<dim> center(0.75, 0.1);
const double r = (p - center).norm();
if (component == 0)
return std::exp(-100. * r * r);
return 0;
}
template <int dim>
void RightHandSide<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = RightHandSide<dim>::value(p, c);
}
template <class MatrixType, class PreconditionerType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix(const MatrixType & m,
const PreconditionerType &preconditioner,
const IndexSet & locally_owned,
const MPI_Comm mpi_communicator);
const TrilinosWrappers::MPI::Vector &src) const;
private:
};
template <class MatrixType, class PreconditionerType>
InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
const MatrixType & m,
const PreconditionerType &preconditioner,
const IndexSet & locally_owned,
const MPI_Comm mpi_communicator)
: matrix(&m)
, preconditioner(&preconditioner)
, tmp(locally_owned, mpi_communicator)
{}
template <class MatrixType, class PreconditionerType>
void InverseMatrix<MatrixType, PreconditionerType>::vmult(
{
SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
TrilinosWrappers::SolverCG cg(solver_control,
tmp = 0.;
cg.solve(*matrix, tmp, src, *preconditioner);
dst = tmp;
}
template <class PreconditionerType>
class SchurComplement : public TrilinosWrappers::SparseMatrix
{
public:
SchurComplement(const TrilinosWrappers::BlockSparseMatrix &system_matrix,
const InverseMatrix<TrilinosWrappers::SparseMatrix,
PreconditionerType> & A_inverse,
const IndexSet & owned_pres,
const MPI_Comm mpi_communicator);
const TrilinosWrappers::MPI::Vector &src) const;
private:
const SmartPointer<
const InverseMatrix<TrilinosWrappers::SparseMatrix, PreconditionerType>>
A_inverse;
mutable TrilinosWrappers::MPI::Vector tmp1, tmp2;
};
template <class PreconditionerType>
SchurComplement<PreconditionerType>::SchurComplement(
const InverseMatrix<TrilinosWrappers::SparseMatrix, PreconditionerType>
& A_inverse,
const IndexSet &owned_vel,
const MPI_Comm mpi_communicator)
: system_matrix(&system_matrix)
, A_inverse(&A_inverse)
, tmp1(owned_vel, mpi_communicator)
, tmp2(tmp1)
{}
template <class PreconditionerType>
void SchurComplement<PreconditionerType>::vmult(
{
system_matrix->block(0, 1).vmult(tmp1, src);
A_inverse->vmult(tmp2, tmp1);
system_matrix->block(1, 0).vmult(dst, tmp2);
}
template <int dim>
StokesProblem<dim>::StokesProblem(const unsigned int degree)
: degree(degree)
, mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator)
, fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1)
, dof_handler(triangulation)
, pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
, mapping(degree + 1)
{}
template <int dim>
void StokesProblem<dim>::create_mesh()
{
const double inner_radius = .5;
const double outer_radius = 1.;
triangulation, center, inner_radius, outer_radius, 0, true);
periodicity_vector;
FullMatrix<double> rotation_matrix(dim);
rotation_matrix[0][1] = 1.;
rotation_matrix[1][0] = -1.;
2,
3,
1,
periodicity_vector,
rotation_matrix);
triangulation.add_periodicity(periodicity_vector);
triangulation.refine_global(4 - dim);
}
template <int dim>
void StokesProblem<dim>::setup_dofs()
{
dof_handler.distribute_dofs(fe);
std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise(dof_handler, block_component);
const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
{
owned_partitioning.clear();
IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u + n_p));
relevant_partitioning.clear();
const IndexSet locally_relevant_dofs =
relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
relevant_partitioning.push_back(
locally_relevant_dofs.get_view(n_u, n_u + n_p));
constraints.clear();
constraints.reinit(locally_relevant_dofs);
const FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dof_handler,
0,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
dof_handler,
1,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
FullMatrix<double> rotation_matrix(dim);
rotation_matrix[0][1] = 1.;
rotation_matrix[1][0] = -1.;
std::vector<
periodicity_vector;
const unsigned int direction = 1;
2,
3,
direction,
periodicity_vector,
offset,
rotation_matrix);
std::vector<unsigned int> first_vector_components;
first_vector_components.push_back(0);
DoFTools::make_periodicity_constraints<dim, dim>(periodicity_vector,
constraints,
fe.component_mask(
velocities),
first_vector_components);
}
constraints.close();
{
owned_partitioning,
relevant_partitioning,
mpi_communicator);
Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if (!((c == dim) && (d == dim)))
coupling[c][d] = DoFTools::always;
else
coupling[c][d] = DoFTools::none;
coupling,
bsp,
constraints,
false,
mpi_communicator));
bsp.compress();
system_matrix.reinit(bsp);
}
{
owned_partitioning,
owned_partitioning,
relevant_partitioning,
mpi_communicator);
Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if ((c == dim) && (d == dim))
preconditioner_coupling[c][d] = DoFTools::always;
else
preconditioner_coupling[c][d] = DoFTools::none;
preconditioner_coupling,
preconditioner_bsp,
constraints,
false,
mpi_communicator));
preconditioner_bsp.compress();
preconditioner_matrix.reinit(preconditioner_bsp);
}
system_rhs.reinit(owned_partitioning, mpi_communicator);
solution.reinit(owned_partitioning,
relevant_partitioning,
mpi_communicator);
}
template <int dim>
void StokesProblem<dim>::assemble_system()
{
system_matrix = 0.;
system_rhs = 0.;
preconditioner_matrix = 0.;
QGauss<dim> quadrature_formula(degree + 2);
FEValues<dim> fe_values(mapping,
fe,
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> local_preconditioner_matrix(dofs_per_cell,
dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
local_matrix = 0;
local_preconditioner_matrix = 0;
local_rhs = 0;
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
symgrad_phi_u[k] =
fe_values[velocities].symmetric_gradient(k, q);
div_phi_u[k] = fe_values[velocities].divergence(k, q);
phi_p[k] = fe_values[pressure].value(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j <= i; ++j)
{
local_matrix(i, j) +=
(symgrad_phi_u[i] * symgrad_phi_u[j] // diffusion
- div_phi_u[i] * phi_p[j] // pressure force
- phi_p[i] * div_phi_u[j]) // divergence
* fe_values.JxW(q);
local_preconditioner_matrix(i, j) +=
(phi_p[i] * phi_p[j]) * fe_values.JxW(q);
}
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) += fe_values.shape_value(i, q)
* rhs_values[q](component_i)
* fe_values.JxW(q);
}
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
{
local_matrix(i, j) = local_matrix(j, i);
local_preconditioner_matrix(i, j) =
local_preconditioner_matrix(j, i);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(local_matrix,
local_rhs,
local_dof_indices,
system_matrix,
system_rhs);
constraints.distribute_local_to_global(local_preconditioner_matrix,
local_dof_indices,
preconditioner_matrix);
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
pcout << " Computing preconditioner..." << std::endl << std::flush;
}
template <int dim>
void StokesProblem<dim>::solve()
{
A_preconditioner.initialize(system_matrix.block(0, 0));
const InverseMatrix<TrilinosWrappers::SparseMatrix,
A_inverse(system_matrix.block(0, 0),
A_preconditioner,
owned_partitioning[0],
mpi_communicator);
TrilinosWrappers::MPI::BlockVector tmp(owned_partitioning,
mpi_communicator);
{
TrilinosWrappers::MPI::Vector schur_rhs(owned_partitioning[1],
mpi_communicator);
A_inverse.vmult(tmp.block(0), system_rhs.block(0));
system_matrix.block(1, 0).vmult(schur_rhs, tmp.block(0));
schur_rhs -= system_rhs.block(1);
SchurComplement<TrilinosWrappers::PreconditionJacobi> schur_complement(
system_matrix, A_inverse, owned_partitioning[0], mpi_communicator);
SolverControl solver_control(solution.block(1).size(),
1e-6 * schur_rhs.l2_norm());
preconditioner.initialize(preconditioner_matrix.block(1, 1));
m_inverse(preconditioner_matrix.block(1, 1),
preconditioner,
owned_partitioning[1],
mpi_communicator);
cg.solve(schur_complement, tmp.block(1), schur_rhs, preconditioner);
constraints.distribute(tmp);
solution.block(1) = tmp.block(1);
}
{
system_matrix.block(0, 1).vmult(tmp.block(0), tmp.block(1));
tmp.block(0) *= -1;
tmp.block(0) += system_rhs.block(0);
A_inverse.vmult(tmp.block(0), tmp.block(0));
constraints.distribute(tmp);
solution.block(0) = tmp.block(0);
}
}
template <int dim>
void
StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
{
std::vector<std::string> solution_names(dim, "velocity");
solution_names.emplace_back("pressure");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
data_component_interpretation.push_back(
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution,
solution_names,
data_component_interpretation);
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches(mapping, degree + 1);
"./", "solution", refinement_cycle, MPI_COMM_WORLD, 2);
}
template <int dim>
void StokesProblem<dim>::refine_mesh()
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
const FEValuesExtractors::Scalar pressure(dim);
dof_handler,
QGauss<dim - 1>(degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell,
fe.component_mask(pressure));
triangulation, estimated_error_per_cell, 0.3, 0.0);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void StokesProblem<dim>::run()
{
create_mesh();
for (unsigned int refinement_cycle = 0; refinement_cycle < 9;
++refinement_cycle)
{
pcout << "Refinement cycle " << refinement_cycle << std::endl;
if (refinement_cycle > 0)
refine_mesh();
setup_dofs();
pcout << " Assembling..." << std::endl << std::flush;
assemble_system();
pcout << " Solving..." << std::flush;
solve();
output_results(refinement_cycle);
pcout << std::endl;
}
}
} // namespace Step45
int main(int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step45;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
StokesProblem<2> flow_problem(1);
flow_problem.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;
}
std::string write_vtu_with_pvtu_record(const std::string &directory, const std::string &filename_without_extension, const unsigned int counter, const MPI_Comm mpi_communicator, const unsigned int n_digits_for_counter=numbers::invalid_unsigned_int, const unsigned int n_groups=0) const
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
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={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1063
Definition fe_q.h:551
IndexSet get_view(const size_type begin, const size_type end) const
Definition index_set.cc:214
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 initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > vmult(VectorType &dst, const VectorType &src) const
#define Assert(cond, exc)
LinearOperator< Range_2, Domain_2, Payload > schur_complement(const LinearOperator< Domain_1, Range_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const LinearOperator< Range_2, Domain_1, Payload > &C, const LinearOperator< Range_2, Domain_2, Payload > &D)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ matrix
Contents is actually a matrix.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)