Reference documentation for deal.II version 9.4.1
\(\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\}}\)
Loading...
Searching...
No Matches
The 'Convection Diffusion Reaction' code gallery program

This program was contributed by David R. Wells wellsd2@rpi.edu.
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):

Annotated version of readme.md

Overview

I started this project with the intent of better understanding adaptive mesh refinement, parallel computing, and CMake. In particular, I started by writing a uniform mesh, single process solver and then ultimately expanded it into solver/cdr.cc. This example program might be useful to look at if you want to see:

A more complex CMake setup, which builds a shared object library and an executable A simple parallel time stepping problem Use of C++11 lambda functions

The other solvers are available here.

Unlike the other tutorial programs, I have split this solver into a number of files in nested directories. In particular, I used the following strategy (more-or-less copied ASPECT):

The common directory, which hosts files common to the four solvers I wrote along the way. Most of the source files in common/source/ are just template specializations*; they compile the template code for specific dimensions and linear algebra (matrix or vector) types. The common/include/deal.II-cdr/ directory contains both templates and plain header files. The solver/ directory contains the actual solver class and strongly resembles a tutorial program. The file solver/cdr.cc just sets up data structures and then calls routines in libdeal.II-cdr-common to populate them and produce output.

Compiling and running

Like the example programs, run

cmake -DDEAL_II_DIR=/path/to/deal.II .
make

in this directory. The solver may be run as

make run

or, for parallelism across 16 processes,

mpirun -np 16 ./cdr

Why use convection-diffusion-reaction?

This equation exhibits very fine boundary layers (usually, from the literature, these layers have width equal to the square root of the diffusion coefficient on internal layers and the diffusion coefficient on boundary layers). A good way to solve it is to use adaptive mesh refinement to refine the mesh only at interior boundary layers. At the same time, this problem is linear (and does not have a pressure term) so it is much simpler to solve than the Navier-Stokes equations with comparable diffusion (Reynolds number).

I use relatively large diffusion values so I can get away without any additional scheme like SUPG.

Recommended Literature

There are many good books and papers on numerical methods for this equation. A good starting point is "Robust Numerical Methods for Singularly Perturbed Problems" by Roos, Stynes, and Tobiska.

Annotated version of common/include/deal.II-cdr/parameters.h

#ifndef dealii__cdr_parameters_h
#define dealii__cdr_parameters_h
#include <deal.II/base/parameter_handler.h>
#include <string>
parameter_handler.h

I prefer to use the ParameterHandler class in a slightly different way than usual: The class Parameters creates, uses, and then destroys a ParameterHandler inside the read_parameter_file method instead of keeping it around. This is nice because now all of the run time parameters are contained in a simple class and it can be copied or passed around very easily.

namespace CDR
{
using namespace dealii;
class Parameters
{
public:
double inner_radius;
double outer_radius;
double diffusion_coefficient;
double reaction_coefficient;
bool time_dependent_forcing;
unsigned int initial_refinement_level;
unsigned int max_refinement_level;
unsigned int fe_order;
double start_time;
double stop_time;
unsigned int n_time_steps;
unsigned int save_interval;
unsigned int patch_level;
void
read_parameter_file(const std::string &file_name);
private:
void
configure_parameter_handler(ParameterHandler &parameter_handler);
};
} // namespace CDR
#endif

Annotated version of common/include/deal.II-cdr/system_matrix.h

#ifndef dealii__cdr_system_matrix_h
#define dealii__cdr_system_matrix_h
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II-cdr/parameters.h>
#include <functional>
affine_constraints.h
dof_handler.h
quadrature_lib.h

One of the goals I had in writing this entry was to split up functions into different compilation units instead of using one large file. This is the header file for a pair of functions (only one of which I ultimately use) which build the system matrix.

namespace CDR
{
using namespace dealii;
template <int dim, typename MatrixType>
void
create_system_matrix(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
MatrixType & system_matrix);
template <int dim, typename MatrixType>
void
create_system_matrix(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
MatrixType & system_matrix);
} // namespace CDR
#endif
Point
Definition: point.h:111
Tensor
Definition: tensor.h:503

Annotated version of common/include/deal.II-cdr/system_matrix.templates.h

