This tutorial depends on step-16, step-37, step-40.
This program was contributed by Thomas C. Clevenger and Timo Heister.
This material is based upon work partly supported by the National Science Foundation Award DMS-2028346, OAC-2015848, EAR-1925575, by the Computational Infrastructure in Geodynamics initiative (CIG), through the NSF under Award EAR-0949446 and EAR-1550901 and The University of California – Davis.
- Note
- As a prerequisite of this program, you need to have both p4est and either the PETSc or Trilinos library installed. The installation of deal.II together with these additional libraries is described in the README file.
Introduction
This example shows the usage of the multilevel functions in deal.II on parallel, distributed meshes and gives a comparison between geometric and algebraic multigrid methods. The algebraic multigrid (AMG) preconditioner is the same used in step-40. Two geometric multigrid (GMG) preconditioners are considered: a matrix-based version similar to that in step-16 (but for parallel computations) and a matrix-free version discussed in step-37. The goal is to find out which approach leads to the best solver for large parallel computations.
Algebraic multigrid methods are obviously the easiest to implement with deal.II since classes such as TrilinosWrappers::PreconditionAMG and PETScWrappers::PreconditionBoomerAMG are, in essence, black box preconditioners that require only a couple of lines to set up even for parallel computations. On the other hand, geometric multigrid methods require changes throughout a code base – not very many, but one has to know what one is doing.
What the results of this program will show is that algebraic and geometric multigrid methods are roughly comparable in performance when using matrix-based formulations, and that matrix-free geometric multigrid methods are vastly better for the problem under consideration here. A secondary conclusion will be that matrix-based geometric multigrid methods really don't scale well strongly when the number of unknowns per processor becomes smaller than 20,000 or so.
The testcase
We consider the variable-coefficient Laplacian weak formulation
\begin{align*} (\epsilon \nabla u, \nabla v) = (f,v) \quad \forall v \in V_h \end{align*}
on the domain \(\Omega = [-1,1]^\text{dim} \setminus [0,1]^\text{dim}\) (an L-shaped domain for 2D and a Fichera corner for 3D) with \(\epsilon = 1\) if \(\min(x,y,z)>-\frac{1}{2}\) and \(\epsilon = 100\) otherwise. In other words, \(\epsilon\) is small along the edges or faces of the domain that run into the reentrant corner, as will be visible in the figure below.
The boundary conditions are \(u=0\) on the whole boundary and the right-hand side is \(f=1\). We use continuous \(Q_2\) elements for the discrete finite element space \(V_h\), and use a residual-based, cell-wise a posteriori error estimator \(e(K) = e_{\text{cell}}(K) + e_{\text{face}}(K)\) from [42] with
\begin{align*} e_{\text{cell}}(K) &= h^2 \| f + \epsilon \triangle u \|_K^2, \\ e_{\text{face}}(K) &= \sum_F h_F \| \jump{ \epsilon \nabla u \cdot n } \|_F^2, \end{align*}
to adaptively refine the mesh. (This is a generalization of the Kelly error estimator used in the KellyErrorEstimator class that drives mesh refinement in most of the other tutorial programs.) The following figure visualizes the solution and refinement for 2D: 
In 3D, the solution looks similar (see below). On the left you can see the solution and on the right we show a slice for \(x\) close to the center of the domain showing the adaptively refined mesh.
Both in 2D and 3D you can see the adaptive refinement picking up the corner singularity and the inner singularity where the viscosity jumps, while the interface along the line that separates the two viscosities is (correctly) not refined as it is resolved adequately. This is because the kink in the solution that results from the jump in the coefficient is aligned with cell interfaces.
Workload imbalance for geometric multigrid methods
As mentioned above, the purpose of this program is to demonstrate the use of algebraic and geometric multigrid methods for this problem, and to do so for parallel computations. An important component of making algorithms scale to large parallel machines is ensuring that every processor has the same amount of work to do. (More precisely, what matters is that there are no small fraction of processors that have substantially more work than the rest since, if that were so, a large fraction of processors will sit idle waiting for the small fraction to finish. Conversely, a small fraction of processors having substantially less work is not a problem because the majority of processors continues to be productive and only the small fraction sits idle once finished with their work.)
For the active mesh, we use the parallel::distributed::Triangulation class as done in step-40 which uses functionality in the external library p4est for the distribution of the active cells among processors. For the non-active cells in the multilevel hierarchy, deal.II implements what we will refer to as the "first-child rule" where, for each cell in the hierarchy, we recursively assign the parent of a cell to the owner of the first child cell. The following figures give an example of such a distribution. Here the left image represents the active cells for a sample 2D mesh partitioned using a space-filling curve (which is what p4est uses to partition cells); the center image gives the tree representation of the active mesh; and the right image gives the multilevel hierarchy of cells. The colors and numbers represent the different processors. The circular nodes in the tree are the non-active cells which are distributed using the "first-child rule".
Included among the output to screen in this example is a value "Partition efficiency" given by one over MGTools::workload_imbalance(). This value, which will be denoted by \(\mathbb{E}\), quantifies the overhead produced by not having a perfect work balance on each level of the multigrid hierarchy. This imbalance is evident from the example above: while level \(\ell=2\) is about as well balanced as is possible with four cells among three processors, the coarse level \(\ell=0\) has work for only one processor, and level \(\ell=1\) has work for only two processors of which one has three times as much work as the other.
For defining \(\mathbb{E}\), it is important to note that, as we are using local smoothing to define the multigrid hierarchy (see the multigrid paper for a description of local smoothing), the refinement level of a cell corresponds to that cell's multigrid level. Now, let \(N_{\ell}\) be the number of cells on level \(\ell\) (both active and non-active cells) and \(N_{\ell,p}\) be the subset owned by process \(p\). We will also denote by \(P\) the total number of processors. Assuming that the workload for any one processor is proportional to the number of cells owned by that processor, the optimal workload per processor is given by
\begin{align*} W_{\text{opt}} = \frac1{P}\sum_{\ell} N_{\ell} = \sum_{\ell}\left(\frac1{P}\sum_{p}N_{\ell,p}\right). \end{align*}
Next, assuming a synchronization of work on each level (i.e., on each level of a V-cycle, work must be completed by all processors before moving on to the next level), the limiting effort on each level is given by
\begin{align*} W_\ell = \max_{p} N_{\ell,p}, \end{align*}
and the total parallel complexity
\begin{align*} W = \sum_{\ell} W_\ell. \end{align*}
Then we define \(\mathbb{E}\) as a ratio of the optimal partition to the parallel complexity of the current partition
\begin{align*} \mathbb{E} = \frac{W_{\text{opt}}}{W}. \end{align*}
For the example distribution above, we have
\begin{align*} W_{\text{opt}}&=\frac{1}{P}\sum_{\ell} N_{\ell} = \frac{1}{3} \left(1+4+4\right)= 3 \qquad \\ W &= \sum_\ell W_\ell = 1 + 2 + 3 = 6 \\ \mathbb{E} &= \frac{W_{\text{opt}}}{W} = \frac12. \end{align*}
The value MGTools::workload_imbalance() \(= 1/\mathbb{E}\) then represents the factor increase in timings we expect for GMG methods (vmults, assembly, etc.) due to the imbalance of the mesh partition compared to a perfectly load-balanced workload. We will report on these in the results section below for a sequence of meshes, and compare with the observed slow-downs as we go to larger and larger processor numbers (where, typically, the load imbalance becomes larger as well).
These sorts of considerations are considered in much greater detail in [16], which contains a full discussion of the partition efficiency model and the effect the imbalance has on the GMG V-cycle timing. In summary, the value of \(\mathbb{E}\) is highly dependent on the degree of local mesh refinement used and has an optimal value \(\mathbb{E} \approx 1\) for globally refined meshes. Typically for adaptively refined meshes, the number of processors used to distribute a single mesh has a negative impact on \(\mathbb{E}\) but only up to a leveling off point, where the imbalance remains relatively constant for an increasing number of processors, and further refinement has very little impact on \(\mathbb{E}\). Finally, \(1/\mathbb{E}\) was shown to give an accurate representation of the slowdown in parallel scaling expected for the timing of a V-cycle.
It should be noted that there is potential for some asynchronous work between multigrid levels, specifically with purely nearest neighbor MPI communication, and an adaptive mesh could be constructed such that the efficiency model would far overestimate the V-cycle slowdown due to the asynchronous work "covering up" the imbalance (which assumes synchronization over levels). However, for most realistic adaptive meshes the expectation is that this asynchronous work will only cover up a very small portion of the imbalance and the efficiency model will describe the slowdown very well.
Workload imbalance for algebraic multigrid methods
The considerations above show that one has to expect certain limits on the scalability of the geometric multigrid algorithm as it is implemented in deal.II because even in cases where the finest levels of a mesh are perfectly load balanced, the coarser levels may not be. At the same time, the coarser levels are weighted less (the contributions of \(W_\ell\) to \(W\) are small) because coarser levels have fewer cells and, consequently, do not contribute to the overall run time as much as finer levels. In other words, imbalances in the coarser levels may not lead to large effects in the big picture.
Algebraic multigrid methods are of course based on an entirely different approach to creating a hierarchy of levels. In particular, they create these purely based on analyzing the system matrix, and very sophisticated algorithms for ensuring that the problem is well load-balanced on every level are implemented in both the hypre and ML/MueLu packages that underly the TrilinosWrappers::PreconditionAMG and PETScWrappers::PreconditionBoomerAMG classes. In some sense, these algorithms are simpler than for geometric multigrid methods because they only deal with the matrix itself, rather than all of the connotations of meshes, neighbors, parents, and other geometric entities. At the same time, much work has also been put into making algebraic multigrid methods scale to very large problems, including questions such as reducing the number of processors that work on a given level of the hierarchy to a subset of all processors, if otherwise processors would spend less time on computations than on communication. (One might note that it is of course possible to implement these same kinds of ideas also in geometric multigrid algorithms where one purposefully idles some processors on coarser levels to reduce the amount of communication. deal.II just doesn't do this at this time.)
These are not considerations we typically have to worry about here, however: For most purposes, we use algebraic multigrid methods as black-box methods.
Running the program
As mentioned above, this program can use three different ways of solving the linear system: matrix-based geometric multigrid ("MB"), matrix-free geometric multigrid ("MF"), and algebraic multigrid ("AMG"). The directory in which this program resides has input files with suffix .prm for all three of these options, and for both 2d and 3d.
You can execute the program as in
and this will take the run-time parameters from the given input file (here, gmg_mb_2d.prm).
The program is intended to be run in parallel, and you can achieve this using a command such as
mpirun -np 4 ./step-50 gmg_mb_2d.prm
if you want to, for example, run on four processors. (That said, the program is also ready to run with, say, -np 28672 if that's how many processors you have available.)
The commented program
Include files
The include files are a combination of step-40, step-16, and step-37:
We use the same strategy as in step-40 to switch between PETSc and Trilinos:
Comment the following preprocessor definition in or out if you have PETSc and Trilinos installed and you prefer using PETSc in this example:
#define FORCE_USE_OF_TRILINOS
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}
The following files are used to assemble the error estimator like in step-12:
Coefficients and helper classes
MatrixFree operators must use the LinearAlgebra::distributed::Vector vector type. Here we define operations which copy to and from Trilinos vectors for compatibility with the matrix-based code. Note that this functionality does not currently exist for PETSc vector types, so Trilinos must be installed to use the MatrixFree solver in this tutorial.
namespace ChangeVectorTypes
{
template <typename number>
const ::LinearAlgebra::distributed::Vector<number> &in)
{
out.locally_owned_elements());
#ifdef USE_PETSC_LA
ExcMessage(
"CopyVectorTypes::copy() not implemented for "
"PETSc vector types."));
#else
#endif
}
template <typename number>
{
#ifdef USE_PETSC_LA
(void)in;
ExcMessage(
"CopyVectorTypes::copy() not implemented for "
"PETSc vector types."));
#else
#endif
}
}
Let's move on to the description of the problem we want to solve. We set the right-hand side function to 1.0. The value function returning a VectorizedArray is used by the matrix-free code path.
template <int dim>
class RightHandSide :
public Function<dim>
{
public:
const unsigned int = 0) const override
{
return 1.0;
}
template <typename number>
const unsigned int = 0) const
{
}
};
This next class represents the diffusion coefficient. We use a variable coefficient which is 100.0 at any point where at least one coordinate is less than -0.5, and 1.0 at all other points. As above, a separate value() returning a VectorizedArray is used for the matrix-free code. An average() function computes the arithmetic average for a set of points.
template <int dim>
class Coefficient :
public Function<dim>
{
public:
const unsigned int = 0) const override;
template <typename number>
const unsigned int = 0) const;
template <typename number>
When using a coefficient in the MatrixFree framework, we also need a function that creates a Table of coefficient values for a set of cells provided by the MatrixFree operator argument here.
template <typename number>
std::shared_ptr<Table<2, VectorizedArray<number>>> make_coefficient_table(
};
template <int dim>
{
for (
int d = 0;
d < dim; ++
d)
{
return 100.0;
}
return 1.0;
}
template <int dim>
template <typename number>
const unsigned int) const
{
for (unsigned int i = 0; i < VectorizedArray<number>::size(); ++i)
{
for (
int d = 0;
d < dim; ++
d)
{
return_value[i] = 100.0;
break;
}
}
return return_value;
}
template <int dim>
template <typename number>
number Coefficient<dim>::average_value(
{
number average(0);
for (unsigned int i = 0; i < points.size(); ++i)
average +=
value(points[i]);
average /= points.size();
return average;
}
template <int dim>
template <typename number>
std::shared_ptr<Table<2, VectorizedArray<number>>>
Coefficient<dim>::make_coefficient_table(
{
auto coefficient_table =
std::make_shared<Table<2, VectorizedArray<number>>>();
const unsigned int n_cells = mf_storage.n_macro_cells();
const unsigned int n_q_points = fe_eval.n_q_points;
coefficient_table->reinit(
n_cells, 1);
for (
unsigned int cell = 0; cell <
n_cells; ++cell)
{
fe_eval.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++q)
average_value +=
value(fe_eval.quadrature_point(q));
average_value /= n_q_points;
(*coefficient_table)(cell, 0) = average_value;
}
return coefficient_table;
}
Run time parameters
We will use ParameterHandler to pass in parameters at runtime. The structure Settings parses and stores these parameters to be queried throughout the program.
{
bool try_parse(const std::string &prm_filename);
{
gmg_mb,
gmg_mf,
amg
};
int dimension;
double smoother_dampen;
unsigned int smoother_steps;
unsigned int n_steps;
bool output;
};
bool Settings::try_parse(const std::string &prm_filename)
{
"10",
"Number of adaptive refinement steps.");
"1.0",
"Dampen factor for the smoother.");
"1",
"Number of smoother steps.");
"MF",
"Switch between matrix-free GMG, "
"matrix-based GMG, and AMG.");
"false",
"Output graphical results.");
if (prm_filename.size() == 0)
{
std::cout << "**** Error: No input file provided!\n"
<< "**** Error: Call this program as './step-50 input.prm\n"
<< "\n"
<< "**** You may want to use one of the input files in this\n"
<< "**** directory, or use the following default values\n"
<< "**** to create an input file:\n";
return false;
}
try
{
}
catch (std::exception &
e)
{
std::cerr <<
e.what() << std::endl;
return false;
}
if (prm.
get(
"solver") ==
"MF")
this->solver = gmg_mf;
else if (prm.
get(
"solver") ==
"MB")
this->solver = gmg_mb;
else if (prm.
get(
"solver") ==
"AMG")
this->solver = amg;
else
this->smoother_dampen = prm.
get_double(
"smoother dampen");
this->smoother_steps = prm.
get_integer(
"smoother steps");
return true;
}
LaplaceProblem class
This is the main class of the program. It looks very similar to step-16, step-37, and step-40. For the MatrixFree setup, we use the MatrixFreeOperators::LaplaceOperator class which defines local_apply(), compute_diagonal(), and set_coefficient() functions internally. Note that the polynomial degree is a template parameter of this class. This is necesary for the matrix-free code.
template <int dim, int degree>
class LaplaceProblem
{
public:
LaplaceProblem(
const Settings &settings);
private:
We will use the following types throughout the program. First the matrix-based types, after that the matrix-free classes. For the matrix-free implementation, we use float for the level operators.
dim,
degree,
degree + 1,
1,
dim,
degree,
degree + 1,
1,
void setup_system();
void setup_multigrid();
void assemble_system();
void assemble_multigrid();
void assemble_rhs();
void solve();
void estimate();
void refine_grid();
void output_results(const unsigned int cycle);
MatrixType system_matrix;
MatrixFreeActiveMatrix mf_system_matrix;
};
The only interesting part about the constructor is that we construct the multigrid hierarchy unless we use AMG. For that, we need to parse the run time parameters before this constructor completes.
template <int dim, int degree>
LaplaceProblem<dim, degree>::LaplaceProblem(
const Settings &settings)
: settings(settings)
, mpi_communicator(MPI_COMM_WORLD)
, mapping()
, fe(degree)
{
}
LaplaceProblem::setup_system()
Unlike step-16 and step-37, we split the set up into two parts, setup_system() and setup_multigrid(). Here is the typical setup_system() function for the active mesh found in most tutorials. For matrix-free, the active mesh set up is similar to step-37; for matrix-based (GMG and AMG solvers), the setup is similar to step-40.
template <int dim, int degree>
void LaplaceProblem<dim, degree>::setup_system()
{
solution.
reinit(locally_owned_dofs, mpi_communicator);
right_hand_side.reinit(locally_owned_dofs, mpi_communicator);
constraints.reinit(locally_relevant_dofs);
constraints.close();
switch (settings.solver)
{
case Settings::gmg_mf:
{
std::shared_ptr<MatrixFree<dim, double>> mf_storage =
std::make_shared<MatrixFree<dim, double>>();
mf_storage->reinit(dof_handler,
constraints,
additional_data);
mf_system_matrix.initialize(mf_storage);
const Coefficient<dim> coefficient;
mf_system_matrix.set_coefficient(
coefficient.make_coefficient_table(*mf_storage));
break;
}
case Settings::gmg_mb:
case Settings::amg:
{
#ifdef USE_PETSC_LA
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
#else
locally_owned_dofs,
locally_relevant_dofs,
MPI_COMM_WORLD);
system_matrix.reinit(dsp);
#endif
break;
}
default:
}
}
LaplaceProblem::setup_multigrid()
This function does the multilevel setup for both matrix-free and matrix-based GMG. The matrix-free setup is similar to that of step-37, and the matrix-based is similar to step-16, except we must use appropriate distributed sparsity patterns.
The function is not called for the AMG approach, but to err on the safe side, the main switch statement of this function nevertheless makes sure that the function only operates on known multigrid settings by throwing an assertion if the function were called for anything other than the two geometric multigrid methods.
template <int dim, int degree>
void LaplaceProblem<dim, degree>::setup_multigrid()
{
mg_constrained_dofs.clear();
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, boundary_ids);
switch (settings.solver)
{
case Settings::gmg_mf:
{
mf_mg_matrix.resize(0, n_levels - 1);
{
relevant_dofs);
level_constraints.
reinit(relevant_dofs);
mg_constrained_dofs.get_boundary_indices(
level));
level_constraints.
close();
std::shared_ptr<MatrixFree<dim, float>> mf_storage_level(
mf_storage_level->reinit(dof_handler,
level_constraints,
additional_data);
mf_mg_matrix[
level].initialize(mf_storage_level,
mg_constrained_dofs,
const Coefficient<dim> coefficient;
mf_mg_matrix[
level].set_coefficient(
coefficient.make_coefficient_table(*mf_storage_level));
mf_mg_matrix[
level].compute_diagonal();
}
break;
}
case Settings::gmg_mb:
{
mg_matrix.resize(0, n_levels - 1);
mg_matrix.clear_elements();
mg_interface_in.resize(0, n_levels - 1);
mg_interface_in.clear_elements();
{
dof_set);
{
#ifdef USE_PETSC_LA
dsp,
mpi_communicator,
dof_set);
dsp,
mpi_communicator);
#else
dof_set,
mpi_communicator);
mg_matrix[
level].reinit(dsp);
#endif
}
{
#ifdef USE_PETSC_LA
mg_constrained_dofs,
dsp,
dsp,
mpi_communicator,
dof_set);
mg_interface_in[
level].reinit(
dsp,
mpi_communicator);
#else
dof_set,
mpi_communicator);
mg_constrained_dofs,
dsp,
mg_interface_in[
level].reinit(dsp);
#endif
}
}
break;
}
default:
}
}
LaplaceProblem::assemble_system()
The assembly is split into three parts: assemble_system(), assemble_multigrid(), and assemble_rhs(). The assemble_system() function here assembles and stores the (global) system matrix and the right-hand side for the matrix-based methods. It is similar to the assembly in step-40.
Note that the matrix-free method does not execute this function as it does not need to assemble a matrix, and it will instead assemble the right-hand side in assemble_rhs().
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_system()
{
quadrature_formula,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.
size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const Coefficient<dim> coefficient;
RightHandSide<dim> rhs;
std::vector<double> rhs_values(n_q_points);
if (cell->is_locally_owned())
{
cell_rhs = 0;
const double coefficient_value =
coefficient.average_value(fe_values.get_quadrature_points());
rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
coefficient_value *
fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) *
fe_values.JxW(q_point);
cell_rhs(i) +=
fe_values.shape_value(i, q_point) *
rhs_values[q_point] *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
cell_rhs,
local_dof_indices,
system_matrix,
right_hand_side);
}
}
LaplaceProblem::assemble_multigrid()
The following function assembles and stores the multilevel matrices for the matrix-based GMG method. This function is similar to the one found in step-16, only here it works for distributed meshes. This difference amounts to adding a condition that we only assemble on locally owned level cells and a call to compress() for each matrix that is built.
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_multigrid()
{
quadrature_formula,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.
size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const Coefficient<dim> coefficient;
std::vector<AffineConstraints<double>> boundary_constraints(
{
dof_set);
boundary_constraints[
level].reinit(dof_set);
boundary_constraints[
level].add_lines(
mg_constrained_dofs.get_refinement_edge_indices(
level));
boundary_constraints[
level].add_lines(
mg_constrained_dofs.get_boundary_indices(
level));
boundary_constraints[
level].close();
}
if (cell->level_subdomain_id() ==
triangulation.locally_owned_subdomain())
{
const double coefficient_value =
coefficient.average_value(fe_values.get_quadrature_points());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
coefficient_value * fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point);
cell->get_mg_dof_indices(local_dof_indices);
boundary_constraints[cell->level()].distribute_local_to_global(
cell_matrix, local_dof_indices, mg_matrix[cell->level()]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
if (mg_constrained_dofs.is_interface_matrix_entry(
cell->level(), local_dof_indices[i], local_dof_indices[j]))
mg_interface_in[cell->level()].add(local_dof_indices[i],
local_dof_indices[j],
}
for (
unsigned int i = 0; i <
triangulation.n_global_levels(); ++i)
{
}
}
LaplaceProblem::assemble_rhs()
The final function in this triptych assembles the right-hand side vector for the matrix-free method – because in the matrix-free framework, we don't have to assemble the matrix and can get away with only assembling the right hand side. We could do this by extracting the code from the assemble_system() function above that deals with the right hand side, but we decide instead to go all in on the matrix-free approach and do the assembly using that way as well.
The result is a function that is similar to the one found in the "Use FEEvaluation::read_dof_values_plain()
to avoid resolving constraints" subsection in the "Possibilities
for extensions" section of step-37.
The reason for this function is that the MatrixFree operators do not take into account non-homogeneous Dirichlet constraints, instead treating all Dirichlet constraints as homogeneous. To account for this, the right-hand side here is assembled as the residual \(r_0 = f-Au_0\), where \(u_0\) is a zero vector except in the Dirichlet values. Then when solving, we have that the solution is \(u = u_0 + A^{-1}r_0\). This can be seen as a Newton iteration on a linear system with initial guess \(u_0\). The CG solve in the solve() function below computes \(A^{-1}r_0\) and the call to constraints.distribute() (which directly follows) adds the \(u_0\).
Obviously, since we are considering a problem with zero Dirichlet boundary, we could have taken a similar approach to step-37 assemble_rhs(), but this additional work allows us to change the problem declaration if we so choose.
This function has two parts in the integration loop: applying the negative of matrix \(A\) to \(u_0\) by submitting the negative of the gradient, and adding the right-hand side contribution by submitting the value \(f\). We must be sure to use read_dof_values_plain() for evaluating \(u_0\) as read_dof_vaues() would set all Dirichlet values to zero.
Finally, the system_rhs vector is of type LA::MPI::Vector, but the MatrixFree class only work for LinearAlgebra::distributed::Vector. Therefore we must compute the right-hand side using MatrixFree funtionality and then use the functions in the ChangeVectorType namespace to copy it to the correct type.
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_rhs()
{
MatrixFreeActiveVector solution_copy;
MatrixFreeActiveVector right_hand_side_copy;
mf_system_matrix.initialize_dof_vector(solution_copy);
mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
solution_copy = 0.;
constraints.distribute(solution_copy);
solution_copy.update_ghost_values();
right_hand_side_copy = 0;
*(mf_system_matrix.get_coefficient());
RightHandSide<dim> right_hand_side_function;
*mf_system_matrix.get_matrix_free());
for (unsigned int cell = 0;
cell < mf_system_matrix.get_matrix_free()->n_macro_cells();
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(solution_copy);
phi.evaluate(false, true, false);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
phi.submit_gradient(-1.0 *
(coefficient(cell, 0) * phi.get_gradient(q)),
q);
phi.submit_value(
right_hand_side_function.value(phi.quadrature_point(q)), q);
}
phi.integrate_scatter(true, true, right_hand_side_copy);
}
}
LaplaceProblem::solve()
Here we set up the multigrid preconditioner, test the timing of a single V-cycle, and solve the linear system. Unsurprisingly, this is one of the places where the three methods differ the most.
template <int dim, int degree>
void LaplaceProblem<dim, degree>::solve()
{
SolverControl solver_control(1000, 1.
e-10 * right_hand_side.l2_norm());
solver_control.enable_history_data();
solution = 0.;
The solver for the matrix-free GMG method is similar to step-37, apart from adding some interface matrices in complete analogy to step-16.
switch (settings.solver)
{
case Settings::gmg_mf:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
mg_transfer.
build(dof_handler);
MatrixFreeLevelMatrix,
coarse_grid_solver(coarse_solver, mf_mg_matrix[0],
identity);
Smoother,
MatrixFreeLevelVector>
smoother;
typename Smoother::AdditionalData(
settings.smoother_dampen));
smoother.set_steps(settings.smoother_steps);
mg_interface_matrices;
mg_interface_matrices[
level].initialize(mf_mg_matrix[
level]);
mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
MatrixFreeLevelVector,
preconditioner(dof_handler,
mg, mg_transfer);
Copy the solution vector and right-hand side from LA::MPI::Vector to LinearAlgebra::distributed::Vector so that we can solve.
MatrixFreeActiveVector solution_copy;
MatrixFreeActiveVector right_hand_side_copy;
mf_system_matrix.initialize_dof_vector(solution_copy);
mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
computing_timer.leave_subsection("Solve: Preconditioner setup");
Timing for 1 V-cycle.
{
"Solve: 1 multigrid V-cycle");
preconditioner.
vmult(solution_copy, right_hand_side_copy);
}
solution_copy = 0.;
Solve the linear system, update the ghost values of the solution, copy back to LA::MPI::Vector and distribute constraints.
{
solver.solve(mf_system_matrix,
solution_copy,
right_hand_side_copy,
preconditioner);
}
solution_copy.update_ghost_values();
constraints.distribute(solution);
break;
}
Solver for the matrix-based GMG method, similar to step-16, only using a Jacobi smoother instead of a SOR smoother (which is not implemented in parallel).
case Settings::gmg_mb:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
mg_transfer.
build(dof_handler);
MatrixType,
coarse_grid_solver(coarse_solver, mg_matrix[0],
identity);
#ifdef USE_PETSC_LA
settings.smoother_dampen == 1.0,
"PETSc's PreconditionJacobi has no support for a damping parameter."));
#else
smoother.
initialize(mg_matrix, settings.smoother_dampen);
#endif
mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
mg.set_edge_matrices(mg_out, mg_in);
preconditioner(dof_handler,
mg, mg_transfer);
computing_timer.leave_subsection("Solve: Preconditioner setup");
Timing for 1 V-cycle.
{
"Solve: 1 multigrid V-cycle");
preconditioner.
vmult(solution, right_hand_side);
}
solution = 0.;
Solve the linear system and distribute constraints.
{
solver.solve(system_matrix,
solution,
right_hand_side,
preconditioner);
}
constraints.distribute(solution);
break;
}
Solver for the AMG method, similar to step-40.
case Settings::amg:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
PreconditionAMG::AdditionalData Amg_data;
#ifdef USE_PETSC_LA
Amg_data.symmetric_operator = true;
#else
Amg_data.smoother_type = "Jacobi";
Amg_data.higher_order_elements = true;
Amg_data.smoother_sweeps = settings.smoother_steps;
Amg_data.aggregation_threshold = 0.02;
#endif
Amg_data.output_details = false;
preconditioner.initialize(system_matrix, Amg_data);
computing_timer.leave_subsection("Solve: Preconditioner setup");
Timing for 1 V-cycle.
{
"Solve: 1 multigrid V-cycle");
preconditioner.vmult(solution, right_hand_side);
}
solution = 0.;
Solve the linear system and distribute constraints.
{
solver.solve(system_matrix,
solution,
right_hand_side,
preconditioner);
}
constraints.distribute(solution);
break;
}
default:
}
pcout << " Number of CG iterations: " << solver_control.last_step()
<< std::endl;
}
The error estimator
We use the FEInterfaceValues class to assemble an error estimator to decide which cells to refine. See the exact definition of the cell and face integrals in the introduction. To use the method, we define Scratch and Copy objects for the MeshWorker::mesh_loop() with much of the following code being in essence as was set up in step-12 already (or at least similar in spirit).
template <int dim>
struct ScratchData
{
const unsigned int quadrature_degree,
: fe_values(mapping, fe,
QGauss<dim>(quadrature_degree), update_flags)
, fe_interface_values(mapping,
fe,
QGauss<dim - 1>(quadrature_degree),
interface_update_flags)
{}
ScratchData(const ScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
scratch_data.fe_values.get_update_flags())
, fe_interface_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_interface_values.get_quadrature(),
scratch_data.fe_interface_values.get_update_flags())
{}
};
struct CopyData
{
CopyData()
{}
CopyData(const CopyData &) = default;
struct FaceData
{
unsigned int cell_indices[2];
double values[2];
};
unsigned int cell_index;
std::vector<FaceData> face_data;
};
template <int dim, int degree>
void LaplaceProblem<dim, degree>::estimate()
{
temp_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
temp_solution = solution;
const Coefficient<dim> coefficient;
Assembler for cell residual \(h \| f + \epsilon \triangle u \|_K\)
auto cell_worker = [&](
const Iterator & cell,
ScratchData<dim> &scratch_data,
CopyData & copy_data) {
RightHandSide<dim> rhs;
const double rhs_value = rhs.value(cell->center());
const double nu = coefficient.value(cell->center());
std::vector<Tensor<2, dim>> hessians(fe_values.n_quadrature_points);
fe_values.get_function_hessians(temp_solution, hessians);
copy_data.cell_index = cell->active_cell_index();
double residual_norm_square = 0.;
for (unsigned k = 0; k < fe_values.n_quadrature_points; ++k)
{
const double residual = (rhs_value + nu *
trace(hessians[k]));
residual_norm_square += residual * residual * fe_values.JxW(k);
}
copy_data.value = cell->diameter() *
std::sqrt(residual_norm_square);
};
Assembler for face term \(\sum_F h_F \| \jump{\epsilon \nabla u \cdot n} \|_F^2\)
auto face_worker = [&](
const Iterator & cell,
const unsigned int &f,
const unsigned int &sf,
const unsigned int &nf,
const unsigned int &nsf,
ScratchData<dim> & scratch_data,
CopyData & copy_data) {
scratch_data.fe_interface_values;
fe_interface_values.
reinit(cell, f, sf, ncell, nf, nsf);
copy_data.face_data.emplace_back();
CopyData::FaceData ©_data_face = copy_data.face_data.back();
copy_data_face.cell_indices[0] = cell->active_cell_index();
copy_data_face.cell_indices[1] = ncell->active_cell_index();
const double coeff1 = coefficient.value(cell->center());
const double coeff2 = coefficient.value(ncell->center());
std::vector<Tensor<1, dim>> grad_u[2];
for (unsigned int i = 0; i < 2; ++i)
{
temp_solution, grad_u[i]);
}
double jump_norm_square = 0.;
for (unsigned int qpoint = 0;
++qpoint)
{
const double jump =
coeff1 * grad_u[0][qpoint] * fe_interface_values.
normal(qpoint) -
coeff2 * grad_u[1][qpoint] * fe_interface_values.
normal(qpoint);
jump_norm_square += jump * jump * fe_interface_values.
JxW(qpoint);
}
const double h = cell->face(f)->measure();
copy_data_face.values[0] = 0.5 * h *
std::sqrt(jump_norm_square);
copy_data_face.values[1] = copy_data_face.values[0];
};
auto copier = [&](const CopyData ©_data) {
estimated_error_per_cell[copy_data.cell_index] += copy_data.value;
for (auto &cdf : copy_data.face_data)
for (unsigned int j = 0; j < 2; ++j)
estimated_error_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
const unsigned int n_gauss_points = degree + 1;
ScratchData<dim> scratch_data(mapping,
fe,
n_gauss_points,
CopyData copy_data;
We need to assemble each interior face once but we need to make sure that both processes assemble the face term between a locally owned and a ghost cell. This is achieved by setting the MeshWorker::assemble_ghost_faces_both flag. We need to do this, because we do not communicate the error estimator contributions here.
cell_worker,
copier,
scratch_data,
copy_data,
nullptr,
face_worker);
}
LaplaceProblem::refine_grid()
We use the cell-wise estimator stored in the vector estimate_vector and refine a fixed number of cells (chosen here to roughly double the number of DoFs in each step).
template <int dim, int degree>
void LaplaceProblem<dim, degree>::refine_grid()
{
const double refinement_fraction = 1. / (
std::pow(2.0, dim) - 1.);
triangulation, estimated_error_per_cell, refinement_fraction, 0.0);
}
LaplaceProblem::output_results()
The output_results() function is similar to the ones found in many of the tutorials (see step-40 for example).
template <int dim, int degree>
void LaplaceProblem<dim, degree>::output_results(const unsigned int cycle)
{
temp_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
temp_solution = solution;
for (unsigned int i = 0; i < subdomain.size(); ++i)
level(cell->active_cell_index()) = cell->level();
if (estimated_error_per_cell.size() > 0)
"estimated_error_per_cell");
"", "solution", cycle, mpi_communicator, 2 , 1 );
pcout << " Wrote " << master << std::endl;
}
As in most tutorials, this function calls the various functions defined above to setup, assemble, solve, and output the results.
template <int dim, int degree>
{
for (unsigned int cycle = 0; cycle < settings.n_steps; ++cycle)
{
pcout << "Cycle " << cycle << ':' << std::endl;
if (cycle > 0)
refine_grid();
pcout << " Number of active cells: "
We only output level cell data for the GMG methods (same with DoF data below). Note that the partition efficiency is irrelevant for AMG since the level hierarchy is not distributed or used during the computation.
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
pcout <<
" (" <<
triangulation.n_global_levels() <<
" global levels)" << std::endl
<< " Partition efficiency: "
pcout << std::endl;
setup_system();
Only set up the multilevel hierarchy for GMG.
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
setup_multigrid();
pcout <<
" Number of degrees of freedom: " << dof_handler.
n_dofs();
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
{
pcout << " (by level: ";
", ");
}
pcout << std::endl;
For the matrix-free method, we only assemble the right-hand side. For both matrix-based methods, we assemble both active matrix and right-hand side, and only assemble the multigrid matrices for matrix-based GMG.
if (settings.solver == Settings::gmg_mf)
assemble_rhs();
else
{
assemble_system();
if (settings.solver == Settings::gmg_mb)
assemble_multigrid();
}
solve();
estimate();
if (settings.output)
output_results(cycle);
computing_timer.print_summary();
computing_timer.reset();
}
}
The main() function
This is a similar main function to step-40, with the exception that we require the user to pass a .prm file as a sole command line argument (see step-29 and the documentation of the ParameterHandler class for a complete discussion of parameter files).
int main(int argc, char *argv[])
{
if (!settings.try_parse((argc > 1) ? (argv[1]) : ""))
return 0;
try
{
constexpr unsigned int fe_degree = 2;
switch (settings.dimension)
{
case 2:
{
LaplaceProblem<2, fe_degree> test(settings);
test.run();
break;
}
case 3:
{
LaplaceProblem<3, fe_degree> test(settings);
test.run();
break;
}
default:
}
}
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;
MPI_Abort(MPI_COMM_WORLD, 1);
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
MPI_Abort(MPI_COMM_WORLD, 2);
return 1;
}
return 0;
}
Results
When you run the program, the screen output should look like the following:
Cycle 0:
Number of active cells: 56 (2 global levels)
Workload imbalance: 1.14286
Number of degrees of freedom: 665 (by
level: 117, 665)
Number of CG iterations: 10
Wrote solution_00.pvtu
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 0.0457s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble right hand side | 1 | 0.000241s | 0.53% |
| Estimate | 1 | 0.0288s | 63% |
| Output results | 1 | 0.00219s | 4.8% |
| Setup | 1 | 0.00264s | 5.8% |
| Setup multigrid | 1 | 0.00261s | 5.7% |
| Solve | 1 | 0.00355s | 7.8% |
| Solve: 1 multigrid
V-cycle | 1 | 0.000315s | 0.69% |
| Solve: CG | 1 | 0.00186s | 4.1% |
| Solve: Preconditioner setup | 1 | 0.000968s | 2.1% |
+---------------------------------+-----------+------------+------------+
Cycle 1:
Number of active cells: 126 (3 global levels)
Workload imbalance: 1.17483
Number of degrees of freedom: 1672 (by
level: 117, 665, 1100)
Number of CG iterations: 11
Wrote solution_01.pvtu
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 0.0433s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble right hand side | 1 | 0.000286s | 0.66% |
| Estimate | 1 | 0.0272s | 63% |
| Output results | 1 | 0.00333s | 7.7% |
| Refine grid | 1 | 0.00196s | 4.5% |
| Setup | 1 | 0.0023s | 5.3% |
| Setup multigrid | 1 | 0.00262s | 6% |
| Solve | 1 | 0.00549s | 13% |
| Solve: 1 multigrid
V-cycle | 1 | 0.000343s | 0.79% |
| Solve: CG | 1 | 0.00293s | 6.8% |
| Solve: Preconditioner setup | 1 | 0.00174s | 4% |
+---------------------------------+-----------+------------+------------+
Cycle 2:
.
.
.
Here, the timing of the solve() function is split up in 3 parts: setting up the multigrid preconditioner, execution of a single multigrid V-cycle, and the CG solver. The V-cycle that is timed is unnecessary for the overall solve and only meant to give an insight at the different costs for AMG and GMG. Also it should be noted that when using the AMG solver, "Workload imbalance" is not included in the output since the hierarchy of coarse meshes is not required.
All results in this section are gathered on Intel Xeon Platinum 8280 (Cascade Lake) nodes which have 56 cores and 192GB per node and support AVX-512 instructions, allowing for vectorization over 8 doubles (vectorization used only in the matrix-free computations). The code is compiled using gcc 7.1.0 with intel-mpi 17.0.3. Trilinos 12.10.1 is used for the matrix-based GMG/AMG computations.
We can then gather a variety of information by calling the program with the input files that are provided in the directory in which step-50 is located. Using these, and adjusting the number of mesh refinement steps, we can produce information about how well the program scales.
The following table gives weak scaling timings for this program on up to 256M DoFs and 7,168 processors. (Recall that weak scaling keeps the number of degrees of freedom per processor constant while increasing the number of processors; i.e., it considers larger and larger problems.) Here, \(\mathbb{E}\) is the partition efficiency from the introduction (also equal to 1.0/workload imbalance), "Setup" is a combination of setup, setup multigrid, assemble, and assemble multigrid from the timing blocks, and "Prec" is the preconditioner setup. Ideally all times would stay constant over each problem size for the individual solvers, but since the partition efficiency decreases from 0.371 to 0.161 from largest to smallest problem size, we expect to see an approximately \(0.371/0.161=2.3\) times increase in timings for GMG. This is, in fact, pretty close to what we really get:
| | MF-GMG | | MB-GMG | | AMG |
| Procs | Cycle | DoFs | \(\mathbb{E}\) | | Setup | Prec | Solve | Total | | Setup | Prec | Solve | Total | | Setup | Prec | Solve | Total |
| 112 | 13 | 4M | 0.37 | | 0.742 | 0.393 | 0.200 | 1.335 | | 1.714 | 2.934 | 0.716 | 5.364 | | 1.544 | 0.456 | 1.150 | 3.150 |
| 448 | 15 | 16M | 0.29 | | 0.884 | 0.535 | 0.253 | 1.672 | | 1.927 | 3.776 | 1.190 | 6.893 | | 1.544 | 0.456 | 1.150 | 3.150 |
| 1,792 | 17 | 65M | 0.22 | | 1.122 | 0.686 | 0.309 | 2.117 | | 2.171 | 4.862 | 1.660 | 8.693 | | 1.654 | 0.546 | 1.460 | 3.660 |
| 7,168 | 19 | 256M | 0.16 | | 1.214 | 0.893 | 0.521 | 2.628 | | 2.386 | 7.260 | 2.560 | 12.206 | | 1.844 | 1.010 | 1.890 | 4.744 |
On the other hand, the algebraic multigrid in the last set of columns is relatively unaffected by the increasing imbalance of the mesh hierarchy (because it doesn't use the mesh hierarchy) and the growth in time is rather driven by other factors that are well documented in the literature (most notably that the algorithmic complexity of some parts of algebraic multigrid methods appears to be \({\cal O}(N \log N)\) instead of \({\cal O}(N)\) for geometric multigrid).
The upshort of the table above is that the matrix-free geometric multigrid method appears to be the fastest approach to solving this equation if not by a huge margin. Matrix-based methods, on the other hand, are consistently the worst.
The following figure provides strong scaling results for each method, i.e., we solve the same problem on more and more processors. Specifically, we consider the problems after 16 mesh refinement cycles (32M DoFs) and 19 cycles (256M DoFs), on between 56 to 28,672 processors:
While the matrix-based GMG solver and AMG scale similarly and have a similar time to solution (at least as long as there is a substantial number of unknowns per processor – say, several 10,000), the matrix-free GMG solver scales much better and solves the finer problem in roughly the same time as the AMG solver for the coarser mesh with only an eighth of the number of processors. Conversely, it can solve the same problem on the same number of processors in about one eighth the time.
Possible extensions
Testing convergence and higher order elements
The finite element degree is currently hard-coded as 2, see the template arguments of the main class. It is easy to change. To test, it would be interesting to switch to a test problem with a reference solution. This way, you can compare error rates.
Coarse solver
A more interesting example would involve a more complicated coarse mesh (see step-49 for inspiration). The issue in that case is that the coarsest level of the mesh hierarchy is actually quite large, and one would have to think about ways to solve the coarse level problem efficiently. (This is not an issue for algebraic multigrid methods because they would just continue to build coarser and coarser levels of the matrix, regardless of their geometric origin.)
In the program here, we simply solve the coarse level problem with a Conjugate Gradient method without any preconditioner. That is acceptable if the coarse problem is really small – for example, if the coarse mesh had a single cell, then the coarse mesh problems has a \(9\times 9\) matrix in 2d, and a \(27\times 27\) matrix in 3d; for the coarse mesh we use on the \(L\)-shaped domain of the current program, these sizes are \(21\times 21\) in 2d and \(117\times 117\) in 3d. But if the coarse mesh consists of hundreds or thousands of cells, this approach will no longer work and might start to dominate the overall run-time of each V-cyle. A common approach is then to solve the coarse mesh problem using an algebraic multigrid preconditioner; this would then, however, require assembling the coarse matrix (even for the matrix-free version) as input to the AMG implementation.
The plain program
#define FORCE_USE_OF_TRILINOS
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}
namespace ChangeVectorTypes
{
template <typename number>
const ::LinearAlgebra::distributed::Vector<number> &in)
{
out.locally_owned_elements());
#ifdef USE_PETSC_LA
ExcMessage(
"CopyVectorTypes::copy() not implemented for "
"PETSc vector types."));
#else
#endif
}
template <typename number>
{
#ifdef USE_PETSC_LA
(void)in;
ExcMessage(
"CopyVectorTypes::copy() not implemented for "
"PETSc vector types."));
#else
#endif
}
}
template <int dim>
class RightHandSide :
public Function<dim>
{
public:
const unsigned int = 0) const override
{
return 1.0;
}
template <typename number>
const unsigned int = 0) const
{
}
};
template <int dim>
class Coefficient :
public Function<dim>
{
public:
const unsigned int = 0) const override;
template <typename number>
const unsigned int = 0) const;
template <typename number>
template <typename number>
std::shared_ptr<Table<2, VectorizedArray<number>>> make_coefficient_table(
};
template <int dim>
{
for (
int d = 0;
d < dim; ++
d)
{
return 100.0;
}
return 1.0;
}
template <int dim>
template <typename number>
const unsigned int) const
{
for (unsigned int i = 0; i < VectorizedArray<number>::size(); ++i)
{
for (
int d = 0;
d < dim; ++
d)
{
return_value[i] = 100.0;
break;
}
}
return return_value;
}
template <int dim>
template <typename number>
number Coefficient<dim>::average_value(
{
number average(0);
for (unsigned int i = 0; i < points.size(); ++i)
average +=
value(points[i]);
average /= points.size();
return average;
}
template <int dim>
template <typename number>
std::shared_ptr<Table<2, VectorizedArray<number>>>
Coefficient<dim>::make_coefficient_table(
{
auto coefficient_table =
std::make_shared<Table<2, VectorizedArray<number>>>();
const unsigned int n_cells = mf_storage.n_macro_cells();
const unsigned int n_q_points = fe_eval.n_q_points;
coefficient_table->reinit(
n_cells, 1);
for (
unsigned int cell = 0; cell <
n_cells; ++cell)
{
fe_eval.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++q)
average_value +=
value(fe_eval.quadrature_point(q));
average_value /= n_q_points;
(*coefficient_table)(cell, 0) = average_value;
}
return coefficient_table;
}
{
bool try_parse(const std::string &prm_filename);
{
gmg_mb,
gmg_mf,
amg
};
int dimension;
double smoother_dampen;
unsigned int smoother_steps;
unsigned int n_steps;
bool output;
};
bool Settings::try_parse(const std::string &prm_filename)
{
"10",
"Number of adaptive refinement steps.");
"1.0",
"Dampen factor for the smoother.");
"1",
"Number of smoother steps.");
"MF",
"Switch between matrix-free GMG, "
"matrix-based GMG, and AMG.");
"false",
"Output graphical results.");
if (prm_filename.size() == 0)
{
std::cout << "**** Error: No input file provided!\n"
<< "**** Error: Call this program as './step-50 input.prm\n"
<< "\n"
<< "**** You may want to use one of the input files in this\n"
<< "**** directory, or use the following default values\n"
<< "**** to create an input file:\n";
return false;
}
try
{
}
catch (std::exception &
e)
{
std::cerr <<
e.what() << std::endl;
return false;
}
if (prm.
get(
"solver") ==
"MF")
this->solver = gmg_mf;
else if (prm.
get(
"solver") ==
"MB")
this->solver = gmg_mb;
else if (prm.
get(
"solver") ==
"AMG")
this->solver = amg;
else
this->smoother_dampen = prm.
get_double(
"smoother dampen");
this->smoother_steps = prm.
get_integer(
"smoother steps");
return true;
}
template <int dim, int degree>
class LaplaceProblem
{
public:
LaplaceProblem(
const Settings &settings);
private:
dim,
degree,
degree + 1,
1,
dim,
degree,
degree + 1,
1,
void setup_system();
void setup_multigrid();
void assemble_system();
void assemble_multigrid();
void assemble_rhs();
void solve();
void estimate();
void refine_grid();
void output_results(const unsigned int cycle);
MatrixType system_matrix;
MatrixFreeActiveMatrix mf_system_matrix;
};
template <int dim, int degree>
LaplaceProblem<dim, degree>::LaplaceProblem(
const Settings &settings)
: settings(settings)
, mpi_communicator(MPI_COMM_WORLD)
, mapping()
, fe(degree)
{
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::setup_system()
{
solution.
reinit(locally_owned_dofs, mpi_communicator);
right_hand_side.reinit(locally_owned_dofs, mpi_communicator);
constraints.reinit(locally_relevant_dofs);
constraints.close();
switch (settings.solver)
{
case Settings::gmg_mf:
{
std::shared_ptr<MatrixFree<dim, double>> mf_storage =
std::make_shared<MatrixFree<dim, double>>();
mf_storage->reinit(dof_handler,
constraints,
additional_data);
mf_system_matrix.initialize(mf_storage);
const Coefficient<dim> coefficient;
mf_system_matrix.set_coefficient(
coefficient.make_coefficient_table(*mf_storage));
break;
}
case Settings::gmg_mb:
case Settings::amg:
{
#ifdef USE_PETSC_LA
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
#else
locally_owned_dofs,
locally_relevant_dofs,
MPI_COMM_WORLD);
system_matrix.reinit(dsp);
#endif
break;
}
default:
}
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::setup_multigrid()
{
mg_constrained_dofs.clear();
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, boundary_ids);
switch (settings.solver)
{
case Settings::gmg_mf:
{
mf_mg_matrix.resize(0, n_levels - 1);
{
relevant_dofs);
level_constraints.
reinit(relevant_dofs);
mg_constrained_dofs.get_boundary_indices(
level));
level_constraints.
close();
std::shared_ptr<MatrixFree<dim, float>> mf_storage_level(
mf_storage_level->reinit(dof_handler,
level_constraints,
additional_data);
mf_mg_matrix[
level].initialize(mf_storage_level,
mg_constrained_dofs,
const Coefficient<dim> coefficient;
mf_mg_matrix[
level].set_coefficient(
coefficient.make_coefficient_table(*mf_storage_level));
mf_mg_matrix[
level].compute_diagonal();
}
break;
}
case Settings::gmg_mb:
{
mg_matrix.resize(0, n_levels - 1);
mg_matrix.clear_elements();
mg_interface_in.resize(0, n_levels - 1);
mg_interface_in.clear_elements();
{
dof_set);
{
#ifdef USE_PETSC_LA
dsp,
mpi_communicator,
dof_set);
dsp,
mpi_communicator);
#else
dof_set,
mpi_communicator);
mg_matrix[
level].reinit(dsp);
#endif
}
{
#ifdef USE_PETSC_LA
mg_constrained_dofs,
dsp,
dsp,
mpi_communicator,
dof_set);
mg_interface_in[
level].reinit(
dsp,
mpi_communicator);
#else
dof_set,
mpi_communicator);
mg_constrained_dofs,
dsp,
mg_interface_in[
level].reinit(dsp);
#endif
}
}
break;
}
default:
}
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_system()
{
quadrature_formula,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.
size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const Coefficient<dim> coefficient;
RightHandSide<dim> rhs;
std::vector<double> rhs_values(n_q_points);
if (cell->is_locally_owned())
{
cell_rhs = 0;
const double coefficient_value =
coefficient.average_value(fe_values.get_quadrature_points());
rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
coefficient_value *
fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) *
fe_values.JxW(q_point);
cell_rhs(i) +=
fe_values.shape_value(i, q_point) *
rhs_values[q_point] *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
cell_rhs,
local_dof_indices,
system_matrix,
right_hand_side);
}
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_multigrid()
{
quadrature_formula,
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.
size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const Coefficient<dim> coefficient;
std::vector<AffineConstraints<double>> boundary_constraints(
{
dof_set);
boundary_constraints[
level].reinit(dof_set);
boundary_constraints[
level].add_lines(
mg_constrained_dofs.get_refinement_edge_indices(
level));
boundary_constraints[
level].add_lines(
mg_constrained_dofs.get_boundary_indices(
level));
boundary_constraints[
level].close();
}
if (cell->level_subdomain_id() ==
triangulation.locally_owned_subdomain())
{
const double coefficient_value =
coefficient.average_value(fe_values.get_quadrature_points());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
coefficient_value * fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point);
cell->get_mg_dof_indices(local_dof_indices);
boundary_constraints[cell->level()].distribute_local_to_global(
cell_matrix, local_dof_indices, mg_matrix[cell->level()]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
if (mg_constrained_dofs.is_interface_matrix_entry(
cell->level(), local_dof_indices[i], local_dof_indices[j]))
mg_interface_in[cell->level()].add(local_dof_indices[i],
local_dof_indices[j],
}
for (
unsigned int i = 0; i <
triangulation.n_global_levels(); ++i)
{
}
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::assemble_rhs()
{
MatrixFreeActiveVector solution_copy;
MatrixFreeActiveVector right_hand_side_copy;
mf_system_matrix.initialize_dof_vector(solution_copy);
mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
solution_copy = 0.;
constraints.distribute(solution_copy);
solution_copy.update_ghost_values();
right_hand_side_copy = 0;
*(mf_system_matrix.get_coefficient());
RightHandSide<dim> right_hand_side_function;
*mf_system_matrix.get_matrix_free());
for (unsigned int cell = 0;
cell < mf_system_matrix.get_matrix_free()->n_macro_cells();
++cell)
{
phi.reinit(cell);
phi.read_dof_values_plain(solution_copy);
phi.evaluate(false, true, false);
for (unsigned int q = 0; q < phi.n_q_points; ++q)
{
phi.submit_gradient(-1.0 *
(coefficient(cell, 0) * phi.get_gradient(q)),
q);
phi.submit_value(
right_hand_side_function.value(phi.quadrature_point(q)), q);
}
phi.integrate_scatter(true, true, right_hand_side_copy);
}
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::solve()
{
SolverControl solver_control(1000, 1.
e-10 * right_hand_side.l2_norm());
solver_control.enable_history_data();
solution = 0.;
switch (settings.solver)
{
case Settings::gmg_mf:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
mg_transfer.
build(dof_handler);
MatrixFreeLevelMatrix,
coarse_grid_solver(coarse_solver, mf_mg_matrix[0],
identity);
Smoother,
MatrixFreeLevelVector>
smoother;
typename Smoother::AdditionalData(
settings.smoother_dampen));
smoother.set_steps(settings.smoother_steps);
mg_interface_matrices;
mg_interface_matrices.resize(0,
triangulation.n_global_levels() - 1);
mg_interface_matrices[
level].initialize(mf_mg_matrix[
level]);
mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
MatrixFreeLevelVector,
preconditioner(dof_handler,
mg, mg_transfer);
MatrixFreeActiveVector solution_copy;
MatrixFreeActiveVector right_hand_side_copy;
mf_system_matrix.initialize_dof_vector(solution_copy);
mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
computing_timer.leave_subsection("Solve: Preconditioner setup");
{
"Solve: 1 multigrid V-cycle");
preconditioner.vmult(solution_copy, right_hand_side_copy);
}
solution_copy = 0.;
{
solver.solve(mf_system_matrix,
solution_copy,
right_hand_side_copy,
preconditioner);
}
solution_copy.update_ghost_values();
constraints.distribute(solution);
break;
}
case Settings::gmg_mb:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
mg_transfer.
build(dof_handler);
MatrixType,
coarse_grid_solver(coarse_solver, mg_matrix[0],
identity);
#ifdef USE_PETSC_LA
settings.smoother_dampen == 1.0,
"PETSc's PreconditionJacobi has no support for a damping parameter."));
#else
smoother.
initialize(mg_matrix, settings.smoother_dampen);
#endif
mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
mg.set_edge_matrices(mg_out, mg_in);
preconditioner(dof_handler,
mg, mg_transfer);
computing_timer.leave_subsection("Solve: Preconditioner setup");
{
"Solve: 1 multigrid V-cycle");
preconditioner.vmult(solution, right_hand_side);
}
solution = 0.;
{
solver.solve(system_matrix,
solution,
right_hand_side,
preconditioner);
}
constraints.distribute(solution);
break;
}
case Settings::amg:
{
computing_timer.enter_subsection("Solve: Preconditioner setup");
PreconditionAMG::AdditionalData Amg_data;
#ifdef USE_PETSC_LA
Amg_data.symmetric_operator = true;
#else
Amg_data.elliptic = true;
Amg_data.smoother_type = "Jacobi";
Amg_data.higher_order_elements = true;
Amg_data.smoother_sweeps = settings.smoother_steps;
Amg_data.aggregation_threshold = 0.02;
#endif
Amg_data.output_details = false;
preconditioner.initialize(system_matrix, Amg_data);
computing_timer.leave_subsection("Solve: Preconditioner setup");
{
"Solve: 1 multigrid V-cycle");
preconditioner.vmult(solution, right_hand_side);
}
solution = 0.;
{
solver.solve(system_matrix,
solution,
right_hand_side,
preconditioner);
}
constraints.distribute(solution);
break;
}
default:
}
pcout << " Number of CG iterations: " << solver_control.last_step()
<< std::endl;
}
template <int dim>
struct ScratchData
{
const unsigned int quadrature_degree,
: fe_values(mapping, fe,
QGauss<dim>(quadrature_degree), update_flags)
, fe_interface_values(mapping,
fe,
QGauss<dim - 1>(quadrature_degree),
interface_update_flags)
{}
ScratchData(const ScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
scratch_data.fe_values.get_update_flags())
, fe_interface_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_interface_values.get_quadrature(),
scratch_data.fe_interface_values.get_update_flags())
{}
};
struct CopyData
{
CopyData()
{}
CopyData(const CopyData &) = default;
struct FaceData
{
unsigned int cell_indices[2];
double values[2];
};
unsigned int cell_index;
std::vector<FaceData> face_data;
};
template <int dim, int degree>
void LaplaceProblem<dim, degree>::estimate()
{
temp_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
temp_solution = solution;
const Coefficient<dim> coefficient;
auto cell_worker = [&](
const Iterator & cell,
ScratchData<dim> &scratch_data,
CopyData & copy_data) {
RightHandSide<dim> rhs;
const double rhs_value = rhs.value(cell->center());
const double nu = coefficient.value(cell->center());
std::vector<Tensor<2, dim>> hessians(fe_values.n_quadrature_points);
fe_values.get_function_hessians(temp_solution, hessians);
copy_data.cell_index = cell->active_cell_index();
double residual_norm_square = 0.;
for (unsigned k = 0; k < fe_values.n_quadrature_points; ++k)
{
const double residual = (rhs_value + nu *
trace(hessians[k]));
residual_norm_square += residual * residual * fe_values.JxW(k);
}
copy_data.value = cell->diameter() *
std::sqrt(residual_norm_square);
};
auto face_worker = [&](
const Iterator & cell,
const unsigned int &f,
const unsigned int &sf,
const unsigned int &nf,
const unsigned int &nsf,
ScratchData<dim> & scratch_data,
CopyData & copy_data) {
scratch_data.fe_interface_values;
fe_interface_values.
reinit(cell, f, sf, ncell, nf, nsf);
copy_data.face_data.emplace_back();
CopyData::FaceData ©_data_face = copy_data.face_data.back();
copy_data_face.cell_indices[0] = cell->active_cell_index();
copy_data_face.cell_indices[1] = ncell->active_cell_index();
const double coeff1 = coefficient.value(cell->center());
const double coeff2 = coefficient.value(ncell->center());
std::vector<Tensor<1, dim>> grad_u[2];
for (unsigned int i = 0; i < 2; ++i)
{
temp_solution, grad_u[i]);
}
double jump_norm_square = 0.;
for (unsigned int qpoint = 0;
++qpoint)
{
const double jump =
coeff1 * grad_u[0][qpoint] * fe_interface_values.
normal(qpoint) -
coeff2 * grad_u[1][qpoint] * fe_interface_values.
normal(qpoint);
jump_norm_square += jump * jump * fe_interface_values.
JxW(qpoint);
}
const double h = cell->face(f)->measure();
copy_data_face.values[0] = 0.5 * h *
std::sqrt(jump_norm_square);
copy_data_face.values[1] = copy_data_face.values[0];
};
auto copier = [&](const CopyData ©_data) {
estimated_error_per_cell[copy_data.cell_index] += copy_data.value;
for (auto &cdf : copy_data.face_data)
for (unsigned int j = 0; j < 2; ++j)
estimated_error_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
const unsigned int n_gauss_points = degree + 1;
ScratchData<dim> scratch_data(mapping,
fe,
n_gauss_points,
CopyData copy_data;
cell_worker,
copier,
scratch_data,
copy_data,
nullptr,
face_worker);
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::refine_grid()
{
const double refinement_fraction = 1. / (
std::pow(2.0, dim) - 1.);
triangulation, estimated_error_per_cell, refinement_fraction, 0.0);
}
template <int dim, int degree>
void LaplaceProblem<dim, degree>::output_results(const unsigned int cycle)
{
temp_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
temp_solution = solution;
for (unsigned int i = 0; i < subdomain.size(); ++i)
level(cell->active_cell_index()) = cell->level();
if (estimated_error_per_cell.size() > 0)
"estimated_error_per_cell");
"", "solution", cycle, mpi_communicator, 2 , 1 );
pcout << " Wrote " << master << std::endl;
}
template <int dim, int degree>
{
for (unsigned int cycle = 0; cycle < settings.n_steps; ++cycle)
{
pcout << "Cycle " << cycle << ':' << std::endl;
if (cycle > 0)
refine_grid();
pcout << " Number of active cells: "
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
pcout <<
" (" <<
triangulation.n_global_levels() <<
" global levels)" << std::endl
<< " Partition efficiency: "
pcout << std::endl;
setup_system();
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
setup_multigrid();
pcout <<
" Number of degrees of freedom: " << dof_handler.
n_dofs();
if (settings.solver == Settings::gmg_mf ||
settings.solver == Settings::gmg_mb)
{
pcout << " (by level: ";
", ");
}
pcout << std::endl;
if (settings.solver == Settings::gmg_mf)
assemble_rhs();
else
{
assemble_system();
if (settings.solver == Settings::gmg_mb)
assemble_multigrid();
}
solve();
estimate();
if (settings.output)
output_results(cycle);
computing_timer.print_summary();
computing_timer.reset();
}
}
int main(int argc, char *argv[])
{
if (!settings.try_parse((argc > 1) ? (argv[1]) : ""))
return 0;
try
{
constexpr unsigned int fe_degree = 2;
switch (settings.dimension)
{
case 2:
{
LaplaceProblem<2, fe_degree> test(settings);
test.run();
break;
}
case 3:
{
LaplaceProblem<3, fe_degree> test(settings);
test.run();
break;
}
default:
}
}
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;
MPI_Abort(MPI_COMM_WORLD, 1);
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
MPI_Abort(MPI_COMM_WORLD, 2);
return 1;
}
return 0;
}
1.8.17