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 computation can be periodically extended to the other boxes; this requires the solution to have a periodic structure.
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):
The second and third step are necessary for distributed meshes 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 have to be omitted.
This call loops over all faces of the container dof_handler on the opposing boundaries with boundary indicator b_id1 and b_id2, respecitvely. If \(\text{vertices}_{1/2}\) are the vertices of \(\text{face}_{1/2}\), it matches 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. Then,
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:
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:=\{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*}
In order to implement periodic boundary conditions only two functions have to be modified:
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 reuslitng information into is here based on the Triangulation.
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\).
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.
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:
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.
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:
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.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/mapping_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
namespace Step45
{
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int degree);
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;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
};
template <int dim>
class BoundaryValues :
public Function<dim>
{
public:
BoundaryValues () :
Function<dim>(dim+1) {}
const unsigned int component = 0) const;
};
template <int dim>
double
const unsigned int component) const
{
(void) component;
Assert (component < this->n_components,
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);
}
template <int dim>
class RightHandSide :
public Function<dim>
{
public:
RightHandSide () :
Function<dim>(dim+1) {}
const unsigned int component = 0) const;
};
template <int dim>
double
const unsigned int component) const
{
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,
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = RightHandSide<dim>::value (p, c);
}
template <class MatrixType, class PreconditionerType>
{
public:
InverseMatrix (const MatrixType &m,
const PreconditionerType &preconditioner,
const MPI_Comm &mpi_communicator);
private:
const MPI_Comm *mpi_communicator;
};
template <class MatrixType, class PreconditionerType>
InverseMatrix<MatrixType,PreconditionerType>::InverseMatrix
(const MatrixType &m,
const PreconditionerType &preconditioner,
const MPI_Comm &mpi_communicator)
:
preconditioner (&preconditioner),
mpi_communicator (&mpi_communicator),
tmp(locally_owned, mpi_communicator)
{}
template <class MatrixType, class PreconditionerType>
void InverseMatrix<MatrixType,PreconditionerType>::vmult
{
tmp = 0.;
cg.solve (*matrix, tmp, src, *preconditioner);
dst = tmp;
}
template <class PreconditionerType>
{
public:
PreconditionerType> &A_inverse,
const MPI_Comm &mpi_communicator);
private:
PreconditionerType> > A_inverse;
};
template <class PreconditionerType>
SchurComplement<PreconditionerType>::
PreconditionerType> &A_inverse,
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
{
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,
dof_handler (triangulation),
mapping(degree+1)
{}
template <int dim>
void StokesProblem<dim>::create_mesh()
{
const double inner_radius = .5;
const double outer_radius = 1.;
center,
inner_radius,
outer_radius,
0,
true);
std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
periodicity_vector;
rotation_matrix[0][1]=1.;
rotation_matrix[1][0]=-1.;
rotation_matrix);
triangulation.add_periodicity(periodicity_vector);
triangulation.refine_global (4-dim);
}
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
std::vector<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
std::vector<types::global_dof_index> dofs_per_block (2);
const unsigned int n_u = dofs_per_block[0],
n_p = dofs_per_block[1];
{
owned_partitioning.clear();
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();
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);
constraints);
dof_handler,
0,
BoundaryValues<dim>(),
constraints,
dof_handler,
1,
BoundaryValues<dim>(),
constraints,
rotation_matrix[0][1]=1.;
rotation_matrix[1][0]=-1.;
std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
periodicity_vector;
const unsigned int direction = 1;
periodicity_vector, offset,
rotation_matrix);
std::vector<unsigned int> first_vector_components;
first_vector_components.push_back(0);
DoFTools::make_periodicity_constraints<DoFHandler<dim> >
first_vector_components);
dof_handler,
0,
BoundaryValues<dim>(),
constraints,
dof_handler,
1,
BoundaryValues<dim>(),
constraints,
}
constraints.close ();
{
(owned_partitioning, owned_partitioning,
relevant_partitioning, mpi_communicator);
(dof_handler, bsp, constraints, false,
bsp.compress();
system_matrix.reinit (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.;
const unsigned int n_q_points = quadrature_formula.
size();
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,
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);
endc = dof_handler.
end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
local_matrix = 0;
local_rhs = 0;
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]
- div_phi_u[i] * phi_p[j]
- phi_p[i] * div_phi_u[j]
+ phi_p[i] * phi_p[j])
}
const unsigned int component_i =
rhs_values[q](component_i) *
}
}
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);
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
pcout << " Computing preconditioner..." << std::endl << std::flush;
}
template <int dim>
void StokesProblem<dim>::solve ()
{
A_preconditioner.
initialize(system_matrix.block(0,0));
A_inverse (system_matrix.block(0,0),
A_preconditioner,
owned_partitioning[0],
mpi_communicator);
mpi_communicator);
{
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>
owned_partitioning[0],
mpi_communicator);
1
e-6*schur_rhs.l2_norm());
preconditioner.
initialize (system_matrix.block(1,1));
m_inverse (system_matrix.block(1,1), preconditioner,
owned_partitioning[1], mpi_communicator);
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
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();
std::ofstream output ("solution-"
+ "."
+ ".vtu");
{
std::vector<std::string> filenames;
filenames.push_back (std::string("solution-") +
"." +
".vtu");
const std::string
pvtu_master_filename = ("solution-" +
".pvtu");
std::ofstream pvtu_master (pvtu_master_filename.c_str());
}
}
template <int dim>
void
StokesProblem<dim>::refine_mesh ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
solution,
estimated_error_per_cell,
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;
}
}
}
int main (int argc, char *argv[])
{
try
{
using namespace Step45;
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;
}
1.8.14