#ifndef dealii__cdr_system_matrix_templates_h
#define dealii__cdr_system_matrix_templates_h
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_matrix.h>
#include <functional>
#include <vector>
namespace CDR
{
using namespace dealii;
fe_values.h
fe_q.h

This is the actual implementation of the create_system_matrix function described in the header file. It is similar to the system matrix assembly routine in step-40.

template <int dim, typename UpdateFunction>
void
internal_create_system_matrix(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
UpdateFunction update_system_matrix)
{
auto & fe = dof_handler.get_fe();
const auto dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
FEValues<dim> fe_values(fe,
quad,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
cell_matrix = 0.0;
cell->get_dof_indices(local_indices);
for (unsigned int q = 0; q < quad.size(); ++q)
{
const auto current_convection =
convection_function(fe_values.quadrature_point(q));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const auto convection_contribution =
current_convection * fe_values.shape_grad(j, q);
cell_matrix(i, j) +=
fe_values.JxW(q) *
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_data.h:433
Quadrature::size
unsigned int size() const
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const

Here are the time step, mass, and reaction parts:

((1.0 +
time_step / 2.0 * parameters.reaction_coefficient) *
fe_values.shape_value(i, q) *
fe_values.shape_value(j, q) +
time_step / 2.0 *

and the convection part:

(fe_values.shape_value(i, q) *
convection_contribution

and, finally, the diffusion part:

+ parameters.diffusion_coefficient *
(fe_values.shape_grad(i, q) *
fe_values.shape_grad(j, q))));
}
}
}
update_system_matrix(local_indices, cell_matrix);
}
}
}
template <int dim, typename MatrixType>
void
create_system_matrix(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
MatrixType & system_matrix)
{
internal_create_system_matrix<dim>(
dof_handler,
quad,
convection_function,
parameters,
time_step,
[&constraints, &system_matrix](
const std::vector<types::global_dof_index> &local_indices,
const FullMatrix<double> & cell_matrix) {
constraints.distribute_local_to_global(cell_matrix,
local_indices,
system_matrix);
});
}
template <int dim, typename MatrixType>
void
create_system_matrix(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
MatrixType & system_matrix)
{
internal_create_system_matrix<dim>(
dof_handler,
quad,
convection_function,
parameters,
time_step,
[&system_matrix](
const std::vector<types::global_dof_index> &local_indices,
const FullMatrix<double> & cell_matrix) {
system_matrix.add(local_indices, cell_matrix);
});
}
} // namespace CDR
#endif
AffineConstraints::distribute_local_to_global
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
Definition: affine_constraints.h:2334

Annotated version of common/include/deal.II-cdr/system_rhs.h

#ifndef dealii__cdr_system_rhs_h
#define dealii__cdr_system_rhs_h
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II-cdr/parameters.h>
#include <functional>

Similarly to create_system_matrix, I wrote a separate function to compute the right hand side.

namespace CDR
{
using namespace dealii;
template <int dim, typename VectorType>
void
create_system_rhs(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const std::function<double(double, const Point<dim>)> &forcing_function,
const CDR::Parameters & parameters,
const VectorType & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
VectorType & system_rhs);
} // namespace CDR
#endif

Annotated version of common/include/deal.II-cdr/system_rhs.templates.h

