This program demonstrates how to use the cell-based implementation of finite
element operators with the MatrixFree class, first introduced in step-37, to
-solve nonlinear partial differential equations. Moreover, we demonstrate how
-the MatrixFree class handles constraints, an issue shortly mentioned in the
-results section of step-37. Finally, we will use an explicit time-stepping
+solve nonlinear partial differential equations. Moreover, we have another look
+at the handling of constraints within the matrix-free framework.
+Finally, we will use an explicit time-stepping
method to solve the problem and introduce Gauss-Lobatto finite elements that
are very convenient in this case since their mass matrix can be accurately
approximated by a diagonal, and thus trivially invertible, matrix. The two
Lagrange polynomials according to the point distribution of the Gauss-Lobatto
quadrature rule. Secondly, the quadrature is done with the same Gauss-Lobatto
quadrature rule. In this formula, the integrals $\int_K \varphi_i \varphi_j
-dx\approx \sum_q \varphi_i \varphi_j \mathrm{det}(J) \big |_{x_q}$ are
-approximated to zero whenever $i\neq j$, because on the points defining the
-Lagrange polynomials exactly one function $\varphi_j$ is one and all others
-zero. Moreover, the Gauss-Lobatto distribution of nodes of Lagrange
+dx\approx \sum_q \varphi_i \varphi_j \mathrm{det}(J) \big |_{x_q}$ become
+zero whenever $i\neq j$, because exactly one function $\varphi_j$ is one and
+all others zero in the points defining the Lagrange polynomials.
+Moreover, the Gauss-Lobatto distribution of nodes of Lagrange
polynomials clusters the nodes towards the element boundaries. This results in
a well-conditioned polynomial basis for high-order discretization
methods. Indeed, the condition number of an FE_Q elements with equidistant
nodes grows exponentially with the degree, which destroys any benefit for
orders of about five and higher. For this reason, Gauss-Lobatto points are the
-default distribution for FE_Q (but at degrees one and two, those are
-equivalent to the equidistant points).
+default distribution for the FE_Q element (but at degrees one and two, those
+are equivalent to the equidistant points).
<h3> Problem statement and discretization </h3>
that was already introduced in step-25. As a simple explicit time
integration method, we choose leap frog scheme using the second-order
-formulation of the equation. Then, the scheme reads in weak form
+formulation of the equation. With this time stepping, the scheme reads in
+weak form
\f{eqnarray*}
(v,u^{n+1}) = (v,2 u^n-u^{n-1} -
on the left hand side of the equation, making the solution of the
linear system in each time step trivial.
-Using this quadrature rule, for a <i>p</i>th order finite element, we
-use a <i>(2p-1)</i>th order accurate formula to evaluate the
-integrals. Since the product of two <i>p</i>th order basis functions
-when computing a mass matrix gives a function with polynomial degree
-<i>2p</i> in each direction, the integrals are not computed exactly.
-However, considering the fact that the interpolation order
-of finite elements of degree <i>p</i> is <i>p+1</i>, the overall
-convergence properties are not disturbed by the quadrature error, in
-particular not when we use high orders.
+Using this quadrature rule, for a <i>p</i>th order finite element, we use a
+<i>(2p-1)</i>th order accurate formula to evaluate the integrals. Since the
+product of two <i>p</i>th order basis functions when computing a mass matrix
+gives a function with polynomial degree <i>2p</i> in each direction, the
+integrals are not computed exactly. However, the overall convergence
+properties are not disturbed by the quadrature error on meshes with affine
+element shapes with L2 errors proportional to <i>h<sup>p+1</sup></i>. Note
+however that order reduction with sub-optimal convergence rates of the L2
+error of <i>O(h<sup>p</sup>)</i> or even <i>O(h<sup>p-1</sup>)</i> for some 3D
+setups has been reported <a href="https://dx.doi.org/10.1002/num.20353">in
+literature</a> on deformed (non-affine) element shapes for wave equations
+when the integrand is not a polynomial any more.
Apart from the fact that we avoid solving linear systems with this
type of elements when using explicit time-stepping, they come with two
Gauss-Lobatto elements, where quadrature points and node points of the
finite element coincide, this operation is trivial since the value
of the function at the quadrature points is given by its one-dimensional
-coefficients. In this way, the complexity of a finite element operator
-evaluation is further reduced compared to equidistant elements.
-
-The third advantage is the fact that these elements are better conditioned than
-equidistant Lagrange polynomials for increasing order so that we can use
-higher order elements for an accurate solution of the equation. Lagrange
-elements FE_Q with equidistant points should not be used for polynomial
-degrees four and higher.
+coefficients. In this way, the arithmetic work for the finite element operator
+evaluation is reduced by approximately a factor of two compared to the generic
+Gaussian quadrature.
To sum up the discussion, by using the right finite element and
-quadrature rule combination, we end up with a scheme where we in each
-time step need to compute the right hand side vector corresponding
+quadrature rule combination, we end up with a scheme where we
+only need to compute the right hand side vector corresponding
to the formulation above and then multiply it by the inverse of the
-diagonal mass matrix. In practice, of course, we extract the diagonal
-elements and invert them only once at the beginning of the program.
+diagonal mass matrix in each time step. In practice, of course, we extract
+the diagonal elements and invert them only once at the beginning of the
+program.
<h3>Implementation of constraints</h3>
how they are constrained. This format uses an unnecessarily large
amount of memory since there are not so many different types of
constraints: for example, in the case of hanging nodes when using
-linear finite element on every cell, constraints most have the form
+linear finite element on every cell, most constraints have the form
$x_k = \frac 12 x_i + \frac 12 x_j$ where the coefficients $\frac 12$
are always the same and only $i,j,k$ are different. While storing this
redundant information is not a problem in general because it is only
needed once during matrix and right hand side assembly, it becomes a
-problem when we want to use the matrix-free approach since there this
-information has to be accessed every time we apply the operator. Thus,
-instead of an AffineConstraints object, we use a variable that we call
-<code>constraint_pool</code> that collects the weights of the
-different constraints. Then, we only have to store an identifier of
-each constraint in the mesh instead of all the weights. Moreover, we
-do not want to apply the constraints in a pre- and postprocessing step
-but want to take care of constraints as we evaluate the finite element
-operator. Therefore, we embed the constraint information into the
+bottleneck in the matrix-free approach since there this
+information has to be accessed every time we apply the operator, and the
+remaining components of the operator evaluation are so fast. Thus,
+instead of an AffineConstraints object, MatrixFree uses a variable that
+we call <code>constraint_pool</code> that collects the weights of the
+different constraints. Then, only an identifier of each constraint in the
+mesh instead of all the weights have to be stored. Moreover,
+the constraints are not applied in a pre- and postprocessing step
+but rather as we evaluate the finite element
+operator. Therefore, the constraint information is embedded into the
variable <code>indices_local_to_global</code> that is used to extract
the cell information from the global vector. If a DoF is constrained,
the <code>indices_local_to_global</code> variable contains the global
indices of the DoFs that it is constrained to. Then, we have another
variable <code>constraint_indicator</code> at hand that holds, for
each cell, the local indices of DoFs that are constrained as well as
-the identifier of the type of constraint. Actually, you will not see
+the identifier of the type of constraint. Fortunately, you will not see
these data structures in the example program since the class
<code>FEEvaluation</code> takes care of the constraints without user
interaction.
+In the presence of hanging nodes, the diagonal mass matrix obtained on the
+element level via the Gauss-Lobatto quadrature/node point procedure does not
+directly translate to a diagonal global mass matrix, as following the
+constraints on rows and columns would also add off-diagonal entries. As
+explained in <a href="https://dx.doi.org/10.4208/cicp.101214.021015a">Kormann
+(2016)</a>, interpolating constraints on a vector, which maintains the
+diagonal shape of the mass matrix, is consistent with the equations up to an
+error of the same magnitude as the quadrature error. In the program below, we
+will simply assemble the diagonal of the mass matrix as if it were a vector to
+enable this approximation.
+
<h3> Parallelization </h3>
The MatrixFree class comes with the option to be parallelized on three levels:
MPI parallelization on clusters of distributed nodes, thread parallelization
scheduled by the Threading Building Blocks library, and finally with a
-vectorization by clustering of two (or more) cells into a SIMD data type for
-the operator application. As we have already discussed in step-37, you will
+vectorization by working on a batch of two (or more) cells via SIMD data type
+(sometimes called cross-element or external vectorization).
+As we have already discussed in step-37, you will
get best performance by using an instruction set specific to your system,
e.g. with the cmake variable <tt>-DCMAKE_CXX_FLAGS="-march=native"</tt>. The
MPI parallelization was already exploited in step-37. Here, we additionally
need to do is to tell the initialization of the MatrixFree object about the
fact that we want to use a thread parallel scheme through the variable
MatrixFree::AdditionalData::thread_parallel_scheme. During setup, a dependency
-graph similar to the one described in the @ref workstream_paper , which allows
-the code to schedule the work of the @p local_apply function on chunks of cell
-without several threads accessing the same vector indices. As opposed to the
-WorkStream loops, some additional clever tricks to avoid global
+graph is set up similar to the one described in the @ref workstream_paper ,
+which allows to schedule the work of the @p local_apply function on chunks of
+cells without several threads accessing the same vector indices. As opposed to
+the WorkStream loops, some additional clever tricks to avoid global
synchronizations as described in <a
href="https://dx.doi.org/10.1109/eScience.2011.53">Kormann and Kronbichler
(2011)</a> are also applied.
It is apparent that the matrix-free code outperforms the standard assembly
routines in deal.II by far. In 3D and for fourth order elements, one operator
-application is also almost ten times as fast as a sparse matrix-vector
+evaluation is also almost ten times as fast as a sparse matrix-vector
product.
-<h3>Parallel run in 3D</h3>
+<h3>Parallel run in 2D and 3D</h3>
-To demonstrate how the example scales for a parallel run and to demonstrate
-that hanging node constraints can be handled in an efficient way, we run the
-example in 3D with $\mathcal{Q}_4$ elements. First, we run it on a notebook
-with 2 cores (Sandy Bridge CPU) at 2.7 GHz.
+We start with the program output obtained on a workstation with 12 cores / 24
+threads (one Intel Xeon E5-2687W v4 CPU running at 3.2 GHz, hyperthreading
+enabled), running the program in release mode:
@code
-\$ make debug-mode=off run
+\$ make run
+Number of MPI ranks: 1
+Number of threads on each rank: 24
+Vectorization over 4 doubles = 256 bits (AVX), VECTORIZATION_LEVEL=2
+
+ Number of global active cells: 15412
+ Number of degrees of freedom: 249065
+ Time step size: 0.00292997, finest cell: 0.117188
+
+ Time: -10, solution norm: 9.5599
+ Time: -9.41, solution norm: 17.678
+ Time: -8.83, solution norm: 23.504
+ Time: -8.24, solution norm: 27.5
+ Time: -7.66, solution norm: 29.513
+ Time: -7.07, solution norm: 29.364
+ Time: -6.48, solution norm: 27.23
+ Time: -5.9, solution norm: 23.527
+ Time: -5.31, solution norm: 18.439
+ Time: -4.73, solution norm: 11.935
+ Time: -4.14, solution norm: 5.5284
+ Time: -3.55, solution norm: 8.0354
+ Time: -2.97, solution norm: 14.707
+ Time: -2.38, solution norm: 20
+ Time: -1.8, solution norm: 22.834
+ Time: -1.21, solution norm: 22.771
+ Time: -0.624, solution norm: 20.488
+ Time: -0.0381, solution norm: 16.697
+ Time: 0.548, solution norm: 11.221
+ Time: 1.13, solution norm: 5.3912
+ Time: 1.72, solution norm: 8.4528
+ Time: 2.31, solution norm: 14.335
+ Time: 2.89, solution norm: 18.555
+ Time: 3.48, solution norm: 20.894
+ Time: 4.06, solution norm: 21.305
+ Time: 4.65, solution norm: 19.903
+ Time: 5.24, solution norm: 16.864
+ Time: 5.82, solution norm: 12.223
+ Time: 6.41, solution norm: 6.758
+ Time: 6.99, solution norm: 7.2423
+ Time: 7.58, solution norm: 12.888
+ Time: 8.17, solution norm: 17.273
+ Time: 8.75, solution norm: 19.654
+ Time: 9.34, solution norm: 19.838
+ Time: 9.92, solution norm: 17.964
+ Time: 10, solution norm: 17.595
+
+ Performed 6826 time steps.
+ Average wallclock time per time step: 0.0013453s
+ Spent 14.976s on output and 9.1831s on computations.
+@endcode
+
+In 3D, the respective output looks like
+@code
+\$ make run
+Number of MPI ranks: 1
+Number of threads on each rank: 24
+Vectorization over 4 doubles = 256 bits (AVX), VECTORIZATION_LEVEL=2
+
Number of global active cells: 17592
Number of degrees of freedom: 1193881
Time step size: 0.0117233, finest cell: 0.46875
Time: 10, solution norm: 94.115
Performed 1706 time steps.
- Average wallclock time per time step: 0.038261s
- Spent 11.977s on output and 65.273s on computations.
+ Average wallclock time per time step: 0.0084542s
+ Spent 16.766s on output and 14.423s on computations.
@endcode
-It takes 0.04 seconds for one time step on a notebook with more than a million
+It takes 0.008 seconds for one time step with more than a million
degrees of freedom (note that we would need many processors to reach such
-numbers when solving linear systems). If we run the same 3D code on a
-cluster with 2 nodes and each node runs 8 threads, we get the following times:
+numbers when solving linear systems).
+If we replace the thread-parallelization by a pure MPI parallelization, the
+timings change into:
@code
-\$ mpirun --bynode -n 2 ./\step-48
+\$ mpirun -n 24 ./step-48
+Number of MPI ranks: 24
+Number of threads on each rank: 1
+Vectorization over 4 doubles = 256 bits (AVX), VECTORIZATION_LEVEL=2
...
Performed 1706 time steps.
- Average wallclock time per time step: 0.0123188s
- Spent 6.74378s on output and 21.0158s on computations.
+ Average wallclock time per time step: 0.0051747s
+ Spent 2.0535s on output and 8.828s on computations.
@endcode
-We observe a considerable speedup over the notebook (16 cores versus 2 cores;
-nonetheless, one notebook core is considerably faster than one core of the
-cluster because of a newer processor architecture). If we run the same program
-on 4 nodes with 8 threads on each node, we get:
+We observe a dramatic speedup for the output (which makes sense, given that
+most code of the output is not parallelized via threads, whereas it is for
+MPI), but less than the theoretical factor of 12 we would expect from the
+parallelism. More interestingly, the computations also get faster when
+switching from the threads-only variant to the MPI-only variant. This is a
+general observation for the MatrixFree framework (as of updating this data in
+2019). The main reason is that the decisions regarding work on conflicting
+cell batches made to enable execution in parallel are overly pessimistic:
+While they ensure that no work on neighboring cells is done on different
+threads at the same time, this conservative setting implies that data from
+neighboring cells is also evicted from caches by the time neighbors get
+touched. Furthermore, the current scheme is not able to provide a constant
+load for all 24 threads for the given mesh with 17,592 cells.
+
+The current program allows to also mix MPI parallelization with thread
+parallelization. This is most beneficial when running programs on clusters
+with multiple nodes, using MPI for the inter-node parallelization and threads
+for the intra-node parallelization. On the workstation used above, we can run
+threads in the hyperthreading region (i.e., using 2 threads for each of the 12
+MPI ranks). An important setting for mixing MPI with threads is to ensure
+proper binning of tasks to CPUs. On many clusters the placing is either
+automatically via the `mpirun/mpiexec` environment, or there can be manual
+settings. Here, we simply report the run times the plain version of the
+program (noting that things could be improved towards the timings of the
+MPI-only program when proper pinning is done):
@code
-\$ mpirun --bynode -n 4 ./\step-48
+\$ mpirun -n 12 ./step-48
+Number of MPI ranks: 12
+Number of threads on each rank: 2
+Vectorization over 4 doubles = 256 bits (AVX), VECTORIZATION_LEVEL=2
...
Performed 1706 time steps.
- Average wallclock time per time step: 0.00689865s
- Spent 3.54145s on output and 11.7691s on computations.
+ Average wallclock time per time step: 0.0056651s
+ Spent 2.5175s on output and 9.6646s on computations.
@endcode
-By comparing the times for two nodes and four nodes, we observe the nice
-scaling behavior of the implementation. Of course, the code can also be run in
-MPI-mode only by disabling the multithreading flag in the code. If we use the
-same 32 cores as for the hybrid parallelization above, we observe the
-following run-time:
-
-@code
-\$ mpirun -n 32 ./\step-48
-...
- Performed 1706 time steps.
- Average wallclock time per time step: 0.0189041s
- Spent 0.968967s on output and 32.2504s on computations.
-@endcode
-We observe slower speed for computations, but faster output (which makes
-sense, as output is only parallelized by MPI and not threads), whereas the
-computations are faster if we use hybrid parallelism in the given case.
<h3>Possibilities for extensions</h3>
// mass matrix. Since we use Gauss-Lobatto elements, the mass matrix is a
// diagonal matrix and can be stored as a vector. The computation of the
// mass matrix diagonal is simple to achieve with the data structures
- // provided by FEEvaluation: Just loop over all (macro-) cells and integrate
- // over the function that is constant one on all quadrature points by using
- // the <code>integrate</code> function with @p true argument at the slot for
- // values. Finally, we invert the diagonal entries since we have to multiply
- // by the inverse mass matrix in each time step.
+ // provided by FEEvaluation: Just loop over all cell batches, i.e.,
+ // collections of cells due to SIMD vectorization, and integrate over the
+ // function that is constant one on all quadrature points by using the
+ // <code>integrate</code> function with @p true argument at the slot for
+ // values. Finally, we invert the diagonal entries to have the inverse mass
+ // matrix directly available in each time step.
template <int dim, int fe_degree>
SineGordonOperation<dim, fe_degree>::SineGordonOperation(
const MatrixFree<dim, double> &data_in,
: data(data_in)
, delta_t_sqr(make_vectorized_array(time_step * time_step))
{
- VectorizedArray<double> one = make_vectorized_array(1.);
-
data.initialize_dof_vector(inv_mass_matrix);
FEEvaluation<dim, fe_degree> fe_eval(data);
const unsigned int n_q_points = fe_eval.n_q_points;
- for (unsigned int cell = 0; cell < data.n_macro_cells(); ++cell)
+ for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
{
fe_eval.reinit(cell);
for (unsigned int q = 0; q < n_q_points; ++q)
- fe_eval.submit_value(one, q);
+ fe_eval.submit_value(make_vectorized_array(1.), q);
fe_eval.integrate(true, false);
fe_eval.distribute_local_to_global(inv_mass_matrix);
}
inv_mass_matrix.local_element(k) =
1. / inv_mass_matrix.local_element(k);
else
- inv_mass_matrix.local_element(k) = 0;
+ inv_mass_matrix.local_element(k) = 1;
}
// of shape function values on quadrature points which is simply the
// injection of the values of cell degrees of freedom. The MatrixFree class
// detects possible structure of the finite element at quadrature points
- // when initializing, which is then used by FEEvaluation for selecting the
- // most appropriate numerical kernel.
+ // when initializing, which is then automatically used by FEEvaluation for
+ // selecting the most appropriate numerical kernel.
// The nonlinear function that we have to evaluate for the time stepping
// routine includes the value of the function at the present time @p current
// as well as the value at the previous time step @p old. Both values are
- // passed to the operator in the collection of source vectors @p src, which is
- // simply a <tt>std::vector</tt> of pointers to the actual solution
+ // passed to the operator in the collection of source vectors @p src, which
+ // is simply a <tt>std::vector</tt> of pointers to the actual solution
// vectors. This construct of collecting several source vectors into one is
// necessary as the cell loop in @p MatrixFree takes exactly one source and
- // one destination vector, even if we happen to use many vectors like the two
- // in this case. Note that the cell loop accepts any valid class for input and
- // output, which does not only include vectors but general data types.
- // However, only in case it encounters a
+ // one destination vector, even if we happen to use many vectors like the
+ // two in this case. Note that the cell loop accepts any valid class for
+ // input and output, which does not only include vectors but general data
+ // types. However, only in case it encounters a
// LinearAlgebra::distributed::Vector<Number> or a <tt>std::vector</tt>
- // collecting these vectors, it calls functions that exchange data at the
- // beginning and the end of the loop. In the loop over the cells, we first
- // have to read in the values in the vectors related to the local values.
- // Then, we evaluate the value and the gradient of the current solution vector
- // and the values of the old vector at the quadrature points. Then, we combine
- // the terms in the scheme in the loop over the quadrature points. Finally, we
- // integrate the result against the test
- // function and accumulate the result to the global solution vector @p dst.
+ // collecting these vectors, it calls functions that exchange ghost data due
+ // to MPI at the beginning and the end of the loop. In the loop over the
+ // cells, we first have to read in the values in the vectors related to the
+ // local values. Then, we evaluate the value and the gradient of the
+ // current solution vector and the values of the old vector at the
+ // quadrature points. Next, we combine the terms in the scheme in the loop
+ // over the quadrature points. Finally, we integrate the result against the
+ // test function and accumulate the result to the global solution vector @p
+ // dst.
template <int dim, int fe_degree>
void SineGordonOperation<dim, fe_degree>::local_apply(
const MatrixFree<dim> & data,
//@sect4{SineGordonOperation::apply}
// This function performs the time stepping routine based on the cell-local
- // strategy. First the destination vector is set to zero, then the cell-loop
- // is called, and finally the solution is multiplied by the inverse mass
- // matrix. The structure of the cell loop is implemented in the cell finite
+ // strategy. Note that we need to set the destination vector to zero before
+ // we add the integral contributions of the current time step (via the
+ // FEEvaluation::distribute_local_to_global() call). In this tutorial, we
+ // let the cell-loop do the zero operation via the fifth `true` argument
+ // passed to MatrixFree::cell_loop. The loop can schedule the zero operation
+ // closer to the operations on vector entries for supported vector entries,
+ // thereby possibly increasing data locality (the vector entries that first
+ // get zeroed are later re-used in the `distribute_local_to_global()`
+ // call). The structure of the cell loop is implemented in the cell finite
// element operator class. On each cell it applies the routine defined as
// the <code>local_apply()</code> method of the class
// <code>SineGordonOperation</code>, i.e., <code>this</code>. One could also
- // provide a function with the same signature that is not part of a class.
+ // provide a function with the same signature that is not part of a
+ // class. Finally, the result of the integration is multiplied by the
+ // inverse mass matrix.
template <int dim, int fe_degree>
void SineGordonOperation<dim, fe_degree>::apply(
LinearAlgebra::distributed::Vector<double> & dst,
const std::vector<LinearAlgebra::distributed::Vector<double> *> &src) const
{
- dst = 0;
- data.cell_loop(&SineGordonOperation<dim, fe_degree>::local_apply,
- this,
- dst,
- src);
+ data.cell_loop(
+ &SineGordonOperation<dim, fe_degree>::local_apply, this, dst, src, true);
dst.scale(inv_mass_matrix);
}
+
//@sect3{Equation data}
// We define a time-dependent function that is used as initial
// value. Different solutions can be obtained by varying the starting
- // time. This function has already been explained in step-25.
+ // time. This function, taken from step-25, would represent an analytic
+ // solution in 1D for all times, but is merely used for setting some
+ // starting solution of interest here. More elaborate choices that could
+ // test the convergence of this program are given in step-25.
template <int dim>
- class ExactSolution : public Function<dim>
+ class InitialCondition : public Function<dim>
{
public:
- ExactSolution(const unsigned int n_components = 1, const double time = 0.)
+ InitialCondition(const unsigned int n_components = 1,
+ const double time = 0.)
: Function<dim>(n_components, time)
{}
- virtual double value(const Point<dim> & p,
- const unsigned int component = 0) const override;
+ virtual double value(const Point<dim> &p,
+ const unsigned int /*component*/) const override
+ {
+ double t = this->get_time();
+
+ const double m = 0.5;
+ const double c1 = 0.;
+ const double c2 = 0.;
+ const double factor =
+ (m / std::sqrt(1. - m * m) * std::sin(std::sqrt(1. - m * m) * t + c2));
+ double result = 1.;
+ for (unsigned int d = 0; d < dim; ++d)
+ result *= -4. * std::atan(factor / std::cosh(m * p[d] + c1));
+ return result;
+ }
};
- template <int dim>
- double ExactSolution<dim>::value(const Point<dim> &p,
- const unsigned int /* component */) const
- {
- double t = this->get_time();
-
- const double m = 0.5;
- const double c1 = 0.;
- const double c2 = 0.;
- const double factor =
- (m / std::sqrt(1. - m * m) * std::sin(std::sqrt(1. - m * m) * t + c2));
- double result = 1.;
- for (unsigned int d = 0; d < dim; ++d)
- result *= -4. * std::atan(factor / std::cosh(m * p[d] + c1));
- return result;
- }
-
// @sect3{SineGordonProblem class}
// though), see also the discussion in the introduction. Note that FE_Q
// selects the Gauss-Lobatto nodal points by default due to their improved
// conditioning versus equidistant points. To make things more explicit, we
- // choose to state the selection of the nodal points nonetheless.
+ // state the selection of the nodal points nonetheless.
template <int dim>
SineGordonProblem<dim>::SineGordonProblem()
: pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
, n_global_refinements(10 - 2 * dim)
, time(-10)
, time_step(10.)
- , final_time(10)
+ , final_time(10.)
, cfl_number(.1 / fe_degree)
, output_timestep_skip(200)
{}
// solution. As in step-40, we need to equip the constraint matrix with
// the IndexSet of locally relevant degrees of freedom to avoid it to
// consume too much memory for big problems. Next, the <code> MatrixFree
- // </code> for the problem is set up. Note that we specify the MPI
- // communicator which we are going to use, and that we also want to use
- // shared-memory parallelization (hence one would use multithreading for
- // intra-node parallelism and not MPI; note that we here choose the
- // standard option — if we wanted to disable shared memory
- // parallelization, we would choose @p none). Finally, three solution
- // vectors are initialized. MatrixFree stores the layout that is to be
- // used by distributed vectors, so we just ask it to initialize the
- // vectors.
+ // </code> object for the problem is set up. Note that we specify a
+ // particular scheme for shared-memory parallelization (hence one would
+ // use multithreading for intra-node parallelism and not MPI; we here
+ // choose the standard option — if we wanted to disable shared
+ // memory parallelization even in case where there is more than one TBB
+ // thread available in the program, we would choose
+ // MatrixFree::AdditionalData::TasksParallelScheme::none). Also note that,
+ // instead of using the default QGauss quadrature argument, we supply a
+ // QGaussLobatto quadrature formula to enable the desired
+ // behavior. Finally, three solution vectors are initialized. MatrixFree
+ // expects a particular layout of ghost indices (as it handles index
+ // access in MPI-local numbers that need to match between the vector and
+ // MatrixFree), so we just ask it to initialize the vectors to be sure the
+ // ghost exchange is properly handled.
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);
constraints.close();
- QGaussLobatto<1> quadrature(fe_degree + 1);
typename MatrixFree<dim>::AdditionalData additional_data;
additional_data.tasks_parallel_scheme =
- MatrixFree<dim>::AdditionalData::partition_partition;
+ MatrixFree<dim>::AdditionalData::TasksParallelScheme::partition_partition;
matrix_free_data.reinit(dof_handler,
constraints,
- quadrature,
+ QGaussLobatto<1>(fe_degree + 1),
additional_data);
matrix_free_data.initialize_dof_vector(solution);
// This function prints the norm of the solution and writes the solution
// vector to a file. The norm is standard (except for the fact that we need
- // to accumulate the norms over all processors for the parallel grid), and
- // the second is similar to what we did in step-40 or step-37. Note that we
- // can use the same vector for output as the one used during computations:
- // The vectors in the matrix-free framework always provide full information
- // on all locally owned cells (this is what is needed in the local
- // evaluations, too), including ghost vector entries on these cells. This is
- // the only data that is needed in the integrate_difference function as well
- // as in DataOut. The only action to take at this point is then to make sure
- // that the vector updates its ghost values before we read from them. This
- // is a feature present only in the LinearAlgebra::distributed::Vector
- // class. Distributed vectors with PETSc and Trilinos, on the other hand,
- // need to be copied to special vectors including ghost values (see the
- // relevant section in step-40). If we also wanted to access all degrees of
- // freedom on ghost cells (e.g. when computing error estimators that use the
- // jump of solution over cell boundaries), we would need more information
- // and create a vector initialized with locally relevant dofs just as in
- // step-40. Observe also that we need to distribute constraints for output -
- // they are not filled during computations (rather, they are interpolated on
- // the fly in the matrix-free method read_dof_values).
+ // to accumulate the norms over all processors for the parallel grid which
+ // we do via the VectorTools::compute_global_error() function), and the
+ // second is similar to what we did in step-40 or step-37. Note that we can
+ // use the same vector for output as the one used during computations: The
+ // vectors in the matrix-free framework always provide full information on
+ // all locally owned cells (this is what is needed in the local evaluations,
+ // too), including ghost vector entries on these cells. This is the only
+ // data that is needed in the VectorTools::integrate_difference() function
+ // as well as in DataOut. The only action to take at this point is to make
+ // sure that the vector updates its ghost values before we read from
+ // them. This is a feature present only in the
+ // LinearAlgebra::distributed::Vector class. Distributed vectors with PETSc
+ // and Trilinos, on the other hand, need to be copied to special vectors
+ // including ghost values (see the relevant section in step-40). If we also
+ // wanted to access all degrees of freedom on ghost cells (e.g. when
+ // computing error estimators that use the jump of solution over cell
+ // boundaries), we would need more information and create a vector
+ // initialized with locally relevant dofs just as in step-40. Observe also
+ // that we need to distribute constraints for output - they are not filled
+ // during computations (rather, they are interpolated on the fly in the
+ // matrix-free method FEEvaluation::read_dof_values()).
template <int dim>
void
SineGordonProblem<dim>::output_results(const unsigned int timestep_number)
// @sect4{SineGordonProblem::run}
- // This function is called by the main function and calls the subroutines of
- // the class.
+ // This function is called by the main function and steps into the
+ // subroutines of the class.
//
- // The first step is to set up the grid and the cell operator. Then, the
- // time step is computed from the CFL number given in the constructor and
- // the finest mesh size. The finest mesh size is computed as the diameter of
- // the last cell in the triangulation, which is the last cell on the finest
- // level of the mesh. This is only possible for Cartesian meshes, otherwise,
- // one needs to loop over all cells. Note that we need to query all the
- // processors for their finest cell since the not all processors might hold
- // a region where the mesh is at the finest level. Then, we readjust the
- // time step a little to hit the final time exactly.
+ // After printing some information about the parallel setup, the first
+ // action is to set up the grid and the cell operator. Then, the time step
+ // is computed from the CFL number given in the constructor and the finest
+ // mesh size. The finest mesh size is computed as the diameter of the last
+ // cell in the triangulation, which is the last cell on the finest level of
+ // the mesh. This is only possible for meshes where all elements on a level
+ // have the same size, otherwise, one needs to loop over all cells. Note
+ // that we need to query all the processors for their finest cell since
+ // not all processors might hold a region where the mesh is at the finest
+ // level. Then, we readjust the time step a little to hit the final time
+ // exactly.
template <int dim>
void SineGordonProblem<dim>::run()
{
+ {
+ pcout << "Number of MPI ranks: "
+ << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) << std::endl;
+ pcout << "Number of threads on each rank: "
+ << MultithreadInfo::n_threads() << std::endl;
+ const unsigned int n_vect_doubles =
+ VectorizedArray<double>::n_array_elements;
+ const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles;
+ pcout << "Vectorization over " << n_vect_doubles
+ << " doubles = " << n_vect_bits << " bits ("
+ << Utilities::System::get_current_vectorization_level()
+ << "), VECTORIZATION_LEVEL=" << DEAL_II_COMPILER_VECTORIZATION_LEVEL
+ << std::endl
+ << std::endl;
+ }
make_grid_and_dofs();
const double local_min_cell_diameter =
// difficulty and just set it to the initial value function at that
// artificial time.
- // We create an output of the initial value. Then we also need to collect
- // the two starting solutions in a <tt>std::vector</tt> of pointers field
- // and to set up an instance of the <code> SineGordonOperation class </code>
- // based on the finite element degree specified at the top of this file.
+ // We then go on by writing the initial state to file and collecting
+ // the two starting solutions in a <tt>std::vector</tt> of pointers that
+ // get later consumed by the SineGordonOperation::apply() function. Next,
+ // an instance of the <code> SineGordonOperation class </code> based on
+ // the finite element degree specified at the top of this file is set up.
VectorTools::interpolate(dof_handler,
- ExactSolution<dim>(1, time),
+ InitialCondition<dim>(1, time),
solution);
VectorTools::interpolate(dof_handler,
- ExactSolution<dim>(1, time - time_step),
+ InitialCondition<dim>(1, time - time_step),
old_solution);
output_results(0);
std::vector<LinearAlgebra::distributed::Vector<double> *>
- previous_solutions;
- previous_solutions.push_back(&old_solution);
- previous_solutions.push_back(&old_old_solution);
+ previous_solutions({&old_solution, &old_old_solution});
SineGordonOperation<dim, fe_degree> sine_gordon_op(matrix_free_data,
time_step);
// Now loop over the time steps. In each iteration, we shift the solution
- // vectors by one and call the <code> apply </code> function of the <code>
- // SineGordonOperator </code>. Then, we write the solution to a file. We
+ // vectors by one and call the `apply` function of the
+ // `SineGordonOperator` class. Then, we write the solution to a file. We
// clock the wall times for the computational time needed as wall as the
// time needed to create the output and report the numbers when the time
// stepping is finished.
//
// Note how this shift is implemented: We simply call the swap method on
// the two vectors which swaps only some pointers without the need to copy
- // data around. Obviously, this is a more efficient way to update the
- // vectors during time stepping. Let us see what happens in more detail:
- // First, we exchange <code>old_solution</code> with
- // <code>old_old_solution</code>, which means that
- // <code>old_old_solution</code> gets <code>old_solution</code>, which is
- // what we expect. Similarly, <code>old_solution</code> gets the content
- // from <code>solution</code> in the next step. Afterward,
- // <code>solution</code> holds <code>old_old_solution</code>, but that
- // will be overwritten during this step.
+ // data around, a relatively expensive operation within an explicit time
+ // stepping method. Let us see what happens in more detail: First, we
+ // exchange <code>old_solution</code> with <code>old_old_solution</code>,
+ // which means that <code>old_old_solution</code> gets
+ // <code>old_solution</code>, which is what we expect. Similarly,
+ // <code>old_solution</code> gets the content from <code>solution</code>
+ // in the next step. After this, <code>solution</code> holds
+ // <code>old_old_solution</code>, but that will be overwritten during this
+ // step.
unsigned int timestep_number = 1;
Timer timer;
projects / dealii.git / commitdiff