Reference documentation for deal.II version Git 1aa49a030b 2021-12-07 05:12:24 -0500
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
The 'Laplace equation coupled to an external simulation program' code gallery program

This program was contributed by David Schneider <dav.schneider@tum.de>.
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

Laplace equation coupled to an external simulation program

Overview

preCICE allows to couple deal.II to external simulation software, such as OpenFOAM, SU2, or CalculiX. To keep dependencies of this example minimal, we couple deal.II to an external C++ program, which provides a time varying boundary condition. The deal.II code consists mainly of the step-4 tutorial program, where a simple Laplace problem is solved.

Coupling with preCICE is usually carried out along surfaces in order to apply a Dirichlet-Neumann coupling between two domains (volume coupling is also possible). For the sake of simplicity, we couple here an external C++ program in a unidirectional fashion to one side of our quadrilateral domain. The external C++ program generates a parabolic boundary profile with time varying amplitude. The boundary values are then used in the Laplace solver as a Dirichlet boundary condition.

Time discretization

Coupled simulations deal mostly with time-dependent problems. Hence, we make the stationary Laplace problem from step-4 time dependent.

\begin{align*} \frac{\partial u}{\partial t}-\Delta u &= f \qquad\qquad & \text{in}\ \Omega, \\ u &= x^2+y^2 \qquad\qquad & \text{on}\ \partial\Omega_s, \\ u &= g(t) \qquad\qquad & \text{on}\ \partial\Omega_c. \end{align*}

with the fixed Dirichlet boundary , the coupling boundary and the time-dependent coupling data g(t).

The system is consequently discretized by a first-order backward Euler method, resulting in

\begin{align*} \frac{u^{n+1}-u^n}{\Delta t}-\Delta u^{n+1} &= f \qquad\qquad & \text{in}\ \Omega, \end{align*}

at the next time level n+1.

Requirements

deal.II, version 9.2 or greater. Older versions might work as well, but have not been tested.

preCICE, version 2.0 or greater. Have a look at the provided link for an installation guide.

Compiling and running

Similar to the example programs, run

cmake -DDEAL_II_DIR=/path/to/deal.II -Dprecice_DIR=/path/to/precice .

in this directory to configure the problem. The explicit specification of the library locations can be omitted if they are installed globally or via environment variables. You can then switch between debug and release mode by calling either

make debug

or

make release

This command will generate two executables: one for the coupled_laplace_problem and one for the fancy_boundary_condition participant.

./coupled_laplace_problem

executes the coupled_laplace_problem. In order to start the coupled simulation, execute

./fancy_boundary_condition

in the same directory from another terminal window.

preCICE and the deal.II solver create several log files during the simulation. In order to remove all result-associated files and clean-up the simulation directory the Allclean script can be executed.

Results

By default, the results are saved in each time step using vtk files, which can, for example, be visualized with ParaView. You should get something like:

result

Further reading

A complete overview of the preCICE project can be found on the preCICE webpage. The source code of preCICE is hosted on Github. For more coupled deal.II codes, have a look at the deal.II adapter repository. In the preCICE tutorials, a fluid-structure interaction example is given coupling OpenFOAM (fluid dynamics) to deal.II (solid mechanics). The preCICE reference paper

Annotated version of coupled_laplace_problem.cc

The included deal.II header files are the same as in the other example programs:

In addition to the deal.II header files, we include the preCICE API in order to obtain access to preCICE specific functionality

#include <precice/SolverInterface.hpp>
#include <fstream>
#include <iostream>
using namespace dealii;

Configuration parameters

We set up a simple hard-coded struct containing all names we need for external coupling. The struct includes the name of the preCICE configuration file as well as the name of the simulation participant, the name of the coupling mesh and the name of the exchanged data. The last three names you also find in the preCICE configuration file. For real application cases, these names are better handled by a parameter file.

struct CouplingParamters
{
const std::string config_file = "precice-config.xml";
const std::string participant_name = "laplace-solver";
const std::string mesh_name = "dealii-mesh";
const std::string read_data_name = "boundary-data";
};

The Adapter class

The Adapter class handles all functionalities to couple the deal.II solver code to other solvers with preCICE, i.e., data structures are set up and all relevant information is passed to preCICE.