#ifndef dealii__cdr_system_rhs_templates_h
#define dealii__cdr_system_rhs_templates_h
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/vector.h>
#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_rhs.h>
#include <functional>
#include <vector>
namespace CDR
{
using namespace dealii;
template <int dim, typename VectorType>
void
create_system_rhs(
const DoFHandler<dim> & dof_handler,
const QGauss<dim> & quad,
const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
const std::function<double(double, const Point<dim>)> &forcing_function,
const CDR::Parameters & parameters,
const VectorType & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
VectorType & system_rhs)
{
auto & fe = dof_handler.get_fe();
const auto dofs_per_cell = fe.dofs_per_cell;
const double time_step =
(parameters.stop_time - parameters.start_time) / parameters.n_time_steps;
FEValues<dim> fe_values(fe,
quad,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
Vector<double> cell_rhs(dofs_per_cell);
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> current_fe_coefficients(dofs_per_cell);
std::vector<types::global_dof_index> local_indices(dofs_per_cell);
const double previous_time{current_time - time_step};
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
cell_rhs = 0.0;
cell->get_dof_indices(local_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
current_fe_coefficients[i] =
previous_solution[local_indices[i]];
}
for (unsigned int q = 0; q < quad.size(); ++q)
{
const auto current_convection =
convection_function(fe_values.quadrature_point(q));
const double current_forcing =
forcing_function(current_time, fe_values.quadrature_point(q));
const double previous_forcing =
forcing_function(previous_time,
fe_values.quadrature_point(q));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const auto convection_contribution =
current_convection * fe_values.shape_grad(j, q);
cell_rhs(i) +=
fe_values.JxW(q) *
Vector
Definition: vector.h:109

Here are the mass and reaction part:

(((1.0 - time_step / 2.0 *
parameters.reaction_coefficient) *
fe_values.shape_value(i, q) *
fe_values.shape_value(j, q) -
time_step / 2.0 *

the convection part:

(fe_values.shape_value(i, q) *
convection_contribution

the diffusion part:

+ parameters.diffusion_coefficient *
(fe_values.shape_grad(i, q) *
fe_values.shape_grad(j, q)))) *
current_fe_coefficients[j]

and, finally, the forcing function part:

+ time_step / 2.0 *
(current_forcing + previous_forcing) *
fe_values.shape_value(i, q));
}
}
}
constraints.distribute_local_to_global(cell_rhs,
local_indices,
system_rhs);
}
}
}
} // namespace CDR
#endif

Annotated version of common/include/deal.II-cdr/write_pvtu_output.h

#ifndef dealii__cdr_write_pvtu_output_h
#define dealii__cdr_write_pvtu_output_h
#include <deal.II/dofs/dof_handler.h>

This is a small class which handles PVTU output.

namespace CDR
{
using namespace dealii;
class WritePVTUOutput
{
public:
WritePVTUOutput(const unsigned int patch_level);
template <int dim, typename VectorType>
void
write_output(const DoFHandler<dim> &dof_handler,
const VectorType & solution,
const unsigned int time_step_n,
const double current_time);
private:
const unsigned int patch_level;
const unsigned int this_mpi_process;
};
} // namespace CDR
#endif

Annotated version of common/include/deal.II-cdr/write_pvtu_output.templates.h

#ifndef dealii__cdr_write_pvtu_output_templates_h
#define dealii__cdr_write_pvtu_output_templates_h
#include <deal.II/base/data_out_base.h>
#include <deal.II/base/utilities.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II-cdr/write_pvtu_output.h>
#include <fstream>
#include <string>
#include <vector>
data_out_base.h
data_out.h

Here is the implementation of the important function. This is similar to what is presented in step-40.

namespace CDR
{
using namespace dealii;
template <int dim, typename VectorType>
void
WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
const VectorType & solution,
const unsigned int time_step_n,
const double current_time)
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "u");
const auto & triangulation = dof_handler.get_triangulation();
Vector<float> subdomain(triangulation.n_active_cells());
for (auto &domain : subdomain)
{
domain = triangulation.locally_owned_subdomain();
}
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches(patch_level);
DataOutBase::VtkFlags flags;
flags.time = current_time;
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
DataOut_DoFData::add_data_vector
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={})
Definition: data_out_dof_data.h:1082
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1064
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DataOutBase::VtkFlags::time
double time
Definition: data_out_base.h:1102
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69

While the default flag is for the best compression level, using best_speed makes this function much faster.

flags.compression_level =
DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
data_out.set_flags(flags);
unsigned int subdomain_n;
if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1)
{
subdomain_n = 0;
}
else
{
subdomain_n = triangulation.locally_owned_subdomain();
}
std::ofstream output("solution-" + Utilities::int_to_string(time_step_n) +
"." + Utilities::int_to_string(subdomain_n, 4) +
".vtu");
data_out.write_vtu(output);
if (this_mpi_process == 0)
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
++i)
filenames.push_back("solution-" +
Utilities::int_to_string(time_step_n) + "." +
Utilities::int_to_string(i, 4) + ".vtu");
std::ofstream master_output(
"solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
} // namespace CDR
#endif
DataOutInterface::write_pvtu_record
void write_pvtu_record(std::ostream &out, const std::vector< std::string > &piece_names) const
Definition: data_out_base.cc:7820
DataOutBase::VtkFlags::compression_level
ZlibCompressionLevel compression_level
Definition: data_out_base.h:1154
DataOutInterface::write_vtu
void write_vtu(std::ostream &out) const
Definition: data_out_base.cc:7673
DataOutInterface::set_flags
void set_flags(const FlagType &flags)
Definition: data_out_base.cc:8680
DataOutBase::VtkFlags::best_speed
@ best_speed
Definition: data_out_base.h:1137
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:140
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473

