From: Wolfgang Bangerth Date: Fri, 21 Jun 2024 23:20:11 +0000 (-0600) Subject: Edit step-89.cc. X-Git-Tag: v9.6.0-rc1~152^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5da64faaa8476a862679a37e5c6446542c7892ce;p=dealii.git Edit step-89.cc. --- diff --git a/examples/step-89/step-89.cc b/examples/step-89/step-89.cc index 8621149067..9a1322ab01 100644 --- a/examples/step-89/step-89.cc +++ b/examples/step-89/step-89.cc @@ -58,19 +58,21 @@ namespace Step89 // @sect3{Initial conditions for vibrating membrane} // - // Function that provides the initial condition for the vibrating membrane - // test case. + // In the following, let us first define a function that provides + // the initial condition for the vibrating membrane test case. It + // implementes both the initial pressure (component 0) and velocity + // (components 1 to 1 + dim). + // + // There is also a function that computes the duration of one + // oscillation. template class InitialConditionVibratingMembrane : public Function { public: InitialConditionVibratingMembrane(const double modes); - // Function that the gives the initial pressure (comp 0) and velocity (comp - // 1 to 1 + dim). double value(const Point &p, const unsigned int comp) const final; - // Function that calculates the duration of one oscillation. double get_period_duration(const double speed_of_sound) const; private: @@ -107,15 +109,16 @@ namespace Step89 // @sect3{Gauss pulse} // - // Function that provides the values of a pressure Gauss pulse. + // Next up is a function that provides the values of a pressure + // Gauss pulse. As with the previous function, it implements both + // the initial pressure (component 0) and velocity (components 1 to + // 1 + dim). template class GaussPulse : public Function { public: GaussPulse(const double shift_x, const double shift_y); - // Function that the gives the initial pressure (comp 0) and velocity (comp - // 1 to 1 + dim). double value(const Point &p, const unsigned int comp) const final; private: @@ -132,8 +135,6 @@ namespace Step89 static_assert(dim == 2, "Only implemented for dim==2"); } - // Function that the gives the initial pressure (comp 0) and velocity (comp 1 - // to 1 + dim). template double GaussPulse::value(const Point &p, const unsigned int comp) const @@ -153,9 +154,7 @@ namespace Step89 { // Helper function to check if a boundary ID is related to a non-matching // face. A @c std::set that contains all non-matching boundary IDs is - // handed over additionally to the face ID under question. This function - // could certainly also be defined inline but this way the code is more easy - // to read. + // handed over in addition to the face ID under question. bool is_non_matching_face( const std::set &non_matching_face_ids, const types::boundary_id face_id) @@ -213,12 +212,27 @@ namespace Step89 } } // namespace HelperFunctions - //@sect3{Material access} + // @sect3{Material parameter description} + // + // The following class stores the information if the fluid is + // homogeneous as well as the material properties at every cell. + // This class helps access the correct values without accessing a + // large vector of materials in the homogeneous case. // - // This class stores the information if the fluid is homogeneous - // as well as the material properties at every cell. - // This class helps to access the correct values without accessing - // a large vector of materials in the homogeneous case. + // The class is provided a map from material ids to material + // properties (given as a pair of values for the speed of sound and + // the density). If the map has only one entry, the material is + // homogeneous -- using the same values everywhere -- and we can + // remember that fact in the `homogeneous` member variable and use it + // to optimize some code paths below. If the material is not + // homogeneous, we will fill a vector with the correct materials for + // each batch of cells; this information can then be access via + // FEEvaluationData::read_cell_data(). + // + // As is usual when working with the MatrixFree framework, we will + // not only access material parameters at a single site, but for + // whole "batches" of cells. As a consequence, the functions below + // returned VectorizedArray objects for a batch at a time. template class CellwiseMaterialData { @@ -228,24 +242,18 @@ namespace Step89 const MatrixFree> &matrix_free, const std::map> &material_id_map) - // If the map is of size 1, the material is constant in every cell. : homogeneous(material_id_map.size() == 1) { Assert(material_id_map.size() > 0, - ExcMessage("No materials given to CellwiseMaterialData")); + ExcMessage("No materials given to CellwiseMaterialData.")); if (homogeneous) { - // In the homogeneous case we know the materials in the whole domain. speed_of_sound_homogeneous = material_id_map.begin()->second.first; density_homogeneous = material_id_map.begin()->second.second; } else { - // In the in-homogeneous case materials vary between cells. We are - // filling a vector with the correct materials, that can be processed - // via - // @c read_cell_data(). const auto n_cell_batches = matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); @@ -278,47 +286,49 @@ namespace Step89 const AlignedVector> &get_speed_of_sound() const { - Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()")); + Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound().")); return speed_of_sound; } const AlignedVector> &get_density() const { - Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()")); + Assert(!homogeneous, ExcMessage("Use get_homogeneous_density().")); return density; } VectorizedArray get_homogeneous_speed_of_sound() const { - Assert(homogeneous, ExcMessage("Use get_speed_of_sound()")); + Assert(homogeneous, ExcMessage("Use get_speed_of_sound().")); return speed_of_sound_homogeneous; } VectorizedArray get_homogeneous_density() const { - Assert(homogeneous, ExcMessage("Use get_density()")); + Assert(homogeneous, ExcMessage("Use get_density().")); return density_homogeneous; } private: const bool homogeneous; - // Materials in the in-homogeneous case. + /* Materials in the inhomogeneous case. */ AlignedVector> speed_of_sound; AlignedVector> density; - // Materials in the homogeneous case. + /* Materials in the homogeneous case. */ VectorizedArray speed_of_sound_homogeneous; VectorizedArray density_homogeneous; }; - // To be able to access the material data in every cell in a thread safe way - // @c MaterialEvaluation is used. Similar to @c FEEvaluation, every thread - // creates its own instance and thus, there are no race conditions. For - // in-homogeneous materials, a @c reinit_cell() or @c reinit_face() function - // is used to set the correct material at the current cell batch. In the - // homogeneous case the @c _reinit() functions don't have to reset the - // materials. + // To be able to access the material data in every cell in a thread + // safe way, the following class @c MaterialEvaluation is + // used. Similar to @c FEEvaluation, functions below will create + // their own instances of this class; thus, there can be no race + // conditions even if these functions run multiple times in + // parallel. For inhomogeneous materials, a @c reinit_cell() or @c + // reinit_face() function is used to set the correct material at the + // current cell batch. In the homogeneous case the @c _reinit() + // functions don't have to reset the materials. template class MaterialEvaluation { @@ -332,7 +342,6 @@ namespace Step89 { if (material_data.is_homogeneous()) { - // Set the material that is used in every cell. speed_of_sound = material_data.get_homogeneous_speed_of_sound(); density = material_data.get_homogeneous_density(); } @@ -343,13 +352,14 @@ namespace Step89 return material_data.is_homogeneous(); } - // Update the cell data, given a cell batch index. + // The following two functions update the data for the current + // cell or face, given a cell batch index. If the material is + // homogeneous, there is nothing to do. Otherwise, we reinit the + // FEEvaluation object and store the data for the current object. void reinit_cell(const unsigned int cell) { - // In the homogeneous case we do not have to reset the cell data. if (!material_data.is_homogeneous()) { - // Reinit the FEEvaluation object and set the cell data. phi.reinit(cell); speed_of_sound = phi.read_cell_data(material_data.get_speed_of_sound()); @@ -357,13 +367,10 @@ namespace Step89 } } - // Update the cell data, given a face batch index. void reinit_face(const unsigned int face) { - // In the homogeneous case we do not have to reset the cell data. if (!material_data.is_homogeneous()) { - // Reinit the FEFaceEvaluation object and set the cell data. phi_face.reinit(face); speed_of_sound = phi_face.read_cell_data(material_data.get_speed_of_sound()); @@ -371,33 +378,33 @@ namespace Step89 } } - // Return the speed of sound at the current cell batch. + // The following functions then return the speed of sound and + // density for the current cell batch. VectorizedArray get_speed_of_sound() const { return speed_of_sound; } - // Return the density at the current cell batch. VectorizedArray get_density() const { return density; } private: - // Members needed for the in-homogeneous case. + /* Members needed for the inhomogeneous case. */ FEEvaluation phi; FEFaceEvaluation phi_face; - // Material defined at every cell. + /* Material defined at every cell. */ const CellwiseMaterialData &material_data; - // Materials at current cell. + /* Materials at current cell. */ VectorizedArray speed_of_sound; VectorizedArray density; }; - //@sect3{Boundary conditions} + // @sect3{Boundary conditions} // // To be able to use the same kernel, for all face integrals we define // a class that returns the needed values at boundaries. In this tutorial @@ -441,17 +448,23 @@ namespace Step89 const FEFaceEvaluation &velocity_m; }; - //@sect3{Acoustic operator} + // @sect3{Acoustic operator} + // + // The following class then defines the acoustic operator. The class is + // heavily based on matrix-free methods. For a better understanding in + // matrix-free methods please refer to step-67. // - // Class that defines the acoustic operator. The class is heavily based on - // matrix-free methods. For a better understanding in matrix-free methods - // please refer to step-67. + // At the top we define a flag that decides whether we want to use + // mortaring. If the remote evaluators are set up with a + // VectorizedArray we are using point-to-point interpolation; + // otherwise we make use of Nitsche-type mortaring. The decision is + // made using `std::is_floating_point_v`, which is a variable that + // is true if the template argument is a floating point type (such + // as `float` or `double`) and false otherwise (in particular if the + // template argument is a vectorized array type). template class AcousticOperator { - // If the remote evaluators are set up with a VectorizedArray we are - // using point-to-point interpolation. Otherwise we make use of - // Nitsche-type mortaring. static constexpr bool use_mortaring = std::is_floating_point_v; @@ -486,20 +499,26 @@ namespace Step89 has to be provided.")); } - // Function to evaluate the acoustic operator. + // The following function then evaluates the acoustic operator. + // It first updates the precomputed values in corresponding the + // FERemoteEvaluation objects. The material parameters do not change and + // thus, we do not have to update precomputed values in @c c_r_eval and @c + // rho_r_eval. + // + // It then either performs a matrix-free loop with Nitsche-type + // mortaring at non-matching faces (if `use_mortaring` is true) or + // with point-to-point interpolation at non-matching faces (in the + // `else` branch). The difference is only in the third argument to + // the loop function, denoting the function object that is + // executed at boundary faces. template void evaluate(VectorType &dst, const VectorType &src) const { - // Update the precomputed values in corresponding the FERemoteEvaluation - // objects. The material parameters do not change and thus, we do - // not have to update precomputed values in @c c_r_eval and @c rho_r_eval. pressure_r_eval->gather_evaluate(src, EvaluationFlags::values); velocity_r_eval->gather_evaluate(src, EvaluationFlags::values); if constexpr (use_mortaring) { - // Perform matrix free loop with Nitsche-type mortaring at - // non-matching faces. matrix_free.loop( &AcousticOperator::local_apply_cell, &AcousticOperator::local_apply_face, @@ -513,8 +532,6 @@ namespace Step89 } else { - // Perform matrix free loop with point-to-point interpolation at - // non-matching faces. matrix_free.loop( &AcousticOperator::local_apply_cell, &AcousticOperator::local_apply_face, @@ -529,8 +546,14 @@ namespace Step89 } } + // In the `private` section of the class, we then define the + // functions that evaluate volume, interior face, and boundary + // face integrals. The concrete terms these functions evaluate are + // stated in the documentation at the top of this tutorial + // program. Each of these functions has its own object of type + // `MaterialEvaluation` that provides access to the material at + // each cell or face. private: - // This function evaluates the volume integrals. template void local_apply_cell( const MatrixFree &matrix_free, @@ -541,7 +564,6 @@ namespace Step89 FEEvaluation pressure(matrix_free, 0, 0, 0); FEEvaluation velocity(matrix_free, 0, 0, 1); - // Class that gives access to the material at each cell MaterialEvaluation material(matrix_free, *material_data); for (unsigned int cell = cell_range.first; cell < cell_range.second; @@ -553,13 +575,15 @@ namespace Step89 pressure.gather_evaluate(src, EvaluationFlags::gradients); velocity.gather_evaluate(src, EvaluationFlags::gradients); - // Get the materials at the corresponding cell. Since we + // Get the materials on the corresponding cell. Since we // introduced @c MaterialEvaluation we can write the code - // independent if the material is homogeneous or in-homogeneous. + // independent of whether the material is homogeneous or + // inhomogeneous. material.reinit_cell(cell); const auto c = material.get_speed_of_sound(); const auto rho = material.get_density(); - for (unsigned int q : pressure.quadrature_point_indices()) + + for (const unsigned int q : pressure.quadrature_point_indices()) { pressure.submit_value(rho * c * c * velocity.get_divergence(q), q); @@ -571,13 +595,18 @@ namespace Step89 } } - // This function evaluates the fluxes at faces between cells with the same - // material. If boundary faces are under consideration fluxes into + // This next function evaluates the fluxes at faces between cells with the + // same material. If boundary faces are under consideration fluxes into // neighboring faces do not have to be considered which is enforced via // `weight_neighbor=false`. For non-matching faces the fluxes into // neighboring faces are not considered as well. This is because we iterate // over each side of the non-matching face separately (similar to a cell // centric loop). + // + // In this and following functions, we also introduce the factors + // $\tau$ and $\gamma$ that are computed from $\rho$ and $c$ along + // interfaces and that appear in the bilinear forms. Their + // definitions are provided in the introduction. template velocity_p( matrix_free, false, 0, 0, 1); - // Class that gives access to the material at each cell MaterialEvaluation material(matrix_free, *material_data); for (unsigned int face = face_range.first; face < face_range.second; - face++) + ++face) { velocity_m.reinit(face); velocity_p.reinit(face); @@ -728,9 +752,9 @@ namespace Step89 } - //@sect4{Matrix-free boundary function for point-to-point interpolation} + // @sect4{Matrix-free boundary function for point-to-point interpolation} // - // This function evaluates the boundary face integrals and the + // The following function then evaluates the boundary face integrals and the // non-matching face integrals using point-to-point interpolation. template void local_apply_boundary_face_point_to_point( @@ -739,17 +763,14 @@ namespace Step89 const VectorType &src, const std::pair &face_range) const { - // Standard face evaluators. FEFaceEvaluation pressure_m( matrix_free, true, 0, 0, 0); FEFaceEvaluation velocity_m( matrix_free, true, 0, 0, 1); - // Classes that return the correct BC values. BCEvaluationP pressure_bc(pressure_m); BCEvaluationU velocity_bc(velocity_m); - // Class that gives access to the material at each cell MaterialEvaluation material(matrix_free, *material_data); // Remote evaluators. @@ -759,7 +780,7 @@ namespace Step89 auto rho_r = rho_r_eval->get_data_accessor(); for (unsigned int face = face_range.first; face < face_range.second; - face++) + ++face) { velocity_m.reinit(face); pressure_m.reinit(face); @@ -787,12 +808,12 @@ namespace Step89 material.reinit_face(face); + // If we are considering a homogeneous material, do not use the + // inhomogeneous fluxes. While it would be possible + // to use the inhomogeneous fluxes they are more expensive to + // compute. if (material.is_homogeneous()) { - // If homogeneous material is considered do not use the - // inhomogeneous fluxes. While it would be possible - // to use the inhomogeneous fluxes they are more expensive to - // compute. evaluate_face_kernel(pressure_m, velocity_m, pressure_r, @@ -802,8 +823,6 @@ namespace Step89 } else { - // If inhomogeneous material is considered use the - // in-homogeneous fluxes. c_r.reinit(face); rho_r.reinit(face); evaluate_face_kernel_inhomogeneous( @@ -838,7 +857,7 @@ namespace Step89 } } - //@sect4{Matrix-free boundary function for Nitsche-type mortaring} + // @sect4{Matrix-free boundary function for Nitsche-type mortaring} // // This function evaluates the boundary face integrals and the // non-matching face integrals using Nitsche-type mortaring. @@ -849,7 +868,6 @@ namespace Step89 const VectorType &src, const std::pair &face_range) const { - // Standard face evaluators for BCs. FEFaceEvaluation pressure_m( matrix_free, true, 0, 0, 0); FEFaceEvaluation velocity_m( @@ -929,13 +947,12 @@ namespace Step89 velocity_r_mortar.reinit(face * n_lanes + v); pressure_r_mortar.reinit(face * n_lanes + v); + // As above, if we are considering a homogeneous + // material, do not use the inhomogeneous + // fluxes. Since we are operating on face @c v we + // call @c material.get_density()[v]. if (material.is_homogeneous()) { - // If homogeneous material is considered do not use the - // inhomogeneous fluxes. While it would be possible - // to use the inhomogeneous fluxes they are more - // expensive to compute. Since we are operating on face @c - // v we call @c material.get_density()[v]. evaluate_face_kernel( pressure_m_mortar, velocity_m_mortar, @@ -983,7 +1000,6 @@ namespace Step89 } else { - // Same as in @c local_apply_boundary_face_point_to_point(). velocity_m.reinit(face); pressure_m.reinit(face); @@ -1006,8 +1022,6 @@ namespace Step89 const MatrixFree &matrix_free; - // CellwiseMaterialData is stored as shared pointer with the same - // argumentation. const std::shared_ptr> material_data; const std::set remote_face_ids; @@ -1064,9 +1078,12 @@ namespace Step89 } }; - //@sect3{Inverse mass operator} + // @sect3{Inverse mass operator} // - // Class to apply the inverse mass operator. + // For the time stepping methods below, we need the inverse mass + // operator. We apply it via a loop over all (batches of) cells as + // always, where the contribution of each cell is computed in a + // matrix-free way: template class InverseMassOperator { @@ -1075,7 +1092,6 @@ namespace Step89 : matrix_free(matrix_free) {} - // Function to apply the inverse mass operator. template void apply(VectorType &dst, const VectorType &src) const { @@ -1087,7 +1103,6 @@ namespace Step89 } private: - // Apply the inverse mass operator onto every cell batch. template void local_apply_cell( const MatrixFree &mf, @@ -1112,7 +1127,7 @@ namespace Step89 const MatrixFree &matrix_free; }; - //@sect3{Runge-Kutta time-stepping} + // @sect3{Runge-Kutta time-stepping} // // This class implements a Runge-Kutta scheme of order 2. template @@ -1130,7 +1145,14 @@ namespace Step89 , acoustic_operator(acoustic_operator) {} - // Set up and run time loop. + // The `run()` function of this class contains the time loop. It + // starts by initializing some member variables (such as + // short-cuts to objects that describe the finite element, its + // properties, and the mapping; an by initializing vectors). It + // then computes the time step size via minimum element edge + // length. Assuming non-distorted elements, we can compute the + // edge length as the distance between two vertices. From this, + // we can then obtain the time step size via the CFL condition. void run(const MatrixFree &matrix_free, const double cr, const double end_time, @@ -1138,25 +1160,19 @@ namespace Step89 const Function &initial_condition, const std::string &vtk_prefix) { - // Get needed members of matrix free. const auto &dof_handler = matrix_free.get_dof_handler(); const auto &mapping = *matrix_free.get_mapping_info().mapping; const auto degree = dof_handler.get_fe().degree; - // Initialize needed Vectors. VectorType solution; matrix_free.initialize_dof_vector(solution); VectorType solution_temp; matrix_free.initialize_dof_vector(solution_temp); - // Set the initial condition. HelperFunctions::set_initial_condition(matrix_free, initial_condition, solution); - // Compute time step size: Compute minimum element edge length. - // We assume non-distorted elements, therefore we only compute - // the distance between two vertices double h_local_min = std::numeric_limits::max(); for (const auto &cell : dof_handler.active_cell_iterators()) h_local_min = @@ -1166,16 +1182,14 @@ namespace Step89 const double h_min = Utilities::MPI::min(h_local_min, dof_handler.get_communicator()); - // Compute constant time step size via the CFL condition. const double dt = cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound); - // Perform time integration loop. + // We can then perform the time integration loop: double time = 0.0; unsigned int timestep = 0; while (time < end_time) { - // Write output. HelperFunctions::write_vtu(solution, matrix_free.get_dof_handler(), mapping, @@ -1183,40 +1197,48 @@ namespace Step89 "step_89-" + vtk_prefix + std::to_string(timestep)); - // Perform a single time step. std::swap(solution, solution_temp); time += dt; - timestep++; + ++timestep; perform_time_step(dt, solution, solution_temp); } } + // The main work of this class is done by a `private` member + // function that performs one Runge-Kutta 2 time step. Recall that + // this method requires two sub-steps ("stages") computing + // intermediate values `k1` and `k2`, after which the intermediate + // values are summed with weights to obtain the new solution at + // the end of the time step. The RK2 method allows for the + // elimination of the intermediate vector `k2` by utilizing the + // `dst` vector as temporary storage. private: - // Perform one Runge-Kutta 2 time step. void perform_time_step(const double dt, VectorType &dst, const VectorType &src) { VectorType k1 = src; - // First stage. + /* First stage. */ evaluate_stage(k1, src); - // Second stage. + /* Second stage. */ k1.sadd(0.5 * dt, 1.0, src); evaluate_stage(dst, k1); + + /* Summing things into the output vector. */ dst.sadd(dt, 1.0, src); } - // Evaluate a single Runge-Kutta stage. + // Evaluating a single Runge-Kutta stage is a straightforward step + // that really only requires applying the operator once, and then + // applying the inverse of the mass matrix. void evaluate_stage(VectorType &dst, const VectorType &src) { - // Evaluate the stage acoustic_operator->evaluate(dst, src); dst *= -1.0; inverse_mass_operator->apply(dst, dst); } - // Needed operators. const std::shared_ptr> inverse_mass_operator; const std::shared_ptr> @@ -1226,11 +1248,40 @@ namespace Step89 // @sect3{Construction of non-matching triangulations} // - // This function creates a two dimensional squared triangulation - // that spans from (0,0) to (1,1). It consists of two sub-domains. - // The left sub-domain spans from (0,0) to (0.525,1). The right - // sub-domain spans from (0.525,0) to (1,1). The left sub-domain has - // three times smaller elements compared to the right sub-domain. + // Let us now make our way to the higher-level functions of this program. + // + // The first of these functions creates a two dimensional square + // triangulation that spans from $(0,0)$ to $(1,1)$. It consists of + // two sub-domains. The left sub-domain spans from $(0,0)$ to + // $(0.525,1)$. The right sub-domain spans from $(0.525,0)$ to + // $(1,1)$. The left sub-domain has elements that are three times + // smaller compared to the ones for the right sub-domain. + // + // At non-matching interfaces, we need to provide different boundary + // IDs for the cells that make up the two parts (because, while they + // may be physically adjacent, they are not logically neighbors + // given that the faces of cells on both sides do not match, and the + // Triangulation class will therefore treat the interface between + // the two parts as a "boundary"). These boundary IDs have to differ + // because later on RemotePointEvaluation has to search for remote + // points for each face, that are defined in the same mesh (since we + // merge the mesh) but not on the same side of the non-matching + // interface. As a consequence, we declare at the top symbolic names + // for these boundary indicators, and ensure that we return a set + // with these values to the caller for later use. + // + // The actual mesh is then constructed by first constructing the + // left and right parts separately (setting material ids to zero and + // one, respectively), and using the appropriate boundary ids for + // all parts of the mesh. We then use + // GridGenerator::merge_triangulations() to combine them into one + // (non-matching) mesh. We have to pay attention that should the two + // sub-triangulations have vertices at the same locations, that they + // are not merged (connecting the two triangulations logically) + // since we want the interface to be an actual boundary. We achieve + // this by providing a tolerance of zero for the merge, see the + // documentation of the function + // GridGenerator::merge_triangulations(). template void build_non_matching_triangulation( Triangulation &tria, @@ -1239,19 +1290,12 @@ namespace Step89 { const double length = 1.0; - // At non-matching interfaces, we provide different boundary - // IDs. These boundary IDs have to differ because later on - // RemotePointEvaluation has to search for remote points for - // each face, that are defined in the same mesh (since we merge - // the mesh) but not on the same side of the non-matching interface. const types::boundary_id non_matching_id_left = 98; const types::boundary_id non_matching_id_right = 99; - // Provide this information to the caller. - non_matching_faces.insert(non_matching_id_left); - non_matching_faces.insert(non_matching_id_right); + non_matching_faces = {non_matching_id_left, non_matching_id_right}; - // Construct left part of mesh. + /* Construct left part of mesh. */ Triangulation tria_left; const unsigned int subdiv_left = 11; GridGenerator::subdivided_hyper_rectangle(tria_left, @@ -1259,12 +1303,8 @@ namespace Step89 {0.0, 0.0}, {0.525 * length, length}); - // The left part of the mesh has the material ID 0. for (const auto &cell : tria_left.active_cell_iterators()) cell->set_material_id(0); - - // The right face is non-matching. All other boundary IDs - // are set to 0. for (const auto &face : tria_left.active_face_iterators()) if (face->at_boundary()) { @@ -1273,7 +1313,7 @@ namespace Step89 face->set_boundary_id(non_matching_id_left); } - // Construct right part of mesh. + /* Construct right part of mesh. */ Triangulation tria_right; const unsigned int subdiv_right = 4; GridGenerator::subdivided_hyper_rectangle(tria_right, @@ -1281,12 +1321,8 @@ namespace Step89 {0.525 * length, 0.0}, {length, length}); - // The right part of the mesh has the material ID 1. for (const auto &cell : tria_right.active_cell_iterators()) cell->set_material_id(1); - - // The left face is non-matching. All other boundary IDs - // are set to 0. for (const auto &face : tria_right.active_face_iterators()) if (face->at_boundary()) { @@ -1295,26 +1331,29 @@ namespace Step89 face->set_boundary_id(non_matching_id_right); } - // Merge triangulations with tolerance 0 to ensure no vertices - // are merged, see the documentation of the function - // @c merge_triangulations(). + /* Merge triangulations. */ GridGenerator::merge_triangulations(tria_left, tria_right, tria, /*tolerance*/ 0., /*copy_manifold_ids*/ false, /*copy_boundary_ids*/ true); + + /* Refine the result. */ tria.refine_global(refinements); } - // @sect3{Set up and run point-to-point interpolation} - // - // The main purpose of this function is to fill a - // `FERemoteEvaluationCommunicator` object that is needed for point-to-point - // interpolation. Additionally, the corresponding remote evaluators are set up - // using this remote communicator. Eventually, the operators are handed to the - // time integrator that runs the simulation. + // @sect3{Set up and running of the two schemes} + + // @sect4{Setup and running of the point-to-point interpolation scheme} // + // We are now at the two functions that run the overall schemes (the + // point-to-point and the mortaring schemes). The first of these + // functions fills a FERemoteEvaluationCommunicator object that is + // needed for point-to-point interpolation. Additionally, the + // corresponding remote evaluators are set up using this remote + // communicator. Eventually, the operators are handed to the time + // integrator that runs the simulation. template void run_with_point_to_point_interpolation( const MatrixFree &matrix_free, @@ -1327,7 +1366,7 @@ namespace Step89 const auto &dof_handler = matrix_free.get_dof_handler(); const auto &tria = dof_handler.get_triangulation(); - // Communication objects know about the communication pattern. I.e., + // Communication objects know about the communication pattern. That is, // they know about the cells and quadrature points that have to be // evaluated at remote faces. This information is given via // RemotePointEvaluation. Additionally, the communication objects @@ -1346,21 +1385,22 @@ namespace Step89 // // For the standard case of point to point-to-point interpolation without // any heuristic we make use of the utility function - // @c compute_remote_communicator_faces_point_to_point_interpolation(). - // Please refer to this function to see how to manually set up the - // remote communicator from outside. - + // Utilities::compute_remote_communicator_faces_point_to_point_interpolation(). + // Please refer to the documentation of this function to see how to manually + // set up the remote communicator from outside. + // + // Among the inputs for the remote communicator we need a list of + // boundary ids for the non-matching faces, along with a function + // object for each boundary id that returns a vector of true/false + // flags in which exactly the vertices of cells of the + // triangulation are marked that have a face at the boundary id in + // question. std::vector< std::pair()>>> non_matching_faces_marked_vertices; - for (const auto &nm_face : non_matching_faces) { - // Sufficient lambda, that rules out all cells connected to the current - // side of the non-matching interface to avoid self intersections. - auto marked_vertices = [&]() { - // only search points at cells that are not connected to - // @c nm_face + auto marked_vertices = [&nm_face, &tria]() -> std::vector { std::vector mask(tria.n_vertices(), true); for (const auto &cell : tria.active_cell_iterators()) @@ -1383,22 +1423,20 @@ namespace Step89 // We are using point-to-point interpolation and can therefore // easily access the corresponding data at face batches. This - // is why we use a @c VectorizedArray as @c remote_value_type + // is why we use a @c VectorizedArray as @c remote_value_type: using remote_value_type = VectorizedArray; - // Set up FERemoteEvaluation object that accesses the pressure - // at remote faces. + // We then set up FERemoteEvaluation objects that access the + // pressure and velocity at remote faces, along with an object to + // describe cell-wise material data. const auto pressure_r = std::make_shared>( remote_communicator, dof_handler, /*first_selected_component*/ 0); - // Set up FERemoteEvaluation object that accesses the velocity - // at remote faces. const auto velocity_r = std::make_shared>( remote_communicator, dof_handler, /*first_selected_component*/ 1); - // Set up cell-wise material data. const auto material_data = std::make_shared>(matrix_free, materials); @@ -1415,9 +1453,14 @@ namespace Step89 matrix_free.get_dof_handler().get_triangulation(), /*first_selected_component*/ 0); + // If the domain is not homogeneous, i.e., if material parameters + // change from cell to cell, we initialize and fill DoF vectors + // that contain the material properties. Materials do not change + // during the simulation, therefore there is no need to ever + // compute the values after the first @c gather_evaluate() (at the + // end of the following block) again. if (!material_data->is_homogeneous()) { - // Initialize and fill DoF vectors that contain the materials. Vector c( matrix_free.get_dof_handler().get_triangulation().n_active_cells()); Vector rho( @@ -1433,21 +1476,26 @@ namespace Step89 materials.at(cell->material_id()).second; } - // Materials do not change during the simulation, therefore - // there is no need to precompute the values after - // the first @c gather_evaluate() again. c_r->gather_evaluate(c, EvaluationFlags::values); rho_r->gather_evaluate(rho, EvaluationFlags::values); } - // Set up inverse mass operator. + // Next, we set up the inverse mass operator and the acoustic + // operator. Using `remote_value_type=VectorizedArray` + // makes the operator use point-to-point interpolation. These two + // objects are then used to create a `RungeKutta2` object to + // perform the time integration. + // + // We also compute the maximum speed of sound, needed for the + // computation of the time-step size, and then run the time integrator. + // For the examples considered here, we found a limiting Courant number of + // $\mathrm{Cr}\approx 0.36$ to maintain stability. To ensure, the + // error of the temporal discretization is small, we use a considerably + // smaller Courant number of $0.2$. const auto inverse_mass_operator = std::make_shared>(matrix_free); - // Set up the acoustic operator. Using - // `remote_value_type=VectorizedArray` makes the operator use - // point-to-point interpolation. const auto acoustic_operator = std::make_shared>( matrix_free, @@ -1458,20 +1506,14 @@ namespace Step89 c_r, rho_r); - // Compute the the maximum speed of sound, needed for the computation of - // the time-step size. + RungeKutta2 time_integrator( + inverse_mass_operator, acoustic_operator); + double speed_of_sound_max = 0.0; for (const auto &mat : materials) speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first); - // Set up time integrator. - RungeKutta2 time_integrator( - inverse_mass_operator, acoustic_operator); - // For considered examples, we found a limiting Courant number of - // $\mathrm{Cr}\approx 0.36$ to maintain stability. To ensure, the - // error of the temporal discretization is small, we use a considerably - // smaller Courant number of $0.2$. time_integrator.run(matrix_free, /*Cr*/ 0.2, end_time, @@ -1480,14 +1522,13 @@ namespace Step89 vtk_prefix); } - // @sect3{Set up and run Nitsche-type mortaring} - // - // The main purpose of this function is to fill a - // `FERemoteEvaluationCommunicator` object that is needed for Nitsche-type - // mortaring. Additionally, the corresponding remote evaluators are set up - // using this remote communicator. Eventually, the operators are handed to the - // time integrator that runs the simulation. + // @sect4{Setup and running of the Nitsche-type mortaring scheme} // + // The alternative to the previous function is to use the mortaring + // scheme -- implemented in the following function. This function + // can only be run when deal.II is configured using CGAL (but the + // previous function can still be used without CGAL), and so errors + // out if CGAL is not available. template void run_with_nitsche_type_mortaring( const MatrixFree &matrix_free, @@ -1519,29 +1560,13 @@ namespace Step89 const auto &mapping = *matrix_free.get_mapping_info().mapping; const auto n_quadrature_pnts = matrix_free.get_quadrature().size(); - // In case of Nitsche-type mortaring a vector of pairs with cell iterator - // and face number is needed as communication object. - // @c FERemoteCommunicationObjectFaces is a container to store this - // information. - // - // For the standard case of Nitsche-type mortaring without - // any heuristic we make use of the utility function - // @c compute_remote_communicator_faces_nitsche_type_mortaring(). - // Please refer to this function to see how to manually set up the - // remote communicator from outside and how to reinit - // NonMatching::MappingInfo. - std::vector< std::pair()>>> non_matching_faces_marked_vertices; for (const auto &nm_face : non_matching_faces) { - // Sufficient lambda, that rules out all cells connected to the current - // side of the non-matching interface to avoid self intersections. auto marked_vertices = [&]() { - // only search points at cells that are not connected to - // @c nm_face std::vector mask(tria.n_vertices(), true); for (const auto &cell : tria.active_cell_iterators()) @@ -1558,7 +1583,10 @@ namespace Step89 std::make_pair(nm_face, marked_vertices)); } - // Quadrature points are arbitrarily distributed on each non-matching + // The only parts in this function that are functionally different + // from the previous one follow here. + // + // First, quadrature points are arbitrarily distributed on each non-matching // face. Therefore, we have to make use of FEFacePointEvaluation. // FEFacePointEvaluation needs NonMatching::MappingInfo to work at the // correct quadrature points that are in sync with used FERemoteEvaluation @@ -1589,11 +1617,14 @@ namespace Step89 0, nm_mapping_info.get()); - // Same as above but since quadrature points are aribtrarily distributed - // we have to consider each face in a batch separately and can not make - // use of @c VecorizedArray. + // Second, since quadrature points are arbitrarily distributed we + // have to consider each face in a batch separately and can not + // make use of @c VecorizedArray. using remote_value_type = Number; + // The rest of the code is then identical to what we had in the + // previous function (though it functions differently because of + // the difference in `remote_value_type`). const auto pressure_r = std::make_shared>( remote_communicator, dof_handler, /*first_selected_component*/ 0); @@ -1637,12 +1668,9 @@ namespace Step89 rho_r->gather_evaluate(rho, EvaluationFlags::values); } - // Set up inverse mass operator. const auto inverse_mass_operator = std::make_shared>(matrix_free); - // Set up the acoustic operator. Using `remote_value_type=Number` - // makes the operator use Nitsche-type mortaring. const auto acoustic_operator = std::make_shared>( matrix_free, @@ -1654,18 +1682,13 @@ namespace Step89 rho_r, nm_mapping_info); - // Compute the the maximum speed of sound, needed for the computation of - // the time-step size. + RungeKutta2 time_integrator( + inverse_mass_operator, acoustic_operator); + double speed_of_sound_max = 0.0; for (const auto &mat : materials) speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first); - - // Set up time integrator. - RungeKutta2 time_integrator( - inverse_mass_operator, acoustic_operator); - - // Run time loop with Courant number $0.2$. time_integrator.run(matrix_free, /*Cr*/ 0.2, end_time, @@ -1681,6 +1704,10 @@ namespace Step89 // // Finally, the `main()` function executes the different versions of handling // non-matching interfaces. +// +// Similar to step-87, the minimum requirement of this tutorial is MPI. +// The parallel::distributed::Triangulation class is used if deal.II is +// configured with p4est. Otherwise parallel::shared::Triangulation is used. int main(int argc, char *argv[]) { using namespace dealii; @@ -1696,11 +1723,6 @@ int main(int argc, char *argv[]) const unsigned int refinements = 1; const unsigned int degree = 3; - // Construct non-matching triangulation and fill non-matching boundary IDs. - - // Similar to step-87, the minimum requirement of this tutorial is MPI. - // The parallel::distributed::Triangulation class is used if deal.II is - // configured with p4est. Otherwise parallel::shared::Triangulation is used. #ifdef DEAL_II_WITH_P4EST parallel::distributed::Triangulation tria(MPI_COMM_WORLD); #else @@ -1717,8 +1739,6 @@ int main(int argc, char *argv[]) pcout << " - Refinement level: " << refinements << std::endl; pcout << " - Number of cells: " << tria.n_cells() << std::endl; - // Set up MatrixFree. - pcout << "Create DoFHandler..." << std::endl; DoFHandler dof_handler(tria); dof_handler.distribute_dofs(FESystem(FE_DGQ(degree) ^ (dim + 1))); @@ -1739,21 +1759,18 @@ int main(int argc, char *argv[]) MappingQ1(), dof_handler, constraints, QGauss(degree + 1), data); - //@sect4{Run vibrating membrane test case} + // @sect4{Run vibrating membrane test case} Homogeneous pressure} + // Dirichlet boundary conditions are applied for + // simplicity. Therefore, modes can not be chosen arbitrarily. pcout << "Run vibrating membrane test case..." << std::endl; - // Vibrating membrane test case: - // - // Homogeneous pressure DBCs are applied for simplicity. Therefore, - // modes can not be chosen arbitrarily. const double modes = 10.0; std::map> homogeneous_material; homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0); const auto initial_solution_membrane = Step89::InitialConditionVibratingMembrane(modes); + /* Run vibrating membrane test case using point-to-point interpolation: */ pcout << " - Point-to-point interpolation: " << std::endl; - // Run vibrating membrane test case using point-to-point interpolation: - Step89::run_with_point_to_point_interpolation( matrix_free, non_matching_faces, @@ -1763,8 +1780,8 @@ int main(int argc, char *argv[]) initial_solution_membrane, "vm-p2p"); + /* Run vibrating membrane test case using Nitsche-type mortaring: */ pcout << " - Nitsche-type mortaring: " << std::endl; - // Run vibrating membrane test case using Nitsche-type mortaring: Step89::run_with_nitsche_type_mortaring( matrix_free, non_matching_faces, @@ -1774,12 +1791,10 @@ int main(int argc, char *argv[]) initial_solution_membrane, "vm-nitsche"); - //@sect4{Run test case with in-homogeneous material} - pcout << "Run test case with in-homogeneous material..." << std::endl; - // In-homogeneous material test case: - // - // Run simple test case with in-homogeneous material and Nitsche-type + // @sect4{Run test case with inhomogeneous material} + // Run simple test case with inhomogeneous material and Nitsche-type // mortaring: + pcout << "Run test case with inhomogeneous material..." << std::endl; std::map> inhomogeneous_material; inhomogeneous_material[0] = std::make_pair(1.0, 1.0);