template <int dim, typename ParameterClass>
class Adapter
{
public:
Adapter(const ParameterClass & parameters,
const types::boundary_id dealii_boundary_interface_id);
double
initialize(const DoFHandler<dim> & dof_handler,
std::map<types::global_dof_index, double> &boundary_data,
const MappingQGeneric<dim> & mapping);
double
advance(std::map<types::global_dof_index, double> &boundary_data,
const double computed_timestep_length);

public precCICE solver interface

precice::SolverInterface precice;

Boundary ID of the deal.II triangulation, associated with the coupling interface. The variable is defined in the constructor of this class and intentionally public so that it can be used during the grid generation and system assembly. The only thing, one needs to make sure is that this ID is unique for a particular triangulation.

const unsigned int dealii_boundary_interface_id;
private:

preCICE related initializations These variables are specified in and read from a parameter file, which is in this simple tutorial program the CouplingParameter struct already introduced in the beginning.

const std::string mesh_name;
const std::string read_data_name;

These IDs are filled by preCICE during the initialization. We set a default value of -1 in order to detect potential errors more easily.

int mesh_id;
int read_data_id;
int n_interface_nodes;

DoF IndexSet, containing relevant coupling DoF indices at the coupling boundary

IndexSet coupling_dofs;

Data containers which are passed to preCICE in an appropriate preCICE specific format

std::vector<int> interface_nodes_ids;
std::vector<double> read_data;

The MPI rank and total number of MPI ranks is required by preCICE when the SolverInterface is created. Since this tutorial runs only in serial mode we define the variables manually in this class instead of using the regular MPI interface.

static constexpr int this_mpi_process = 0;
static constexpr int n_mpi_processes = 1;

Function to transform the obtained data from preCICE into an appropriate map for Dirichlet boundary conditions

void
format_precice_to_dealii(
std::map<types::global_dof_index, double> &boundary_data) const;
};

In the constructor of the Adapter class, we set up the preCICE SolverInterface. We need to tell preCICE our name as participant of the simulation and the name of the preCICE configuration file. Both have already been specified in the CouplingParameter class above. Thus, we pass the class directly to the constructor and read out all relevant information. As a second parameter, we need to specify the boundary ID of our triangulation, which is associated with the coupling interface.

template <int dim, typename ParameterClass>
Adapter<dim, ParameterClass>::Adapter(
const ParameterClass & parameters,
const types::boundary_id deal_boundary_interface_id)
: precice(parameters.participant_name,
parameters.config_file,
this_mpi_process,
n_mpi_processes)
, dealii_boundary_interface_id(deal_boundary_interface_id)
, mesh_name(parameters.mesh_name)
, read_data_name(parameters.read_data_name)
{}

This function initializes preCICE (e.g. establishes communication channels and allocates memory) and passes all relevant data to preCICE. For surface coupling, relevant data is in particular the location of the data points at the associated interface(s). The boundary_data is an empty map, which is filled by preCICE, i.e., information of the other participant. Throughout the system assembly, the map can directly be used in order to apply the Dirichlet boundary conditions in the linear system. preCICE returns the maximum admissible time-step size during the initialization.

template <int dim, typename ParameterClass>
double
Adapter<dim, ParameterClass>::initialize(
const DoFHandler<dim> & dof_handler,
std::map<types::global_dof_index, double> &boundary_data,
const MappingQGeneric<dim> & mapping)
{
AssertDimension(dim, precice.getDimensions());

In a first step, we get preCICE specific IDs from preCICE and store them in the respective variables. Later, they are used for data transfer.

mesh_id = precice.getMeshID(mesh_name);
read_data_id = precice.getDataID(read_data_name, mesh_id);

Afterwards, we extract the number of interface nodes and the coupling DoFs at the coupling interface from our deal.II solver via extract_boundary_dofs()

std::set<types::boundary_id> couplingBoundary;
couplingBoundary.insert(dealii_boundary_interface_id);

The ComponentMask() might be important in case we deal with vector valued problems, because vector valued problems have a DoF for each component.

coupling_dofs,
couplingBoundary);

The coupling DoFs are used to set up the boundary_data map. At the end, we associate here each DoF with a respective boundary value.

for (const auto i : coupling_dofs)
boundary_data[i] = 0.0;

Since we deal with a scalar problem, the number of DoFs at the particular interface corresponds to the number of interface nodes.

n_interface_nodes = coupling_dofs.n_elements();
std::cout << "\t Number of coupling nodes: " << n_interface_nodes
<< std::endl;