Annotated version of common/source/parameters.cc

#include <deal.II-cdr/parameters.h>
#include <fstream>
#include <string>
namespace CDR
{
void
Parameters::configure_parameter_handler(ParameterHandler &parameter_handler)
{
parameter_handler.enter_subsection("Geometry");
{
parameter_handler.declare_entry("inner_radius",
"1.0",
Patterns::Double(0.0),
"Inner radius.");
parameter_handler.declare_entry("outer_radius",
"2.0",
Patterns::Double(0.0),
"Outer radius.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Physical Parameters");
{
parameter_handler.declare_entry("diffusion_coefficient",
"1.0",
Patterns::Double(0.0),
"Diffusion coefficient.");
parameter_handler.declare_entry("reaction_coefficient",
"1.0",
Patterns::Double(0.0),
"Reaction coefficient.");
parameter_handler.declare_entry("time_dependent_forcing",
"true",
Patterns::Bool(),
"Whether or not "
"the forcing function depends on time.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Finite Element");
{
parameter_handler.declare_entry("initial_refinement_level",
"1",
Patterns::Integer(1),
"Initial number of levels in the mesh.");
parameter_handler.declare_entry("max_refinement_level",
"1",
Patterns::Integer(1),
"Maximum number of levels in the mesh.");
parameter_handler.declare_entry("fe_order",
"1",
Patterns::Integer(1),
"Finite element order.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Time Step");
{
parameter_handler.declare_entry("start_time",
"0.0",
Patterns::Double(0.0),
"Start time.");
parameter_handler.declare_entry("stop_time",
"1.0",
Patterns::Double(1.0),
"Stop time.");
parameter_handler.declare_entry("n_time_steps",
"1",
Patterns::Integer(1),
"Number of time steps.");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Output");
{
parameter_handler.declare_entry("save_interval",
"10",
Patterns::Integer(1),
"Save interval.");
parameter_handler.declare_entry("patch_level",
"2",
Patterns::Integer(0),
"Patch level.");
}
parameter_handler.leave_subsection();
}
void
Parameters::read_parameter_file(const std::string &file_name)
{
ParameterHandler parameter_handler;
{
std::ifstream file(file_name);
configure_parameter_handler(parameter_handler);
parameter_handler.parse_input(file);
}
parameter_handler.enter_subsection("Geometry");
{
inner_radius = parameter_handler.get_double("inner_radius");
outer_radius = parameter_handler.get_double("outer_radius");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Physical Parameters");
{
diffusion_coefficient =
parameter_handler.get_double("diffusion_coefficient");
reaction_coefficient =
parameter_handler.get_double("reaction_coefficient");
time_dependent_forcing =
parameter_handler.get_bool("time_dependent_forcing");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Finite Element");
{
initial_refinement_level =
parameter_handler.get_integer("initial_refinement_level");
max_refinement_level =
parameter_handler.get_integer("max_refinement_level");
fe_order = parameter_handler.get_integer("fe_order");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Time Step");
{
start_time = parameter_handler.get_double("start_time");
stop_time = parameter_handler.get_double("stop_time");
n_time_steps = parameter_handler.get_integer("n_time_steps");
}
parameter_handler.leave_subsection();
parameter_handler.enter_subsection("Output");
{
save_interval = parameter_handler.get_integer("save_interval");
patch_level = parameter_handler.get_integer("patch_level");
}
parameter_handler.leave_subsection();
}
} // namespace CDR
ParameterHandler::parse_input
virtual void parse_input(std::istream &input, const std::string &filename="input file", const std::string &last_line="", const bool skip_undefined=false)
Definition: parameter_handler.cc:432
ParameterHandler::get_integer
long int get_integer(const std::string &entry_string) const
Definition: parameter_handler.cc:1058
ParameterHandler::get_bool
bool get_bool(const std::string &entry_name) const
Definition: parameter_handler.cc:1146
ParameterHandler::declare_entry
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
Definition: parameter_handler.cc:829
ParameterHandler::get_double
double get_double(const std::string &entry_name) const
Definition: parameter_handler.cc:1101
ParameterHandler::enter_subsection
void enter_subsection(const std::string &subsection)
Definition: parameter_handler.cc:972

Annotated version of common/source/system_matrix.cc

#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_matrix.h>
#include <deal.II-cdr/system_matrix.templates.h>
sparse_matrix.h
trilinos_sparse_matrix.h

This file exists just to build template specializations of create_system_matrix. Even though the solver is run in parallel with Trilinos objects, other serial solvers can use the same function without recompilation by compiling everything here just one time.

namespace CDR
{
using namespace dealii;
template void
create_system_matrix<2, SparseMatrix<double>>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
SparseMatrix<double> & system_matrix);
template void
create_system_matrix<3, SparseMatrix<double>>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
SparseMatrix<double> & system_matrix);
template void
create_system_matrix<2, SparseMatrix<double>>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
SparseMatrix<double> & system_matrix);
template void
create_system_matrix<3, SparseMatrix<double>>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
SparseMatrix<double> & system_matrix);
template void
create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
TrilinosWrappers::SparseMatrix & system_matrix);
template void
create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
TrilinosWrappers::SparseMatrix & system_matrix);
template void
create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
TrilinosWrappers::SparseMatrix & system_matrix);
template void
create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
const CDR::Parameters & parameters,
const double time_step,
const AffineConstraints<double> & constraints,
TrilinosWrappers::SparseMatrix & system_matrix);
} // namespace CDR

