From: Ignacio Tomas (-EXP) Date: Wed, 1 Jan 2020 02:47:25 +0000 (-0700) Subject: update documentation, part III X-Git-Tag: v9.2.0-rc1~443^2~29 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8c598e9fb2018d9464728a1db09111abaeaaa634;p=dealii.git update documentation, part III --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 65aa4f51b3..08f2890226 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -531,7 +531,7 @@ MRREVIEWER = {Jose Luis Gracia}, DOI = {10.2307/2008211}, } -@inbook{Rainald2008, +@book{Rainald2008, author = {Lohner, Rainald}, publisher = {John Wiley & Sons, Ltd}, isbn = {9780470989746}, @@ -615,8 +615,6 @@ year = {2008}, journal = {Computer Methods in Applied Mechanics and Engineering} } - - % ------------------------------------ % References used elsewhere % ------------------------------------ diff --git a/examples/step-69/doc/intro.dox b/examples/step-69/doc/intro.dox index 1989403057..8af56dbf2a 100644 --- a/examples/step-69/doc/intro.dox +++ b/examples/step-69/doc/intro.dox @@ -265,7 +265,8 @@ Where - $d_{ij} := \max \{ \lambda_{\text{max}} (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}), \lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n}, - \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2} $ + \textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|_{\ell^2}$ if $i \not = j$ + - $d_{ii} = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}} d_{ij}$ - $\textbf{n}_{ij} = \frac{\mathbf{c}_{ij}}{ \|\mathbf{c}_{ij}\|_{\ell^2} }$ The definition of $\lambda_{\text{max}} (\mathbf{U},\mathbf{V}, @@ -359,8 +360,8 @@ cells but rather over all edges of the sparsity graph. no bilinear forms, no cell loops, and no quadrature) outside of the finite element community in the wieder CFD community. There is a rich history of application of this kind of schemes, also called edge-based or -graph-based finite element schemes (see for instance @cite -Rainald2008 for a historical overview). +graph-based finite element schemes (see for instance +@cite Rainald2008 for a historical overview). @todo Explain what to do for slip, dirichlet and do-nothing boundary conditions. diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index b169390033..95367d3c8d 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -276,17 +276,17 @@ namespace Step69 // @sect4{class InitialValues} // - // The class InitialValues only public data type is a - // function initial_state that computes the initial state of - // a given point and time. For the purpose of this example step we simply - // implement a homogeneous uniform flow field for which the direction and - // a 1D primitive state (density, velocity, pressure) are read from the - // parameter file. + // The class InitialValues's only public data attribute is a + // std::function initial_state that computes the initial + // state of a given point and time. For the purpose of this example step + // we simply implement a homogeneous uniform flow field for which the + // direction and a 1D primitive state (density, velocity, pressure) are + // read from the parameter file. // // Instead of implementing yet another setup() function we // use a callback function parse_parameters_callback that - // can be hooked up to corresponding signal of the ParameterAcceptor, - // ParameterAcceptor::parse_parameters_call_back.connect(...). + // can be hooked up to the corresponding signal + // ParameterAcceptor::parse_parameters_call_back. template class InitialValues : public ParameterAcceptor @@ -306,11 +306,18 @@ namespace Step69 }; // @sect4{class TimeStep} - - // Placeholder here + // + // With the OfflineData and ProblemDescription + // classes at hand we can now implement the explicit time-stepping scheme + // that was introduced in the discussion above. The main method of the + // TimeStep class is step(vector_type &U, double + // t). That takes a reference to a state vector U and + // a time point t as arguments, updates the state vector in + // place and returns the chosen step-size $\tau$. + // template - class TimeStep : public dealii::ParameterAcceptor + class TimeStep : public ParameterAcceptor { public: static constexpr unsigned int problem_dimension = @@ -319,12 +326,12 @@ namespace Step69 using rank1_type = typename ProblemDescription::rank1_type; using rank2_type = typename ProblemDescription::rank2_type; - typedef std::array, + typedef std::array, problem_dimension> vector_type; TimeStep(const MPI_Comm & mpi_communicator, - dealii::TimerOutput & computing_timer, + TimerOutput & computing_timer, const OfflineData & offline_data, const InitialValues &initial_values, const std::string & subsection = "TimeStep"); @@ -334,35 +341,35 @@ namespace Step69 double step(vector_type &U, double t); private: - const MPI_Comm & mpi_communicator; - dealii::TimerOutput &computing_timer; + const MPI_Comm &mpi_communicator; + TimerOutput & computing_timer; - dealii::SmartPointer> offline_data; - dealii::SmartPointer> initial_values; + SmartPointer> offline_data; + SmartPointer> initial_values; - dealii::SparseMatrix dij_matrix; + SparseMatrix dij_matrix; vector_type temp; double cfl_update; }; - // @sect4{Declaration of SchlierenPostprocessor class template} - - // At its core, the Schilieren class implements the class member + // @sect4{class SchlierenPostprocessor} + // + // At its core, the Schlieren class implements the class member // compute_schlieren. The main purpose of this class member - // is to compute auxiliary finite element field schlieren - // at each node, defined as + // is to compute an auxiliary finite element field + // schlieren, that is defined at each node by // \f[ \text{schlieren}[i] = e^{\beta \frac{ |\nabla r_i| - // - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } } \f] - // where $r$ in principle could be any scalar finite element field. - // The natural candidate is choosing $r := \rho$. Schlieren postprocessing - // is a standard methodology to enhance the contrast of the visualization - // inspired in actual X-ray and shadowgraphy experimental techniques of - // visualization. + // - \min_j |\nabla r_j| }{\max_j |\nabla r_j| - \min_j |\nabla r_j| } }, \f] + // where $r$ can in principle be any scalar quantitiy, in practice + // though, the density is a natural candidate, viz. $r := \rho$. + // Schlieren postprocessing is a standard method for enhancing the + // contrast of a visualization inspired by actual experimental X-ray and + // shadowgraphy techniques of visualization. template - class SchlierenPostprocessor : public dealii::ParameterAcceptor + class SchlierenPostprocessor : public ParameterAcceptor { public: static constexpr unsigned int problem_dimension = @@ -371,12 +378,11 @@ namespace Step69 using rank1_type = typename ProblemDescription::rank1_type; using vector_type = - std::array, - problem_dimension>; + std::array, problem_dimension>; SchlierenPostprocessor( const MPI_Comm & mpi_communicator, - dealii::TimerOutput & computing_timer, + TimerOutput & computing_timer, const OfflineData &offline_data, const std::string & subsection = "SchlierenPostprocessor"); @@ -384,26 +390,31 @@ namespace Step69 void compute_schlieren(const vector_type &U); - dealii::LinearAlgebra::distributed::Vector schlieren; + LinearAlgebra::distributed::Vector schlieren; private: - const MPI_Comm & mpi_communicator; - dealii::TimerOutput &computing_timer; + const MPI_Comm &mpi_communicator; + TimerOutput & computing_timer; - dealii::SmartPointer> offline_data; + SmartPointer> offline_data; - dealii::Vector r; + Vector r; unsigned int schlieren_index; double schlieren_beta; }; - // @sect4{Declaration of TimeLoop class template} - - // Placeholder here + // @sect4{class TimeLoop} + // + // Now, all that is left to do is to chain the methods implemented in the + // TimeStep, InitialValues, and + // SchlierenPostprocessor classes together. We do this in a + // separate class TimeLoop that contains an object of every + // class and again reads in a number of parameters with the help of the + // ParameterAcceptor class. template - class TimeLoop : public dealii::ParameterAcceptor + class TimeLoop : public ParameterAcceptor { public: using vector_type = typename TimeStep::vector_type; @@ -421,11 +432,11 @@ namespace Step69 unsigned int cycle, bool checkpoint = false); - const MPI_Comm & mpi_communicator; - std::ostringstream timer_output; - dealii::TimerOutput computing_timer; + const MPI_Comm & mpi_communicator; + std::ostringstream timer_output; + TimerOutput computing_timer; - dealii::ConditionalOStream pcout; + ConditionalOStream pcout; std::string base_name; double t_final; @@ -446,7 +457,7 @@ namespace Step69 vector_type output_vector; }; - // @sect3{Implementation of the classes in namespace Step69} + // @sect3{Class template implementations} // @sect4{Implementation of the members of the class Discretization} @@ -744,9 +755,9 @@ namespace Step69 // Now we define a collection of assembly utilities: // - CopyData: This will only be used to compute the off-line // data using WorkStream. It acts as a container: it is just a - // struct where WorkStream stores the local cell contributions. Note + // struct where WorkStream stores the local cell contributions. Note // it also contains a class member - // local_boundary_normal_map used to store the local + // local_boundary_normal_map used to store the local // contributions required to compute the normals at the boundary. // - get_entry: it reads the value stored at the entry // pointed by the iterator it of matrix. Here is @@ -775,10 +786,10 @@ namespace Step69 struct CopyData { bool is_artificial; - std::vector local_dof_indices; + std::vector local_dof_indices; typename OfflineData::BoundaryNormalMap local_boundary_normal_map; - dealii::FullMatrix cell_lumped_mass_matrix; - std::array, dim> cell_cij_matrix; + FullMatrix cell_lumped_mass_matrix; + std::array, dim> cell_cij_matrix; }; @@ -806,10 +817,10 @@ namespace Step69 template - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> + DEAL_II_ALWAYS_INLINE inline Tensor<1, k> gather_get_entry(const std::array &U, const T2 it) { - dealii::Tensor<1, k> result; + Tensor<1, k> result; for (unsigned int j = 0; j < k; ++j) result[j] = get_entry(U[j], it); return result; @@ -817,10 +828,10 @@ namespace Step69 template - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> + DEAL_II_ALWAYS_INLINE inline Tensor<1, k> gather(const std::array &U, const T2 i, const T3 l) { - dealii::Tensor<1, k> result; + Tensor<1, k> result; for (unsigned int j = 0; j < k; ++j) result[j] = U[j](i, l); return result; @@ -828,10 +839,10 @@ namespace Step69 template - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, k> - gather(const std::array &U, const T2 i) + DEAL_II_ALWAYS_INLINE inline Tensor<1, k> gather(const std::array &U, + const T2 i) { - dealii::Tensor<1, k> result; + Tensor<1, k> result; for (unsigned int j = 0; j < k; ++j) result[j] = U[j].local_element(i); return result; @@ -879,7 +890,7 @@ namespace Step69 // $\boldsymbol{\nu}_i := \sum_{T \in \text{supp}(\phi_i)} // \sum_{F \subset \partial T \cap \partial \Omega} // \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F}) - // \phi_i(\mathbf{x}_{q,F})$, here: $T$ denotes elements, + // \phi_i(\mathbf{x}_{q,F})$, here: $T$ denotes elements, // $\text{supp}(\phi_i)$ the support of the shape function $\phi_i$, // $F$ are faces of the element $T$, and $\mathbf{x}_{q,F}$ // are quadrature points on such face. @@ -994,8 +1005,8 @@ namespace Step69 /* Note that "normal" will only represent the contributions from one of the faces in the support of the shape function phi_j. So we cannot normalize this local - contribution right here, we have to take it "as is", store - it and pass it to the copy data routine. The proper + contribution right here, we have to take it "as is", store + it and pass it to the copy data routine. The proper normalization requires an additional loop on nodes.*/ Tensor<1, dim> normal; if (id == Boundary::slip) @@ -1094,7 +1105,7 @@ namespace Step69 // // We have the thread paralellization capability // parallel::apply_to_subranges that is somehow more general than the - // WorkStream framework. In particular, parallel::apply_to_subranges can + // WorkStream framework. In particular, parallel::apply_to_subranges can // be used for our node-loops. // This functionality requires four input arguments: // - A begin iterator: indices.begin() @@ -1108,15 +1119,15 @@ namespace Step69 // - Grainsize: minimum number of "elements" (in this case rows) processed // by each thread. We decided for a minimum of 4096 rows. // - // Here the indices.begin() and indices.end() + // Here the indices.begin() and indices.end() // iterators will represent an interval of "rows" - // in the sparsity graph/matrix. A minor caveat here is that the + // in the sparsity graph/matrix. A minor caveat here is that the // iterators supplied to // parallel::apply_to_subranges have to be random access iterators: // internally, apply_to_subranges will break the range defined by the - // indices.begin() and indices.end() iterators - // into subranges (we want to be able to read any entry in those - // subranges with constant complexity). In order to provide such + // indices.begin() and indices.end() iterators + // into subranges (we want to be able to read any entry in those + // subranges with constant complexity). In order to provide such // iterators we resort to boost::irange. // // We define the operation on_subranges to be @@ -1133,10 +1144,10 @@ namespace Step69 // argument required by std::for_each is the operation applied at each // column (a lambda expression in this case) of such row. We note that // because of the nature of the data that we want to modify (we want to - // modify entries of a entire row at a time) threads cannot collide - // attempting to write the same entry (we do not need a scheduler). This - // advantage appears to be a particular characteristic of edge-based finite - // element schemes when they are properly implemented. + // modify entries of a entire row at a time) threads cannot conflict + // attempting to read/write the same entry (we do not need a scheduler). + // This advantage appears to be a particular characteristic of + // edge-based finite element schemes when they are properly implemented. // // Finally, we normalize the vector stored in // OfflineData::BoundaryNormalMap. This operation has @@ -1187,8 +1198,8 @@ namespace Step69 on_subranges, 4096); - /* We normalize the normals at the boundary. This is not thread - parallelized. It just loops over the very few nodes that happen + /* We normalize the normals at the boundary. This is not thread + parallelized. It just loops over the very few nodes that happen to be at the boundary */ for (auto &it : boundary_normal_map) { @@ -1311,7 +1322,7 @@ namespace Step69 // At this point we are very much done with anything related to offline data. // - // Now we define the implementation of the utility + // Now we define the implementation of the utility // functions momentum, // internal_energy, pressure, // speed_of_sound, and f (the flux of the system). @@ -1319,10 +1330,10 @@ namespace Step69 // their names. template - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim> + DEAL_II_ALWAYS_INLINE inline Tensor<1, dim> ProblemDescription::momentum(const rank1_type &U) { - dealii::Tensor<1, dim> result; + Tensor<1, dim> result; std::copy(&U[1], &U[1 + dim], &result[0]); return result; } @@ -1377,40 +1388,40 @@ namespace Step69 } // Now we discuss the computation of $\lambda_{\text{max}} - // (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})$. Let's start - // by mentioning a thing or two about the actual computation of an estimate - // for maximum wavespeed of Riemann problem. In general, obtaining a sharp - // guaranteed upper-bound on the maximum wavespeed requires solving a - // quite expensive scalar nonlinear problem. In order to simplify the - // presentation we decided not to include such iterative scheme. Here we have - // taken the following shortcut: formulas (2.11) (3.7), (3.8) and (4.3) from - // - J-L Guermond, B. Popov, Fast estimation of - // the maximum wave speed in the Riemann problem for the Euler equations, + // (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})$. Let's start + // by mentioning a thing or two about the actual computation of an estimate + // for maximum wavespeed of Riemann problem. In general, obtaining a sharp + // guaranteed upper-bound on the maximum wavespeed requires solving a + // quite expensive scalar nonlinear problem. In order to simplify the + // presentation we decided not to include such iterative scheme. Here we have + // taken the following shortcut: formulas (2.11) (3.7), (3.8) and (4.3) from + // - J-L Guermond, B. Popov, Fast estimation of + // the maximum wave speed in the Riemann problem for the Euler equations, // JCP, 2016, // - // are enough to define a guaranteed upper bound on the maximum + // are enough to define a guaranteed upper bound on the maximum // wavespeed. This estimate is returned by the a call to the function // lambda_max_two_rarefaction. - // At its core the construction of such upper bound uses the - // so-called two-rarefaction approximation + // At its core the construction of such upper bound uses the + // so-called two-rarefaction approximation // for the intermediate pressure $p^*$, see for instance - // - Formula (4.46), page 128 in: E.Toro, Riemann Solvers and Numerical + // - Formula (4.46), page 128 in: E.Toro, Riemann Solvers and Numerical // Methods for Fluid Dynamics, 2009. // - // This estimate is in general very sharp and it would be enough to - // for this code. However, for some specific situations (in - // particular when one of states is close to vacuum conditions) such - // estimate will be very overly pessimistic. That's why we used a second - // estimate to avoid this degeneracy that will be invoked by a call to - // the function lambda_max_expansion. Finally we take the minimum + // This estimate is in general very sharp and it would be enough to + // for this code. However, for some specific situations (in + // particular when one of states is close to vacuum conditions) such + // estimate will be overly pessimistic. That's why we used a second + // estimate to avoid this degeneracy that will be invoked by a call to + // the function lambda_max_expansion. Finally we take the minimum // between both estimates inside the call to compute_lambda_max. // - // The analysis and derivation of sharp upper-bounds of maximum wavespeeds of - // Riemann problems is a very technical endeavor and we cannot include an + // The analysis and derivation of sharp upper-bounds of maximum wavespeeds of + // Riemann problems is a very technical endeavor and we cannot include an // advanced discussion about it in this tutorial. In this portion of the - // documentation we will limit ourselves to sketch the main functionality of - // these auxiliary functions and point to specific references/formulas in - // order to help the interested reader trace the + // documentation we will limit ourselves to sketch the main functionality of + // these auxiliary functions and point to specific references/formulas in + // order to help the interested reader trace the // source (and proper mathematical justification) of these ideas. // // The most important function here is compute_lambda_max @@ -1418,7 +1429,7 @@ namespace Step69 // - lambda_max_two_rarefaction // - lambda_max_expansion // - // The remaining functions + // The remaining functions // - riemann_data_from_state // - positive_part // - negative_part @@ -1432,9 +1443,9 @@ namespace Step69 template DEAL_II_ALWAYS_INLINE inline std::array riemann_data_from_state( const typename ProblemDescription::rank1_type U, - const dealii::Tensor<1, dim> & n_ij) + const Tensor<1, dim> & n_ij) { - dealii::Tensor<1, 3> projected_U; + Tensor<1, 3> projected_U; projected_U[0] = U[0]; const auto m = ProblemDescription::momentum(U); @@ -1518,9 +1529,9 @@ namespace Step69 }; - /* This estimate is, in general, not as sharp as the two-rarefaction - estimate. But it will save the day in the context of near vacuum - conditions when the two-rarefaction approximation will tend to + /* This estimate is, in general, not as sharp as the two-rarefaction + estimate. But it will save the day in the context of near vacuum + conditions when the two-rarefaction approximation will tend to exaggerate the maximum wave speed. */ DEAL_II_ALWAYS_INLINE inline double lambda_max_expansion(const std::array &riemann_data_i, @@ -1529,7 +1540,7 @@ namespace Step69 const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i; const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j; - /* Here the constant 5.0 multiplying the soundspeeds is NOT + /* Here the constant 5.0 multiplying the soundspeeds is NOT an ad-hoc constant. Do not play with it.*/ return std::max(std::abs(u_i), std::abs(u_j)) + 5. * std::max(a_i, a_j); } @@ -1540,10 +1551,9 @@ namespace Step69 // (\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij})$. template DEAL_II_ALWAYS_INLINE inline double - ProblemDescription::compute_lambda_max( - const rank1_type & U_i, - const rank1_type & U_j, - const dealii::Tensor<1, dim> &n_ij) + ProblemDescription::compute_lambda_max(const rank1_type & U_i, + const rank1_type & U_j, + const Tensor<1, dim> &n_ij) { const auto riemann_data_i = riemann_data_from_state(U_i, n_ij); const auto riemann_data_j = riemann_data_from_state(U_j, n_ij); @@ -1557,21 +1567,28 @@ namespace Step69 return std::min(lambda_1, lambda_2); } - // Placeholder here. + // Here component_names are just tags + // that we will use for the output. template <> - const std::array // - ProblemDescription<1>::component_names{"rho", "m", "E"}; + const std::array ProblemDescription<1>::component_names{"rho", + "m", + "E"}; template <> - const std::array // - ProblemDescription<2>::component_names{"rho", "m_1", "m_2", "E"}; + const std::array ProblemDescription<2>::component_names{"rho", + "m_1", + "m_2", + "E"}; template <> - const std::array // - ProblemDescription<3>::component_names{"rho", "m_1", "m_2", "m_3", "E"}; + const std::array ProblemDescription<3>::component_names{"rho", + "m_1", + "m_2", + "m_3", + "E"}; - // Placeholder here. + // Implementation of the constructor for the class InitialValues. template InitialValues::InitialValues(const std::string &subsection) @@ -1607,7 +1624,7 @@ namespace Step69 static constexpr auto gamma = ProblemDescription::gamma; const auto from_1d_state = - [=](const dealii::Tensor<1, 3, double> &state_1d) -> rank1_type { + [=](const Tensor<1, 3, double> &state_1d) -> rank1_type { const auto &rho = state_1d[0]; const auto &u = state_1d[1]; const auto &p = state_1d[2]; @@ -1622,16 +1639,16 @@ namespace Step69 return state; }; - initial_state = [=](const dealii::Point & /*point*/, double /*t*/) { + initial_state = [=](const Point & /*point*/, double /*t*/) { return from_1d_state(initial_1d_state); }; } - // Placeholder here. + // Implementation of the constructor for the class TimeStep. template TimeStep::TimeStep(const MPI_Comm & mpi_communicator, - dealii::TimerOutput & computing_timer, + TimerOutput & computing_timer, const OfflineData & offline_data, const InitialValues &initial_values, const std::string & subsection /*= "TimeStep"*/) @@ -1663,7 +1680,39 @@ namespace Step69 dij_matrix.reinit(sparsity); } - // Placeholder here. + // Implementation of "step" (to be called be + // TimeLoop::run()). We Start by computing the matrix + // $d_{ij}$. Pretty much all the ideas used to compute/store the entries + // of the matrix + // norm_matrix and the normalization of nij_matrix + // (described a few hundreds of lines above) are used here again. We use + // thread-parallel node-loops (again) via + // parallel::apply_to_subranges: therefore we have to + // define a "worker" on_subranges for this new task. + // + // We note here that $\int_{\Omega} \nabla \phi_j + // \phi_i \, \mathrm{d}\mathbf{x}= - \int_{\Omega} \nabla \phi_i \phi_j + // \, \mathrm{d}\mathbf{x}$ (or equivanlently $\mathbf{c}_{ij} = + // - \mathbf{c}_{ji}$) provided either $\mathbf{x}_i$ or $\mathbf{x}_j$ is a + // support point at the boundary. In such case we can check that: + // + // $\lambda_{\text{max}} (\mathbf{U}_i^{n}, \mathbf{U}_j^{n}, + // \textbf{n}_{ij}) = \lambda_{\text{max}} (\mathbf{U}_j^{n}, + // \mathbf{U}_i^{n}, + // \textbf{n}_{ji})$. + // + // However, if both support points $\mathbf{x}_i$ or $\mathbf{x}_j$ happen to + // lie on the boundary then the equality $\lambda_{\text{max}} + // (\mathbf{U}_i^{n}, \mathbf{U}_j^{n}, + // \textbf{n}_{ij}) = \lambda_{\text{max}} (\mathbf{U}_j^{n}, + // \mathbf{U}_i^{n}, + // \textbf{n}_{ji})$ is not necessarily true. The only mathematically + // safe solution for this dilemma is to compute both of them and take the + // largest one. + // + // The matrix $d_{ij}$ has to be symmetric by construction. Exploiting this + // natural constraint of the scheme we only compute the upper-triangular + // portion of it and then copy the result to the lower-triangular side. template double TimeStep::step(vector_type &U, double t) @@ -1687,6 +1736,7 @@ namespace Step69 { TimerOutput::Scope time(computing_timer, "time_step - 1 compute d_ij"); + /* Definition of the "worker" that computes the viscosity d_{ij} */ const auto on_subranges = [&](auto i1, const auto i2) { for (const auto i : boost::make_iterator_range(i1, i2)) { @@ -1697,6 +1747,8 @@ namespace Step69 { const auto j = jt->column(); + /* We compute only dij and later we copy this + entry into dji. */ if (j >= i) continue; @@ -1710,6 +1762,8 @@ namespace Step69 double d = norm * lambda_max; + /* If both support points happen to be at the boundary + we have to compute dji too and then take max(dij,dji) */ if (boundary_normal_map.count(i) != 0 && boundary_normal_map.count(j) != 0) { @@ -1735,6 +1789,33 @@ namespace Step69 4096); } /* End of the computation of the off-diagonal entries of dij_matrix */ + // So far the matrix dij_matrix contains the off-diagonal + // components. We still have to fill its diagonal entries defined as + // $d_{ii}^n = - \sum_{j \in \mathcal{I}(i)\backslash \{i\}} d_{ij}^n$. We + // use parallel::apply_to_subranges again in order to speed-up + // its computation. + + // While computing the $d_{ii}$'s we also record the largest admissible + // time-step, which is defined as + // + // $\tau_n := c_{\text{cfl}}\,\min_{ + // i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right)$ . + // + // We note that the operation $\min_{i \in \mathcal{V}}$ is intrinsically + // global, it operates on all nodes: first we would have to first take the + // $\min$ among all threads and finally take the $\min$ among all MPI + // processes. In the current implementation: + // - We do not take the $\min$ among threads: we simply define + // tau_max as + // std::atomic . The internal implementation of std::atomic + // will take care of resolving any possible conflict when more than + // one thread attempts read or write tau_max at the same time. + // - In order to take the min among all MPI process we use the utility + // Utilities::MPI::min. + + /* Atomic double in order to avoid any read/write conflict + * between threads */ std::atomic tau_max{std::numeric_limits::infinity()}; { @@ -1748,6 +1829,7 @@ namespace Step69 { double d_sum = 0.; + /* See the definition of dii */ for (auto jt = sparsity.begin(i); jt != sparsity.end(i); ++jt) { const auto j = jt->column(); @@ -1760,7 +1842,8 @@ namespace Step69 dij_matrix.diag_element(i) = d_sum; - const double mass = lumped_mass_matrix.diag_element(i); + const double mass = lumped_mass_matrix.diag_element(i); + /* See the definition of time-step constraint (CFL) */ const double tau = cfl_update * mass / (-2. * d_sum); tau_max_on_subrange = std::min(tau_max_on_subrange, tau); } @@ -1770,7 +1853,7 @@ namespace Step69 current_tau_max > tau_max_on_subrange && !tau_max.compare_exchange_weak(current_tau_max, tau_max_on_subrange)) ; - }; + }; /* End of definition of on_subranges */ parallel::apply_to_subranges(indices_relevant.begin(), indices_relevant.end(), @@ -1784,6 +1867,8 @@ namespace Step69 "do that. - We crashed.")); } /* End of the computation of the diagonal entries of dij_matrix */ + // Placeholder Here + { TimerOutput::Scope time(computing_timer, "time_step - 3 perform update"); @@ -1883,12 +1968,14 @@ namespace Step69 return tau_max; } /* End of TimeStep::step */ + + // Placeholder here. template SchlierenPostprocessor::SchlierenPostprocessor( const MPI_Comm & mpi_communicator, - dealii::TimerOutput & computing_timer, + TimerOutput & computing_timer, const OfflineData &offline_data, const std::string & subsection /*= "SchlierenPostprocessor"*/) : ParameterAcceptor(subsection) @@ -2074,9 +2161,9 @@ namespace Step69 namespace { - void print_head(dealii::ConditionalOStream &pcout, - std::string header, - std::string secondary = "") + void print_head(ConditionalOStream &pcout, + std::string header, + std::string secondary = "") { const auto header_size = header.size(); const auto padded_header = std::string((34 - header_size) / 2, ' ') + @@ -2136,7 +2223,7 @@ namespace Step69 const auto & triangulation = discretization.triangulation; const unsigned int i = triangulation.locally_owned_subdomain(); std::string name = base_name + "-checkpoint-" + - dealii::Utilities::int_to_string(i, 4) + ".archive"; + Utilities::int_to_string(i, 4) + ".archive"; std::ifstream file(name, std::ios::binary); boost::archive::binary_iarchive ia(file); @@ -2253,8 +2340,7 @@ namespace Step69 { const unsigned int i = triangulation.locally_owned_subdomain(); std::string name = base_name + "-checkpoint-" + - dealii::Utilities::int_to_string(i, 4) + - ".archive"; + Utilities::int_to_string(i, 4) + ".archive"; // FIXME: Refactor to Boost (this is C++17) // if (std::filesystem::exists(name)) @@ -2269,7 +2355,7 @@ namespace Step69 oa << it2; } - dealii::DataOut data_out; + DataOut data_out; data_out.attach_dof_handler(dof_handler); for (unsigned int i = 0; i < problem_dimension; ++i) @@ -2309,4 +2395,3 @@ int main(int argc, char *argv[]) time_loop.run(); } -