From 5c2cd2f0359654899fb2e6696629e23cbc97fd85 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 3 Jan 2020 10:48:11 -0600 Subject: [PATCH] update documentation, part IV --- examples/step-69/step-69.cc | 284 +++++++++++++++++++----------------- 1 file changed, 150 insertions(+), 134 deletions(-) diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 95367d3c8d..57efdf5dd8 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -64,7 +64,7 @@ // In addition to above deal.II specific includes, we also include four // boost headers. The first two are for binary archives that we will use -// for implementing a checkpointing and restart mechanism. +// for implementing a check-pointing and restart mechanism. #include #include @@ -74,9 +74,24 @@ #include // @sect3{Class template declarations} +// +// We begin our actual implementation by declaring all classes with its +// data structures and methods upfront. In contrast to previous example +// steps we use a more fine-grained encapsulation of concepts, data +// structures, and parameters into individual classes. A single class thus +// usually centers around either a single data structure (such as the +// Triangulation) in the Discretization class, or a single +// method (such as the step() function of the +// TimeStep class). We typically declare parameter variables +// and scratch data object private and make methods and data structures +// used by other classes public. +// +// @note: A cleaner approach would be to guard access to all data +// structures by getter/setter +// functions. For the sake of brevity, we refrain from that approach, +// though. -// Next we declare all data structures and class templates of the example -// step. namespace Step69 { using namespace dealii; @@ -88,18 +103,18 @@ namespace Step69 enum Boundary : types::boundary_id { do_nothing = 0, - slip = 2, - dirichlet = 3, + slip = 1, + dirichlet = 2, }; - // @sect4{class Discretization} + // @sect4{The Discretization class} // // The class Discretization contains all data structures // concerning the mesh (triangulation) and discretization (mapping, // finite element, quadrature) of the problem. We use the // ParameterAcceptor class to automatically populate problem-specific // parameters, such as the geometry information - // (immersed_disc_length, etc.) or the refinement level + // (length, etc.) or the refinement level // (refinement) from a parameter file. This requires us to // split the initialization of data structures into two functions: We // initialize everything that does not depend on parameters in the @@ -129,15 +144,15 @@ namespace Step69 private: TimerOutput &computing_timer; - double immersed_disc_length; - double immersed_disc_height; - double immersed_disc_object_position; - double immersed_disc_object_diameter; + double length; + double height; + double disc_position; + double disc_diameter; unsigned int refinement; }; - // @sect4{class OfflineData} + // @sect4{The OfflineData class} // // The class OfflineData contains pretty much all components // of the discretization that do not evolve in time, in particular, the @@ -208,7 +223,7 @@ namespace Step69 SmartPointer> discretization; }; - // @sect4{class ProblemDescription} + // @sect4{The ProblemDescription class} // // The member functions of this class are utility functions specific to // Euler's equations: @@ -274,7 +289,7 @@ namespace Step69 const Tensor<1, dim> &n_ij); }; - // @sect4{class InitialValues} + // @sect4{The InitialValues class} // // The class InitialValues's only public data attribute is a // std::function initial_state that computes the initial @@ -305,7 +320,7 @@ namespace Step69 Tensor<1, 3> initial_1d_state; }; - // @sect4{class TimeStep} + // @sect4{The TimeStep class} // // With the OfflineData and ProblemDescription // classes at hand we can now implement the explicit time-stepping scheme @@ -354,7 +369,7 @@ namespace Step69 double cfl_update; }; - // @sect4{class SchlierenPostprocessor} + // @sect4{The SchlierenPostprocessor class} // // At its core, the Schlieren class implements the class member // compute_schlieren. The main purpose of this class member @@ -404,7 +419,7 @@ namespace Step69 double schlieren_beta; }; - // @sect4{class TimeLoop} + // @sect4{The TimeLoop class} // // Now, all that is left to do is to chain the methods implemented in the // TimeStep, InitialValues, and @@ -457,12 +472,18 @@ namespace Step69 vector_type output_vector; }; - // @sect3{Class template implementations} + // @sect3{Grid generation and assembly} - // @sect4{Implementation of the members of the class Discretization} - - // Not much is done here other that initializing the corresponding - // class members in the initialization list. + // The first major task at hand is the typical triplet of grid + // generation, setup of data structures, and assembly. A notable novelty + // in this example step is the use of the ParameterAcceptor class that we + // use to populate parameter values: We first initialize the + // ParameterAcceptor class by calling its constructor with a string + // subsection denoting the correct subsection in the + // parameter file. Then, in the constructor body every parameter value is + // initialized to a sensible default value and registered with the + // ParameterAcceptor class with a call to + // ParameterAcceptor::add_parameter. template Discretization::Discretization(const MPI_Comm & mpi_communicator, @@ -477,24 +498,24 @@ namespace Step69 , face_quadrature(3) , computing_timer(computing_timer) { - immersed_disc_length = 4.; + length = 4.; add_parameter("immersed disc - length", - immersed_disc_length, + length, "Immersed disc: length of computational domain"); - immersed_disc_height = 2.; + height = 2.; add_parameter("immersed disc - height", - immersed_disc_height, + height, "Immersed disc: height of computational domain"); - immersed_disc_object_position = 0.6; + disc_position = 0.6; add_parameter("immersed disc - object position", - immersed_disc_object_position, + disc_position, "Immersed disc: x position of immersed disc center point"); - immersed_disc_object_diameter = 0.5; + disc_diameter = 0.5; add_parameter("immersed disc - object diameter", - immersed_disc_object_diameter, + disc_diameter, "Immersed disc: diameter of immersed disc"); refinement = 5; @@ -505,38 +526,35 @@ namespace Step69 // Note that in the previous constructor we only passed the MPI // communicator to the triangulationbut we still have not - // initialized the underlying geometry/mesh. In order to define the geometry - // we will use the class create_immersed_disc_geometry - // that uses the tools in GridGenerator in order to create a - // rectangular domain with a whole. - - // The following is just a dummy implementation/placeholder that does - // nothing other than throwing an exception if we want to run this program - // with a space dimension that is not 2. + // initialized the underlying geometry/mesh. As mentioned earlier, we + // have to postpone this task to the setup() function that + // gets called after the ParameterAcceptor::initialize() function has + // populated all parameter variables with the final values read from the + // parameter file. + // + // The setup() function is the last class member that has to + // be implemented. It creates the actual triangulation that is a + // benchmark configuration consisting of a channel with a disc obstacle, see + // @cite GuermondEtAl2018. We construct the geometry by modifying the + // mesh generated by GridGenerator::hyper_cube_with_cylindrical_hole(). + // We refer to Step-49, Step-53, and Step-54 for an overview how to + // create advanced meshes. template - void - create_immersed_disc_geometry(parallel::distributed::Triangulation &, - const double /*length*/, - const double /*height*/, - const double /*step_position*/, - const double /*step_height*/) + void Discretization::setup() { - AssertThrow(false, ExcNotImplemented()); - } + TimerOutput::Scope t(computing_timer, "discretization - setup"); - // For the two-dimensional case we have the following template - // specialization that creates the geometry. + triangulation.clear(); - template <> - void create_immersed_disc_geometry<2>( - parallel::distributed::Triangulation<2> &triangulation, - const double length, - const double height, - const double disc_position, - const double disc_diameter) - { - constexpr int dim = 2; + // We first create 4 temporary (non distributed) coarse triangulations + // that we stitch together with the + // GridGenerator::merge_triangulation() function. We center the disc at + // $(0,0)$ with a diameter of disc_diameter. The lower + // left corner of the channel has coordinates + // (-disc_position, -height/2) and the upper + // right corner has (length-disc_position, + // height/2). Triangulation tria1, tria2, tria3, tria4; @@ -568,6 +586,10 @@ namespace Step69 triangulation.set_manifold(0, PolarManifold<2>(Point<2>())); + // We have to fix up the left edge that is currently located at + // $x=-$disc_diameter and has to be shifted to + // $x=-$disc_position: + for (auto cell : triangulation.active_cell_iterators()) for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) { @@ -576,6 +598,12 @@ namespace Step69 vertex[0] = -disc_position; } + // As a last step the boundary has to be colorized with + // Boundary::do_nothing on the right, + // Boundary::dirichlet on the left and + // Boundary::slip on the upper and lower outer boundaries + // and the obstacle: + for (auto cell : triangulation.active_cell_iterators()) { for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) @@ -595,34 +623,13 @@ namespace Step69 face->set_boundary_id(Boundary::slip); } } - } - - // This the the last class member to be implemented of the class. - // Discretization: it initializes the actual mesh - // of the triangulation by calling create_immersed_disc_geometry. - - template - void Discretization::setup() - { - TimerOutput::Scope t(computing_timer, "discretization - setup"); - - triangulation.clear(); - - create_immersed_disc_geometry(triangulation, - immersed_disc_length, - immersed_disc_height, - immersed_disc_object_position, - immersed_disc_object_diameter); triangulation.refine_global(refinement); } - - // @sect4{Implementation of the members of the class OfflineData} - - // Not much is done here other that initializing the corresponding - // class members in the initialization list. - // Constructor of the class OfflineData. + // Not much is done in the constructor of OfflineData other + // than initializing the corresponding class members in the + // initialization list. template OfflineData::OfflineData(const MPI_Comm & mpi_communicator, @@ -635,11 +642,10 @@ namespace Step69 , discretization(&discretization) {} - // Now the class member OfflineData::setup() will take care - // of initializating - // - The dof_handler. - // - The IndexSets corresponding to locally owned and locally relevant DOFs. - // - The partitioner. + // Now we can initialize the DoFHandler, extract the IndexSet objects for + // locally owned and locally relevant DOFs, and initialize a + // Utilities::MPI::Partitioner object that is needed for distributed + // vectors. template void OfflineData::setup() @@ -652,18 +658,12 @@ namespace Step69 dof_handler.initialize(discretization->triangulation, discretization->finite_element); + DoFRenumbering::Cuthill_McKee(dof_handler); locally_owned = dof_handler.locally_owned_dofs(); n_locally_owned = locally_owned.n_elements(); - } - - { - TimerOutput::Scope t( - computing_timer, - "offline_data - create partitioner and affine constraints"); - locally_relevant.clear(); DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant); n_locally_relevant = locally_relevant.n_elements(); @@ -674,41 +674,60 @@ namespace Step69 const auto dofs_per_cell = discretization->finite_element.dofs_per_cell; - // Here we create the sparsity patterns for the off-line data. There are - // quite a few peculiarities that deserve our attention: - // - Our "local" dynamic sparsity pattern (dsp) - // will be of dimensions - // n_locally_relevant $\times$ - // n_locally_relevant (this choice is definitely unusual). - // The goal behind such choice is to reduce communication: we do - // not want to request (to another mpi process) ghosted offline data (such - // as the vectors $\mathbf{c}_{ij}$ when $j$ is not locally owned) for - // every time step. It is more efficient to simply take more memory by - // storing (locally) all relevant off-line data. - // - We loop on all locally owned and ghosted cells (see @ref - // GlossArtificialCell "this glossary entry") order to: - //
    - //
  • Extract the dof_indices associated to the cell DOFs - // (having global numbering) and renumber them using - // partitioner->global_to_local(index). For the case - // of locally owned DOFs: such renumbering consist in applying a - // shift (i.e. we subtract a number) such that now they will - // become a number in the integer interval - // $[0,$n_locally_owned$)$. However, for the case of - // "ghosted DOFs" (i.e. not locally owned) the situation is quite - // different, since the global indices associated to ghosted DOFs - // will not be (in general) a contiguous set of integers. - //
    • - //
  • Once, we are done with that, we add the corresponding entries to - // the rows of the dynamic sparsity pattern with - // dsp.add_entries
    • - //