Annotated version of common/source/system_rhs.cc

#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II-cdr/system_rhs.templates.h>
trilinos_vector.h

Like system_matrix.cc, this file just compiles template specializations.

namespace CDR
{
using namespace dealii;
template void
create_system_rhs<2, Vector<double>>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
const std::function<double(double, const Point<2>)> &forcing_function,
const CDR::Parameters & parameters,
const Vector<double> & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
Vector<double> & system_rhs);
template void
create_system_rhs<3, Vector<double>>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
const std::function<double(double, const Point<3>)> &forcing_function,
const CDR::Parameters & parameters,
const Vector<double> & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
Vector<double> & system_rhs);
template void
create_system_rhs<2, TrilinosWrappers::MPI::Vector>(
const DoFHandler<2> & dof_handler,
const QGauss<2> & quad,
const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
const std::function<double(double, const Point<2>)> &forcing_function,
const CDR::Parameters & parameters,
const TrilinosWrappers::MPI::Vector & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
TrilinosWrappers::MPI::Vector & system_rhs);
template void
create_system_rhs<3, TrilinosWrappers::MPI::Vector>(
const DoFHandler<3> & dof_handler,
const QGauss<3> & quad,
const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
const std::function<double(double, const Point<3>)> &forcing_function,
const CDR::Parameters & parameters,
const TrilinosWrappers::MPI::Vector & previous_solution,
const AffineConstraints<double> & constraints,
const double current_time,
TrilinosWrappers::MPI::Vector & system_rhs);
} // namespace CDR

Annotated version of common/source/write_pvtu_output.cc

#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II-cdr/write_pvtu_output.templates.h>

Again, this file just compiles the constructor and also the templated functions.

namespace CDR
{
using namespace dealii;
WritePVTUOutput::WritePVTUOutput(const unsigned int patch_level)
: patch_level{patch_level}
, this_mpi_process{Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
{}
template void
WritePVTUOutput::write_output(const DoFHandler<2> & dof_handler,
const Vector<double> &solution,
const unsigned int time_step_n,
const double current_time);
template void
WritePVTUOutput::write_output(const DoFHandler<3> & dof_handler,
const Vector<double> &solution,
const unsigned int time_step_n,
const double current_time);
template void
WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
const TrilinosWrappers::MPI::Vector &solution,
const unsigned int time_step_n,
const double current_time);
template void
WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
const TrilinosWrappers::MPI::Vector &solution,
const unsigned int time_step_n,
const double current_time);
} // namespace CDR

Annotated version of solver/cdr.cc

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/numerics/error_estimator.h>
conditional_ostream.h
dof_tools.h
dynamic_sparsity_pattern.h
error_estimator.h
manifold_lib.h
grid_generator.h

These headers are for distributed computations:

#include <deal.II/base/index_set.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II-cdr/parameters.h>
#include <deal.II-cdr/system_matrix.h>
#include <deal.II-cdr/system_rhs.h>
#include <deal.II-cdr/write_pvtu_output.h>
#include <chrono>
#include <functional>
#include <iostream>
using namespace dealii;
constexpr int manifold_id{0};
grid_refinement.h
solution_transfer.h
index_set.h
sparsity_tools.h
trilinos_precondition.h
trilinos_solver.h

This is the actual solver class which performs time iteration and calls the appropriate library functions to do it.

template <int dim>
class CDRProblem
{
public:
CDRProblem(const CDR::Parameters &parameters);
void
run();
private:
const CDR::Parameters parameters;
const double time_step;
double current_time;
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
FE_Q<dim> fe;
QGauss<dim> quad;
const SphericalManifold<dim> boundary_description;
parallel::distributed::Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
const std::function<Tensor<1, dim>(const Point<dim>)> convection_function;
const std::function<double(double, const Point<dim>)> forcing_function;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
bool first_run;
FE_Q
Definition: fe_q.h:549
MPI_Comm

As is usual in parallel programs, I keep two copies of parts of the complete solution: locally_relevant_solution contains both the locally calculated solution as well as the layer of cells at its boundary (the ghost cells) while completely_distributed_solution only contains the parts of the solution computed on the current MPI process.

TrilinosWrappers::MPI::Vector locally_relevant_solution;
TrilinosWrappers::MPI::Vector completely_distributed_solution;
TrilinosWrappers::MPI::Vector system_rhs;
TrilinosWrappers::SparseMatrix system_matrix;
TrilinosWrappers::PreconditionAMG preconditioner;
ConditionalOStream pcout;
void
setup_geometry();
void
setup_system();
void
setup_dofs();
void
refine_mesh();
void
time_iterate();
};
template <int dim>
CDRProblem<dim>::CDRProblem(const CDR::Parameters &parameters)
: parameters(parameters)
, time_step{(parameters.stop_time - parameters.start_time) /
parameters.n_time_steps}
, current_time{parameters.start_time}
, mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes{Utilities::MPI::n_mpi_processes(mpi_communicator)}
, this_mpi_process{Utilities::MPI::this_mpi_process(mpi_communicator)}
, fe(parameters.fe_order)
, quad(parameters.fe_order + 2)
, boundary_description(Point<dim>())
, triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening))
, dof_handler(triangulation)
, convection_function{[](const Point<dim> p) -> Tensor<1, dim> {
Tensor<1, dim> v;
v[0] = -p[1];
v[1] = p[0];
return v;
}}
, forcing_function{[](double t, const Point<dim> p) -> double {
return std::exp(-8 * t) *
std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) *
std::exp(-40 * Utilities::fixed_power<6>(p[1]));
}}
, first_run{true}
, pcout(std::cout, this_mpi_process == 0)
{
Assert(dim == 2, ExcNotImplemented());
}
template <int dim>
void
CDRProblem<dim>::setup_geometry()
{
const Point<dim> center;
GridGenerator::hyper_shell(triangulation,
center,
parameters.inner_radius,
parameters.outer_radius);
triangulation.set_manifold(manifold_id, boundary_description);
for (const auto &cell : triangulation.active_cell_iterators())
{
cell->set_all_manifold_ids(manifold_id);
}
triangulation.refine_global(parameters.initial_refinement_level);
}
template <int dim>
void
CDRProblem<dim>::setup_dofs()
{
dof_handler.distribute_dofs(fe);
pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
DoFTools::make_zero_boundary_constraints(dof_handler,
manifold_id,
constraints);
constraints.close();
completely_distributed_solution.reinit(locally_owned_dofs, mpi_communicator);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
}
template <int dim>
void
CDRProblem<dim>::setup_system()
{
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dynamic_sparsity_pattern,
constraints,
/*keep_constrained_dofs*/ true);
SparsityTools::distribute_sparsity_pattern(dynamic_sparsity_pattern,
dof_handler.locally_owned_dofs(),
mpi_communicator,
locally_relevant_dofs);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
system_matrix.reinit(locally_owned_dofs,
dynamic_sparsity_pattern,
mpi_communicator);
CDR::create_system_matrix<dim>(dof_handler,
quad,
convection_function,
parameters,
time_step,
constraints,
system_matrix);
system_matrix.compress(VectorOperation::add);
preconditioner.initialize(system_matrix);
}
template <int dim>
void
CDRProblem<dim>::time_iterate()
{
double current_time = parameters.start_time;
CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
++time_step_n)
{
current_time += time_step;
system_rhs = 0.0;
CDR::create_system_rhs<dim>(dof_handler,
quad,
convection_function,
forcing_function,
parameters,
locally_relevant_solution,
constraints,
current_time,
system_rhs);
system_rhs.compress(VectorOperation::add);
SolverControl solver_control(dof_handler.n_dofs(),
1e-6 * system_rhs.l2_norm(),
/*log_history = */ false,
/*log_result = */ false);
TrilinosWrappers::SolverGMRES solver(solver_control);
solver.solve(system_matrix,
completely_distributed_solution,
system_rhs,
preconditioner);
constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
if (time_step_n % parameters.save_interval == 0)
{
pvtu_output.write_output(dof_handler,
locally_relevant_solution,
time_step_n,
current_time);
}
refine_mesh();
}
}
template <int dim>
void
CDRProblem<dim>::refine_mesh()
{
using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(dof_handler,
QGauss<dim - 1>(fe.degree + 1),
FunctionMap(),
locally_relevant_solution,
estimated_error_per_cell);
DoFHandler::distribute_dofs
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
DoFHandler::locally_owned_dofs
const IndexSet & locally_owned_dofs() const
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
KellyErrorEstimator::estimate
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)
center
Point< 3 > center
Definition: data_out_base.cc:187
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
VectorOperation::add
@ add
Definition: vector_operation.h:53
DoFTools::make_zero_boundary_constraints
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask=ComponentMask())
Definition: dof_tools_constraints.cc:3451
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1788
DoFTools::make_sparsity_pattern
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)
Definition: dof_tools_sparsity.cc:62
DoFTools::extract_locally_relevant_dofs
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1144
GridGenerator::hyper_shell
void hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
Definition: sparsity_tools.cc:1030
std::exp
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5466