Now, we need to tell preCICE the coordinates of the interface nodes. Hence, we set up a std::vector to pass the node positions to preCICE. Each node is specified only once.

std::vector<double> interface_nodes_positions;
interface_nodes_positions.reserve(dim * n_interface_nodes);

Set up the appropriate size of the data container needed for data exchange. Here, we deal with a scalar problem, so that only a scalar value is read/written per interface node.

read_data.resize(n_interface_nodes);

The IDs are again filled by preCICE during the initializations.

interface_nodes_ids.resize(n_interface_nodes);

The node location is obtained using map_dofs_to_support_points().

std::map<types::global_dof_index, Point<dim>> support_points;
DoFTools::map_dofs_to_support_points(mapping, dof_handler, support_points);

support_points contains now the coordinates of all DoFs. In the next step, the relevant coordinates are extracted using the IndexSet with the extracted coupling_dofs.

for (const auto element : coupling_dofs)
for (int i = 0; i < dim; ++i)
interface_nodes_positions.push_back(support_points[element][i]);

Now we have all information to define the coupling mesh and pass the information to preCICE.

precice.setMeshVertices(mesh_id,
n_interface_nodes,
interface_nodes_positions.data(),
interface_nodes_ids.data());

Then, we initialize preCICE internally calling the API function initialize()

const double max_delta_t = precice.initialize();

read first coupling data from preCICE if available (i.e. deal.II is the second participant in a serial coupling scheme)

if (precice.isReadDataAvailable())
{
precice.readBlockScalarData(read_data_id,
n_interface_nodes,
interface_nodes_ids.data(),
read_data.data());

After receiving the coupling data in read_data, we convert it to the std::map boundary_data which is later needed in order to apply Dirichlet boundary conditions

format_precice_to_dealii(boundary_data);
}
return max_delta_t;
}

The function advance() is called in the main time loop after the computation in each time step. Here, coupling data is passed to and obtained from preCICE.

template <int dim, typename ParameterClass>
double
std::map<types::global_dof_index, double> &boundary_data,
const double computed_timestep_length)
{

We specify the computed time-step length and pass it to preCICE. In return, preCICE tells us the maximum admissible time-step size our participant is allowed to compute in order to not exceed the next coupling time step.

const double max_delta_t = precice.advance(computed_timestep_length);

As a next step, we obtain data, i.e. the boundary condition, from another participant. We have already all IDs and just need to convert our obtained data to the deal.II compatible 'boundary map' , which is done in the format_deal_to_precice function.

precice.readBlockScalarData(read_data_id,
n_interface_nodes,
interface_nodes_ids.data(),
read_data.data());
format_precice_to_dealii(boundary_data);
return max_delta_t;
}

This function takes the std::vector obtained by preCICE in read_data and inserts the values to the right position in the boundary map used throughout our deal.II solver for Dirichlet boundary conditions. The function is only used internally in the Adapter class and not called in the solver itself. The order, in which preCICE sorts the data in the read_data vector is exactly the same as the order of the initially passed vertices coordinates.

template <int dim, typename ParameterClass>
void
Adapter<dim, ParameterClass>::format_precice_to_dealii(
std::map<types::global_dof_index, double> &boundary_data) const
{

We already stored the coupling DoF indices in the boundary_data map, so that we can simply iterate over all keys in the map.

auto dof_component = boundary_data.begin();
for (int i = 0; i < n_interface_nodes; ++i)
{
AssertIndexRange(i, read_data.size());
boundary_data[dof_component->first] = read_data[i];
++dof_component;
}
}

The solver class is essentially the same as in step-4. We only extend the stationary problem to a time-dependent problem and introduced the coupling. Comments are added at any point, where the workflow differs from step-4.