- // Finally we use dsp to initialize the actual sparsity - // pattern sparsity_pattern. + // @sect4{Translation to local index ranges} + + // We are now in a position to create the sparsity pattern for our + // matrices. There are quite a few peculiarities that need a detailed + // explanation. We avoid using a distributed matrix class (as for + // example provided by Trilinos or PETSc) and instead rely on deal.II's + // own SparseMatrix object to store the local part of all matrices. + // This design decision is motivated by the fact that we actually never + // perform a matrix-vector multiplication. Instead, we will have to + // perform nonlinear updates while iterating over (the local part) of a + // connectivity stencil; a task for which deal.II own SparsityPattern + // is specificially optimized for. + // + // This design consideration has a caveat, though. What makes the + // deal.II SparseMatrix class fast is the compressed row + // storage (CSR) used in the SparsityPattern (see @ref Sparsity). + // This, unfortunately, does not play nicely with a global distributed + // index range because a sparsity pattern with CSR cannot contain + // "holes" in the index range. The distributed matrices offered by + // deal.II avoid this by translating from a global index range into a + // contiguous local index range. But this is the precisely the type of + // index manipulation we want to avoid. + // + // Lucky enough, the Utilities::MPI::Partitioner used for distributed + // vectors provides exactly what we need: It manages a translation from + // a global index range to a contiguous local (per MPI rank) index + // range. We therefore simply create a "local" sparsity pattern for the + // contiguous index range $[0,$n_locally_relevant$)$ and + // translate between global dof indices and the above local range with + // the help of the Utilities::MPI::Partitioner::global_to_local() + // function. All that is left to do is to ensure that we always access + // elements of a distributed vector by a call to + // LinearAlgebra::distributed::Vector::local_element(). That way we + // avoid index translations altogether. { - TimerOutput::Scope t(computing_timer, - "offline_data - create sparsity pattern"); + TimerOutput::Scope t( + computing_timer, + "offline_data - create sparsity pattern and set up matrices"); + + // We have to create the "local" sparsity pattern by hand. We + // therefore loop over all locally owned and ghosted cells (see @ref + // GlossArtificialCell) and extract the (global) + // dof_indices associated to the cell DOFs and renumber + // them using partitioner->global_to_local(index). + // + // @note In the case of a locally owned dof, such renumbering consist + // of applying a shift (i.e. we subtract an offset) such that now they + // will become a number in the integer interval + // $[0,$n_locally_owned$)$. However, in the case of a + // ghosted dof (i.e. not locally owned) the situation is quite + // different, since the global indices associated to ghosted DOFs will + // not be (in general) a contiguous set of integers. DynamicSparsityPattern dsp(n_locally_relevant, n_locally_relevant); @@ -732,13 +751,6 @@ namespace Step69 } sparsity_pattern.copy_from(dsp); - } - - // We initialize the off-line data matrices. Note that these matrices do - // not require an mpi communicator (that's the idea). - - { - TimerOutput::Scope t(computing_timer, "offline_data - set up matrices"); lumped_mass_matrix.reinit(sparsity_pattern); norm_matrix.reinit(sparsity_pattern); @@ -1320,6 +1332,8 @@ namespace Step69 } } /* assemble() */ + // @sect3{Problem specific setup and approximate Riemann solver} + // At this point we are very much done with anything related to offline data. // // Now we define the implementation of the utility @@ -1588,6 +1602,8 @@ namespace Step69 "m_3", "E"}; + // @sect3{Initial values and time stepping} + // Implementation of the constructor for the class InitialValues. template @@ -2115,7 +2131,7 @@ namespace Step69 schlieren.update_ghost_values(); } - // Placeholder here. + // @sect3{The Timeloop::run() function} template TimeLoop::TimeLoop(const MPI_Comm &mpi_comm) -- 2.39.5