This solver uses a crude refinement strategy where cells with relatively high errors are refined and cells with relatively low errors are coarsened. The maximum refinement level is capped to prevent run-away refinement.

for (const auto &cell : triangulation.active_cell_iterators())
{
if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) >= 1e-3)
{
cell->set_refine_flag();
}
else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <=
1e-5)
{
cell->set_coarsen_flag();
}
}
if (triangulation.n_levels() > parameters.max_refinement_level)
{
for (const auto &cell : triangulation.cell_iterators_on_level(
parameters.max_refinement_level))
{
cell->clear_refine_flag();
}
}
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5572

Transferring the solution between different grids is ultimately just a few function calls but they must be made in exactly the right order.

parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
solution_transfer(dof_handler);
triangulation.prepare_coarsening_and_refinement();
solution_transfer.prepare_for_coarsening_and_refinement(
locally_relevant_solution);
triangulation.execute_coarsening_and_refinement();
setup_dofs();

The solution_transfer object stores a pointer to locally_relevant_solution, so when parallel::distributed::SolutionTransfer::interpolate is called it uses those values to populate temporary.

TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator);
solution_transfer.interpolate(temporary);

After temporary has the correct value, this call correctly populates completely_distributed_solution, which had its index set updated above with the call to setup_dofs.

completely_distributed_solution = temporary;

Constraints cannot be applied to vectors with ghost entries since the ghost entries are write only, so this first goes through the completely distributed vector.

constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
setup_system();
}
template <int dim>
void
CDRProblem<dim>::run()
{
setup_geometry();
setup_dofs();
setup_system();
time_iterate();
}
constexpr int dim{2};
int
main(int argc, char *argv[])
{

One of the new features in C++11 is the chrono component of the standard library. This gives us an easy way to time the output.

auto t0 = std::chrono::high_resolution_clock::now();
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
CDR::Parameters parameters;
parameters.read_parameter_file("parameters.prm");
CDRProblem<dim> cdr_problem(parameters);
cdr_problem.run();
auto t1 = std::chrono::high_resolution_clock::now();
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
std::cout << "time elapsed: "
<< std::chrono::duration_cast<std::chrono::milliseconds>(t1 -
t0)
.count()
<< " milliseconds." << std::endl;
}
return 0;
}
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:980
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:151
Generated by doxygen 1.9.5