template <int dim>
class CoupledLaplaceProblem
{
public:
CoupledLaplaceProblem();
void
run();
private:
void
make_grid();
void
setup_system();
void
assemble_system();
void
solve();
void
output_results() const;
DoFHandler<dim> dof_handler;
MappingQ1<dim> mapping;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;

We allocate all structures required for the preCICE coupling: The map is used to apply Dirichlet boundary conditions and filled in the Adapter class with data from the other participant. The CouplingParameters hold the preCICE configuration as described above. The interface boundary ID is the ID associated to our coupling interface and needs to be specified, when we set up the Adapter class object, because we pass it directly to the Constructor of this class.

std::map<types::global_dof_index, double> boundary_data;
CouplingParamters parameters;
const types::boundary_id interface_boundary_id;
Adapter<dim, CouplingParamters> adapter;

The time-step size delta_t is the acutual time-step size used for all computations. The preCICE time-step size is obtained by preCICE in order to ensure a synchronization at all coupling time steps. The solver time step-size is the desired time-step size of our individual solver. In more sophisticated computations, it might be determined adaptively. The time_step counter is just used for the time-step number.

double delta_t;
double precice_delta_t;
const double solver_delta_t = 0.1;
unsigned int time_step = 0;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double
value(const Point<dim> &p, const unsigned int component = 0) const override;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double
value(const Point<dim> &p, const unsigned int component = 0) const override;
};
template <int dim>
double
RightHandSide<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
double return_value = 0.0;
for (unsigned int i = 0; i < dim; ++i)
return_value += 4.0 * std::pow(p(i), 4.0);
return return_value;
}
template <int dim>
double
BoundaryValues<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
return p.square();
}
template <int dim>
CoupledLaplaceProblem<dim>::CoupledLaplaceProblem()
: fe(1)
, dof_handler(triangulation)
, interface_boundary_id(1)
, adapter(parameters, interface_boundary_id)
{}
template <int dim>
void
CoupledLaplaceProblem<dim>::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(4);
for (const auto &cell : triangulation.active_cell_iterators())
for (const auto face : cell->face_iterators())
{

We choose the boundary in positive x direction for the interface coupling.

if (face->at_boundary() && (face->center()[0] == 1))
face->set_boundary_id(interface_boundary_id);
}
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: " << triangulation.n_cells()
<< std::endl;
}
template <int dim>
void
CoupledLaplaceProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void
CoupledLaplaceProblem<dim>::assemble_system()
{

Reset global structures

system_rhs = 0;
system_matrix = 0;

Update old solution values

old_solution = solution;
QGauss<dim> quadrature_formula(fe.degree + 1);
RightHandSide<dim> right_hand_side;
FEValues<dim> fe_values(fe,
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

The solution values from previous time steps are stored for each quadrature point

std::vector<double> local_values_old_solution(fe_values.n_quadrature_points);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;

Get the local values from the `fe_values' object

fe_values.get_function_values(old_solution, local_values_old_solution);

The system matrix contains additionally a mass matrix due to the time discretization. The RHS has contributions from the old solution values.

for (const unsigned int q_index : fe_values.quadrature_point_indices())
for (const unsigned int i : fe_values.dof_indices())
{
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) +=
((fe_values.shape_value(i, q_index) * // phi_i(x_q)
fe_values.shape_value(j, q_index)) + // phi_j(x_q)
(delta_t * // delta t
fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index))) * // grad phi_j(x_q)
fe_values.JxW(q_index); // dx
const auto x_q = fe_values.quadrature_point(q_index);
const auto &local_value = local_values_old_solution[q_index];
cell_rhs(i) += ((delta_t * // delta t
fe_values.shape_value(i, q_index) * // phi_i(x_q)
right_hand_side.value(x_q)) + // f(x_q)
fe_values.shape_value(i, q_index) *
local_value) * // phi_i(x_q)*val
fe_values.JxW(q_index); // dx
}

Copy local to global

cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
{
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
{

At first, we apply the Dirichlet boundary condition from step-4, as usual.

std::map<types::global_dof_index, double> boundary_values;
0,
BoundaryValues<dim>(),
boundary_values);
system_matrix,
solution,
system_rhs);
}
{

Afterwards, we apply the coupling boundary condition. The boundary_data has already been filled by preCICE.

system_matrix,
solution,
system_rhs);
}
}
template <int dim>
void
CoupledLaplaceProblem<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence." << std::endl;
}
template <int dim>
void
CoupledLaplaceProblem<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches(mapping);
std::ofstream output("solution-" + std::to_string(time_step) + ".vtk");
data_out.write_vtk(output);
}
template <int dim>
void
{
std::cout << "Solving problem in " << dim << " space dimensions."
<< std::endl;
make_grid();
setup_system();

After we set up the system, we initialize preCICE using the functionalities of the Adapter. preCICE returns the maximum admissible time-step size, which needs to be compared to our desired solver time-step size.

precice_delta_t = adapter.initialize(dof_handler, boundary_data, mapping);
delta_t = std::min(precice_delta_t, solver_delta_t);

preCICE steers the coupled simulation: isCouplingOngoing is used to synchronize the end of the simulation with the coupling partner

while (adapter.precice.isCouplingOngoing())
{

The time step number is solely used to generate unique output files

++time_step;

In the time loop, we assemble the coupled system and solve it as usual.

assemble_system();
solve();

After we solved the system, we advance the coupling to the next time level. In a bi-directional coupled simulation, we would pass our calculated data to and obtain new data from preCICE. Here, we simply obtain new data from preCICE, so from the other participant. As before, we obtain a maximum time-step size and compare it against the desired solver time-step size.

precice_delta_t = adapter.advance(boundary_data, delta_t);
delta_t = std::min(precice_delta_t, solver_delta_t);

Write an output file if the time step is completed. In case of an implicit coupling, where individual time steps are computed more than once, the function isTimeWindowCompleted prevents unnecessary result writing. For this simple tutorial configuration (explicit coupling), the function returns always true.

if (adapter.precice.isTimeWindowComplete())
output_results();
}
}
int
main()
{
CoupledLaplaceProblem<2> laplace_problem;
laplace_problem.run();
return 0;
}

Annotated version of fancy_boundary_condition.cc

This program does not use any deal.II functionality and depends only on preCICE and the standard libraries.

#include <precice/SolverInterface.hpp>
#include <iostream>
#include <sstream>

The program computes a time-varying parabolic boundary condition, which is passed to preCICE and serves as Dirichlet boundary condition for the other coupling participant.

Function to generate boundary values in each time step

void
define_boundary_values(std::vector<double> &boundary_data,
const double time,
const double end_time)
{

Scale the current time value

const double relative_time = time / end_time;

Define the amplitude. Values run from -0.5 to 0.5

const double amplitude = (relative_time - 0.5);

Specify the actual data we want to pass to the other participant. Here, we choose a parabola with boundary values 2 in order to enforce continuity to adjacent boundaries.

const double n_elements = boundary_data.size();
const double right_zero = boundary_data.size() - 1;
const double left_zero = 0;
const double offset = 2;
for (uint i = 0; i < n_elements; ++i)
boundary_data[i] =
-amplitude * ((i - left_zero) * (i - right_zero)) + offset;
}
int
main()
{
std::cout << "Boundary participant: starting... \n";

Configuration

const std::string configFileName("precice-config.xml");
const std::string solverName("boundary-participant");
const std::string meshName("boundary-mesh");
const std::string dataWriteName("boundary-data");

Adjust to MPI rank and size for parallel computation

const int commRank = 0;
const int commSize = 1;
precice::SolverInterface precice(solverName,
configFileName,
commRank,
commSize);
const int meshID = precice.getMeshID(meshName);
const int dimensions = precice.getDimensions();
const int numberOfVertices = 6;
const int dataID = precice.getDataID(dataWriteName, meshID);

Set up data structures

std::vector<double> writeData(numberOfVertices);
std::vector<int> vertexIDs(numberOfVertices);
std::vector<double> vertices(numberOfVertices * dimensions);

Define a boundary mesh

std::cout << "Boundary participant: defining boundary mesh \n";
const double length = 2;
const double xCoord = 1;
const double deltaY = length / (numberOfVertices - 1);
for (int i = 0; i < numberOfVertices; ++i)
for (int j = 0; j < dimensions; ++j)
{
const unsigned int index = dimensions * i + j;

The x-coordinate is always 1, i.e., the boundary is parallel to the y-axis. The y-coordinate is descending from 1 to -1.

if (j == 0)
vertices[index] = xCoord;
else
vertices[index] = 1 - deltaY * i;
}

Pass the vertices to preCICE

precice.setMeshVertices(meshID,
numberOfVertices,
vertices.data(),
vertexIDs.data());

initialize the Solverinterface

double dt = precice.initialize();

Start time loop

const double end_time = 1;
double time = 0;
while (precice.isCouplingOngoing())
{

Generate new boundary data

define_boundary_values(writeData, time, end_time);
{
std::cout << "Boundary participant: writing coupling data \n";
precice.writeBlockScalarData(dataID,
numberOfVertices,
vertexIDs.data(),
writeData.data());
}
dt = precice.advance(dt);
std::cout << "Boundary participant: advancing in time\n";
time += dt;
}
std::cout << "Boundary participant: closing...\n";
return 0;
}