DgOperator.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <string>
8 #include <tuple>
9 #include <unordered_map>
10 #include <utility>
11 
12 #include "DataStructures/DataVector.hpp"
13 #include "DataStructures/SliceVariables.hpp"
17 #include "Domain/Structure/DirectionMap.hpp"
20 #include "Domain/Structure/IndexToSliceAt.hpp"
21 #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
22 #include "Elliptic/Systems/GetSourcesComputer.hpp"
23 #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
25 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
26 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
27 #include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleBoundaryData.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleMortarData.hpp"
30 #include "NumericalAlgorithms/LinearOperators/Divergence.tpp"
32 #include "NumericalAlgorithms/Spectral/Projection.hpp"
35 #include "Utilities/Gsl.hpp"
36 #include "Utilities/TMPL.hpp"
37 #include "Utilities/TaggedTuple.hpp"
38 
39 /*!
40  * \brief Functionality related to discontinuous Galerkin discretizations of
41  * elliptic equations
42  *
43  * The following is a brief overview of the elliptic DG schemes that are
44  * implemented here. A publication that describes the schemes in detail is in
45  * preparation and will be referenced here.
46  *
47  * The DG schemes apply to any elliptic PDE that can be formulated in
48  * first-order flux-form
49  *
50  * \f{equation}
51  * -\partial_i F^i_\alpha + S_\alpha = f_\alpha(x)
52  * \f}
53  *
54  * in terms of fluxes \f$F^i_\alpha\f$, sources \f$S_\alpha\f$ and fixed-sources
55  * \f$f_\alpha(x)\f$. The fluxes and sources are functionals of the "primal"
56  * system variables \f$u_A(x)\f$ and their corresponding "auxiliary" variables
57  * \f$v_A(x)\f$. The uppercase index enumerates the variables such that
58  * \f$v_A\f$ is the auxiliary variable corresponding to \f$u_A\f$. Greek indices
59  * enumerate _all_ variables. In documentation related to particular elliptic
60  * systems we generally use the canonical system-specific symbols for the fields
61  * in place of these indices. See the `Poisson::FirstOrderSystem` and the
62  * `Elasticity::FirstOrderSystem` for examples.
63  *
64  * The DG discretization of equations in this first-order form amounts to
65  * projecting the equations on the set of basis functions that we also use to
66  * represent the fields on the computational grid. The currently implemented DG
67  * operator uses Lagrange interpolating polynomials w.r.t.
68  * Legendre-Gauss-Lobatto collocation points as basis functions. Support for
69  * Legendre-Gauss collocation points can be added if needed. Skipping all
70  * further details here, the discretization results in a linear equation
71  * \f$A(u)=b\f$ over all grid points and primal variables. Solving the elliptic
72  * equations amounts to numerically inverting the DG operator \f$A\f$, typically
73  * without ever constructing the full matrix but by employing an iterative
74  * linear solver that repeatedly applies the DG operator to "test data". Note
75  * that the DG operator applies directly to the primal variables. Auxiliary
76  * variables are only computed temporarily and don't inflate the size of the
77  * operator. This means the DG operator essentially computes second derivatives
78  * of the primal variables, modified by the fluxes and sources of the system
79  * as well as by DG boundary corrections that couple grid points across element
80  * boundaries.
81  *
82  * \par Boundary corrections:
83  * In this implementation we employ the "internal penalty" DG scheme that
84  * couples grid points across nearest-neighbor elements through the fluxes:
85  *
86  * \f{align}
87  * (n_i F^i_v)^* &= \frac{1}{2} n_i \left(
88  * F^i_v(u^\mathrm{int}) +
89  * F^i_v(u^\mathrm{ext})\right) \\
90  * (n_i F^i_u)^* &= \frac{1}{2} n_i \left(
91  * F^i_u(\partial_j F^j_v(u^\mathrm{int}) - S_v(u^\mathrm{int})) +
92  * F^i_u(\partial_j F^j_v(u^\mathrm{ext}) - S_v(u^\mathrm{ext}))\right)
93  * - \sigma n_i \left(
94  * F^i_u(n_j F^j_v(u^\mathrm{int})) -
95  * F^i_u(n_j F^j_v(u^\mathrm{ext}))\right)
96  * \f}
97  *
98  * Note that \f$n_i\f$ denotes the face normal on the "interior" side of the
99  * element under consideration. We assume \f$n^\mathrm{ext}_i=-n_i\f$ in the
100  * implementation, i.e. face normals don't depend on the dynamic variables
101  * (which may be discontinuous on element faces). This is the case for the
102  * problems we are expecting to solve, because those will be on fixed background
103  * metrics (e.g. a conformal metric for the XCTS system). Numerically, the face
104  * normals on either side of a mortar may nonetheless be different because the
105  * two faces adjacent to the mortar may resolve them at different resolutions.
106  *
107  * Also note that the numerical fluxes intentionally don't depend on the
108  * auxiliary field values \f$v\f$. This property allows us to communicate data
109  * for both the primal and auxiliary boundary corrections together, instead of
110  * communicating them in two steps. If we were to resort to a two-step
111  * communication we could replace the divergence in \f$(n_i F^i_u)^*\f$ with
112  * \f$v\f$, which would result in a generalized "stabilized central flux" that
113  * is slightly less sparse than the internal penalty flux (see e.g.
114  * \cite HesthavenWarburton, section 7.2). We could also choose to ignore the
115  * fluxes in the penalty term, but preliminary tests suggest that this may hurt
116  * convergence.
117  *
118  * For a Poisson system (see `Poisson::FirstOrderSystem`) this numerical flux
119  * reduces to the standard internal penalty flux (see e.g.
120  * \cite HesthavenWarburton, section 7.2, or \cite Arnold2002):
121  *
122  * \f{align}
123  * (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left(
124  * u^\mathrm{int} + u^\mathrm{ext}\right) \\
125  * (n_i F^i_u)^* &= n_i v_i^* = \frac{1}{2} n_i \left(
126  * \partial_i u^\mathrm{int} + \partial_i u^\mathrm{ext}\right)
127  * - \sigma \left(u^\mathrm{int} - u^\mathrm{ext}\right)
128  * \f}
129  *
130  * where a sum over repeated indices is assumed, since the equation is
131  * formulated on a Euclidean geometry.
132  *
133  * The penalty factor \f$\sigma\f$ is responsible for removing zero eigenmodes
134  * and impacts the conditioning of the linear operator to be solved. See
135  * `elliptic::dg::penalty` for details. For the element size that goes into
136  * computing the penalty we choose
137  * \f$h=\frac{J_\mathrm{volume}}{J_\mathrm{face}}\f$, i.e. the ratio of Jacobi
138  * determinants from logical to inertial coordinates in the element volume and
139  * on the element face, both evaluated on the face (see \cite Vincent2019qpd).
140  * Since both \f$N_\mathrm{points}\f$ and \f$h\f$ can be different on either
141  * side of the element boundary we take the maximum of \f$N_\mathrm{points}\f$
142  * and the pointwise minimum of \f$h\f$ across the element boundary as is done
143  * in \cite Vincent2019qpd. Note that we use the number of points
144  * \f$N_\mathrm{points}\f$ where \cite Vincent2019qpd uses the polynomial degree
145  * \f$N_\mathrm{points} - 1\f$ because we found unstable configurations on
146  * curved meshes when using the polynomial degree. Optimizing the penalty on
147  * curved meshes is subject to further investigation.
148  */
149 namespace elliptic::dg {
150 
151 namespace Tags {
152 /// Number of grid points perpendicular to an element face. Used to compute
153 /// the penalty (see `elliptic::dg::penalty`).
155  using type = size_t;
156 };
157 
158 /// A measure of element size perpendicular to an element face. Used to compute
159 /// the penalty (see `elliptic::dg::penalty`).
160 struct ElementSize {
161  using type = Scalar<DataVector>;
162 };
163 
164 /// The quantity \f$n_i F^i_u(n_j F^j_v(u))\f$ where \f$F^i\f$ is the system
165 /// flux for primal variables \f$u\f$ and auxiliary variables \f$v\f$, and
166 /// \f$n_i\f$ is the face normal. This quantity is projected to mortars to
167 /// compute the jump term of the numerical flux.
168 template <typename Tag>
170  using type = typename Tag::type;
171  using tag = Tag;
172 };
173 } // namespace Tags
174 
175 /// Data that is projected to mortars and communicated across element
176 /// boundaries
177 template <typename PrimalFields, typename PrimalFluxes>
179  tmpl::append<db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFields>,
182  tmpl::list<Tags::ElementSize>>,
183  tmpl::list<Tags::PerpendicularNumPoints>>;
184 
185 /// Construct `elliptic::dg::BoundaryData` assuming the variable data on the
186 /// element is zero, and project it to the mortar.
187 template <typename PrimalMortarFields, typename PrimalMortarFluxes, size_t Dim>
190  const Direction<Dim>& direction, const Mesh<Dim>& mesh,
191  const Scalar<DataVector>& face_normal_magnitude,
192  const Mesh<Dim - 1>& mortar_mesh,
193  const ::dg::MortarSize<Dim - 1>& mortar_size) noexcept {
194  const auto face_mesh = mesh.slice_away(direction.dimension());
195  const size_t face_num_points = face_mesh.number_of_grid_points();
197  face_num_points};
198  boundary_data.field_data.initialize(face_num_points, 0.);
199  get<Tags::PerpendicularNumPoints>(boundary_data.extra_data) =
200  mesh.extents(direction.dimension());
201  get(get<Tags::ElementSize>(boundary_data.field_data)) =
202  2. / get(face_normal_magnitude);
204  ? boundary_data.project_to_mortar(face_mesh, mortar_mesh,
205  mortar_size)
206  : std::move(boundary_data);
207 }
208 
209 /// Boundary data on both sides of a mortar.
210 ///
211 /// \note This is a struct (as opposed to a type alias) so it can be used to
212 /// deduce the template parameters
213 template <typename TemporalId, typename PrimalFields, typename PrimalFluxes>
215  : ::dg::SimpleMortarData<TemporalId,
216  BoundaryData<PrimalFields, PrimalFluxes>,
217  BoundaryData<PrimalFields, PrimalFluxes>> {};
218 
219 namespace Tags {
220 /// Holds `elliptic::dg::MortarData`, i.e. boundary data on both sides of a
221 /// mortar
222 template <typename TemporalId, typename PrimalFields, typename PrimalFluxes>
225 };
226 } // namespace Tags
227 
228 namespace detail {
229 template <typename System, bool Linearized,
230  typename PrimalFields = typename System::primal_fields,
231  typename AuxiliaryFields = typename System::auxiliary_fields,
232  typename PrimalFluxes = typename System::primal_fluxes,
233  typename AuxiliaryFluxes = typename System::auxiliary_fluxes>
234 struct DgOperatorImpl;
235 
236 template <typename System, bool Linearized, typename... PrimalFields,
237  typename... AuxiliaryFields, typename... PrimalFluxes,
238  typename... AuxiliaryFluxes>
239 struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
240  tmpl::list<AuxiliaryFields...>,
241  tmpl::list<PrimalFluxes...>,
242  tmpl::list<AuxiliaryFluxes...>> {
243  static constexpr size_t Dim = System::volume_dim;
244  using FluxesComputer = typename System::fluxes_computer;
246 
247  struct AllDirections {
248  bool operator()(const Direction<Dim>& /*unused*/) const noexcept {
249  return true;
250  }
251  };
252 
253  static constexpr auto full_mortar_size =
254  make_array<Dim - 1>(Spectral::MortarSize::Full);
255 
256  template <bool DataIsZero, typename... PrimalVars, typename... AuxiliaryVars,
257  typename... PrimalFluxesVars, typename... AuxiliaryFluxesVars,
258  typename... PrimalMortarVars, typename... PrimalMortarFluxes,
259  typename TemporalId, typename ApplyBoundaryCondition,
260  typename... FluxesArgs, typename... SourcesArgs,
261  typename DirectionsPredicate = AllDirections>
262  static void prepare_mortar_data(
263  const gsl::not_null<Variables<tmpl::list<AuxiliaryVars...>>*>
264  auxiliary_vars,
265  const gsl::not_null<Variables<tmpl::list<AuxiliaryFluxesVars...>>*>
266  auxiliary_fluxes,
267  const gsl::not_null<Variables<tmpl::list<PrimalFluxesVars...>>*>
268  primal_fluxes,
270  Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
271  tmpl::list<PrimalMortarFluxes...>>>*>
272  all_mortar_data,
273  const Variables<tmpl::list<PrimalVars...>>& primal_vars,
274  const Element<Dim>& element, const Mesh<Dim>& mesh,
275  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
276  inv_jacobian,
277  const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
278  const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
279  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
280  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
281  const TemporalId& temporal_id,
282  const ApplyBoundaryCondition& apply_boundary_condition,
283  const std::tuple<FluxesArgs...>& fluxes_args,
284  const std::tuple<SourcesArgs...>& sources_args,
285  const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
286  const DirectionsPredicate& directions_predicate =
287  AllDirections{}) noexcept {
288  static_assert(
289  sizeof...(PrimalVars) == sizeof...(PrimalFields) and
290  sizeof...(AuxiliaryVars) == sizeof...(AuxiliaryFields) and
291  sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
292  sizeof...(AuxiliaryFluxesVars) == sizeof...(AuxiliaryFluxes),
293  "The number of variables must match the number of system fields.");
294  static_assert(
295  (std::is_same_v<typename PrimalVars::type,
296  typename PrimalFields::type> and
297  ...) and
298  (std::is_same_v<typename AuxiliaryVars::type,
299  typename AuxiliaryFields::type> and
300  ...) and
301  (std::is_same_v<typename PrimalFluxesVars::type,
302  typename PrimalFluxes::type> and
303  ...) and
304  (std::is_same_v<typename AuxiliaryFluxesVars::type,
305  typename AuxiliaryFluxes::type> and
306  ...),
307  "The variables must have the same tensor types as the system fields.");
308 #ifdef SPECTRE_DEBUG
309  for (size_t d = 0; d < Dim; ++d) {
310  ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
311  mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
312  "The elliptic DG operator is currently only implemented for "
313  "Legendre-Gauss-Lobatto grids. Found basis '"
314  << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
315  << "' in dimension " << d << ".");
316  }
317 #endif // SPECTRE_DEBUG
318  const size_t num_points = mesh.number_of_grid_points();
319 
320  // This function and the one below allocate various Variables to compute
321  // intermediate quantities. It could be a performance optimization to reduce
322  // the number of these allocations and/or move some of the memory buffers
323  // into the DataBox to keep them around permanently. The latter should be
324  // informed by profiling.
325 
326  // Compute the auxiliary variables, and from those the primal fluxes. The
327  // auxiliary variables are the variables `v` in the auxiliary equations
328  // `-div(F_v(u)) + S_v + v = 0` where `F_v` and `S_v` are the auxiliary
329  // fluxes and sources. From the auxiliary variables we compute the primal
330  // fluxes `F_u(v)` for the primal equation `-div(F_u(v)) + S_u = f(x)`. Note
331  // that before taking the second derivative, boundary corrections from the
332  // first derivative have to be added to `F_u(v)`. Therefore the second
333  // derivative is taken after a communication break in the `apply_operator`
334  // function below.
335  if constexpr (DataIsZero) {
336  (void)auxiliary_vars;
337  (void)auxiliary_fluxes;
338  primal_fluxes->initialize(num_points, 0.);
339  } else {
340  // Compute the auxiliary fluxes
341  auxiliary_fluxes->initialize(num_points);
342  std::apply(
343  [&auxiliary_fluxes,
344  &primal_vars](const auto&... expanded_fluxes_args) noexcept {
345  FluxesComputer::apply(
346  make_not_null(&get<AuxiliaryFluxesVars>(*auxiliary_fluxes))...,
347  expanded_fluxes_args..., get<PrimalVars>(primal_vars)...);
348  },
349  fluxes_args);
350  divergence(auxiliary_vars, *auxiliary_fluxes, mesh, inv_jacobian);
351  // Subtract the sources. Invert signs because the sources computers _add_
352  // the sources, i.e. they compute the `+ S_v` term defined above.
353  *auxiliary_vars *= -1.;
354  std::apply(
355  [&auxiliary_vars,
356  &primal_vars](const auto&... expanded_sources_args) noexcept {
357  SourcesComputer::apply(
358  make_not_null(&get<AuxiliaryVars>(*auxiliary_vars))...,
359  expanded_sources_args..., get<PrimalVars>(primal_vars)...);
360  },
361  sources_args);
362  *auxiliary_vars *= -1.;
363  // Compute the primal fluxes
364  primal_fluxes->initialize(num_points);
365  std::apply(
366  [&primal_fluxes,
367  &auxiliary_vars](const auto&... expanded_fluxes_args) noexcept {
368  FluxesComputer::apply(
369  make_not_null(&get<PrimalFluxesVars>(*primal_fluxes))...,
370  expanded_fluxes_args...,
371  get<AuxiliaryVars>(*auxiliary_vars)...);
372  },
373  fluxes_args);
374  }
375 
376  // Populate the mortar data on this element's side of the boundary so it's
377  // ready to be sent to neighbors.
378  for (const auto& direction : [&element]() noexcept -> const auto& {
379  if constexpr (DataIsZero) {
380  // Skipping internal boundaries for zero data because they won't
381  // contribute boundary corrections anyway.
382  return element.external_boundaries();
383  } else {
384  (void)element;
386  };
387  }()) {
388  if (not directions_predicate(direction)) {
389  continue;
390  }
391  const bool is_internal = element.neighbors().contains(direction);
392  const auto face_mesh = mesh.slice_away(direction.dimension());
393  const size_t face_num_points = face_mesh.number_of_grid_points();
394  const auto& face_normal = face_normals.at(direction);
395  const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
396  const auto& fluxes_args_on_face = fluxes_args_on_faces.at(direction);
397  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
398  Variables<tmpl::list<PrimalFluxesVars..., AuxiliaryFluxesVars...>>
399  fluxes_on_face{face_num_points};
400  Variables<tmpl::list<::Tags::NormalDotFlux<AuxiliaryFields>...>>
401  n_dot_aux_fluxes{face_num_points};
402  BoundaryData<tmpl::list<PrimalMortarVars...>,
403  tmpl::list<PrimalMortarFluxes...>>
404  boundary_data{face_num_points};
405  if constexpr (DataIsZero) {
406  // Just setting all boundary field data to zero. Variable-independent
407  // data such as the element size will be set below.
408  boundary_data.field_data.initialize(face_num_points, 0.);
409  } else {
410  // Compute F_u(n.F_v)
411  fluxes_on_face.assign_subset(
412  data_on_slice(*auxiliary_fluxes, mesh.extents(),
413  direction.dimension(), slice_index));
414  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
416  &get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
417  face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
418  std::apply(
419  [&boundary_data, &n_dot_aux_fluxes](
420  const auto&... expanded_fluxes_args_on_face) noexcept {
421  FluxesComputer::apply(
423  boundary_data.field_data))...,
424  expanded_fluxes_args_on_face...,
426  n_dot_aux_fluxes)...);
427  },
428  fluxes_args_on_face);
429  // Compute n.F_u
430  //
431  // For the internal penalty flux we can already slice the n.F_u to the
432  // faces at this point, before the boundary corrections have been added
433  // to the auxiliary variables. The reason is essentially that the
434  // internal penalty flux uses average(grad(u)) in place of average(v),
435  // i.e. the raw primal field derivatives without boundary corrections.
436  fluxes_on_face.assign_subset(
437  data_on_slice(*primal_fluxes, mesh.extents(), direction.dimension(),
438  slice_index));
439  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
441  boundary_data.field_data)),
442  face_normal, get<PrimalFluxesVars>(fluxes_on_face)));
443  // Compute n.F_u(n.F_v) for jump term, re-using the memory buffer from
444  // above
445  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
446  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
447  boundary_data.field_data)),
448  face_normal,
450  boundary_data.field_data)));
451  }
452 
453  // Collect the remaining data that's needed on both sides of the boundary
454  // These are actually constant throughout the solve, so a performance
455  // optimization could be to store them in the DataBox.
456  get<Tags::PerpendicularNumPoints>(boundary_data.extra_data) =
457  mesh.extents(direction.dimension());
458  get(get<Tags::ElementSize>(boundary_data.field_data)) =
459  2. / get(face_normal_magnitude);
460 
461  if (is_internal) {
462  if constexpr (not DataIsZero) {
463  // Project boundary data on internal faces to mortars
464  for (const auto& neighbor_id : element.neighbors().at(direction)) {
465  const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
466  const auto& mortar_mesh = all_mortar_meshes.at(mortar_id);
467  const auto& mortar_size = all_mortar_sizes.at(mortar_id);
468  auto projected_boundary_data =
470  ? boundary_data.project_to_mortar(face_mesh, mortar_mesh,
471  mortar_size)
472  : std::move(boundary_data);
473  (*all_mortar_data)[mortar_id].local_insert(
474  temporal_id, std::move(projected_boundary_data));
475  }
476  }
477  } else {
478  // No need to do projections on external boundaries
479  const ::dg::MortarId<Dim> mortar_id{
481  (*all_mortar_data)[mortar_id].local_insert(temporal_id, boundary_data);
482 
483  // -------------------------
484  // Apply boundary conditions
485  // -------------------------
486  //
487  // To apply boundary conditions we fill the boundary data with
488  // "exterior" or "ghost" data and set it as remote mortar data, so
489  // external boundaries behave just like internal boundaries when
490  // applying boundary corrections.
491  //
492  // The `apply_boundary_conditions` invocable is expected to impose
493  // boundary conditions by modifying the fields and fluxes that are
494  // passed by reference. Dirichlet-type boundary conditions are imposed
495  // by modifying the fields, and Neumann-type boundary conditions are
496  // imposed by modifying the interior n.F_u. Note that all data passed to
497  // the boundary conditions is taken from the "interior" side of the
498  // boundary, i.e. with a normal vector that points _out_ of the
499  // computational domain.
500  auto dirichlet_vars = data_on_slice(primal_vars, mesh.extents(),
501  direction.dimension(), slice_index);
503  direction, make_not_null(&get<PrimalVars>(dirichlet_vars))...,
505  boundary_data.field_data))...);
506 
507  // The n.F_u (Neumann-type conditions) are done, but we have to compute
508  // the n.F_v (Dirichlet-type conditions) from the Dirichlet fields. We
509  // re-use the memory buffer from above.
510  std::apply(
511  [&fluxes_on_face, &dirichlet_vars](
512  const auto&... expanded_fluxes_args_on_face) noexcept {
513  FluxesComputer::apply(
514  make_not_null(&get<AuxiliaryFluxesVars>(fluxes_on_face))...,
515  expanded_fluxes_args_on_face...,
516  get<PrimalVars>(dirichlet_vars)...);
517  },
518  fluxes_args_on_face);
519  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
521  &get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
522  face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
523  std::apply(
524  [&boundary_data, &n_dot_aux_fluxes](
525  const auto&... expanded_fluxes_args_on_face) noexcept {
526  FluxesComputer::apply(
528  boundary_data.field_data))...,
529  expanded_fluxes_args_on_face...,
531  n_dot_aux_fluxes)...);
532  },
533  fluxes_args_on_face);
534 
535  // Invert the sign of the fluxes to account for the inverted normal on
536  // exterior faces. Also multiply by 2 and add the interior fluxes to
537  // impose the boundary conditions on the _average_ instead of just
538  // setting the fields on the exterior.
539  const auto invert_sign_and_impose_on_average =
540  [](const auto exterior_n_dot_flux,
541  const auto& interior_n_dot_flux) {
542  for (size_t i = 0; i < interior_n_dot_flux.size(); ++i) {
543  (*exterior_n_dot_flux)[i] *= -2.;
544  (*exterior_n_dot_flux)[i] += interior_n_dot_flux[i];
545  }
546  };
547  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
549  boundary_data.field_data)),
551  all_mortar_data->at(mortar_id)
552  .local_data(temporal_id)
553  .field_data)));
554  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
556  boundary_data.field_data)),
558  all_mortar_data->at(mortar_id)
559  .local_data(temporal_id)
560  .field_data)));
561 
562  // Compute n.F_u(n.F_v) for jump term
563  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
564  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
565  boundary_data.field_data)),
566  face_normal,
568  boundary_data.field_data)));
569  const auto invert_sign = [](const auto exterior_n_dot_flux) {
570  for (size_t i = 0; i < exterior_n_dot_flux->size(); ++i) {
571  (*exterior_n_dot_flux)[i] *= -1.;
572  }
573  };
574  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign(
575  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
576  boundary_data.field_data))));
577 
578  // Store the exterior boundary data on the mortar
579  all_mortar_data->at(mortar_id).remote_insert(temporal_id,
580  std::move(boundary_data));
581  } // if (is_internal)
582  } // loop directions
583  }
584 
585  // --- This is essentially a break to communicate the mortar data ---
586 
587  template <bool DataIsZero, typename... OperatorTags, typename... PrimalVars,
588  typename... PrimalFluxesVars, typename... PrimalMortarVars,
589  typename... PrimalMortarFluxes, typename TemporalId,
590  typename... FluxesArgs, typename... SourcesArgs,
591  typename DirectionsPredicate = AllDirections>
592  static void apply_operator(
593  const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
594  operator_applied_to_vars,
596  Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
597  tmpl::list<PrimalMortarFluxes...>>>*>
598  all_mortar_data,
599  const Variables<tmpl::list<PrimalVars...>>& primal_vars,
600  // Taking the primal fluxes computed in the `prepare_mortar_data` function
601  // by const-ref here because other code might use them and so we don't
602  // want to modify them by adding boundary corrections. E.g. linearized
603  // sources use the nonlinear fields and fluxes as background fields.
604  const Variables<tmpl::list<PrimalFluxesVars...>>& primal_fluxes,
605  const Mesh<Dim>& mesh,
606  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
607  inv_jacobian,
608  const Scalar<DataVector>& det_inv_jacobian,
609  const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
610  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
611  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
612  const double penalty_parameter, const bool massive,
613  const TemporalId& temporal_id,
614  const std::tuple<SourcesArgs...>& sources_args,
615  const DirectionsPredicate& directions_predicate =
616  AllDirections{}) noexcept {
617  static_assert(
618  sizeof...(PrimalVars) == sizeof...(PrimalFields) and
619  sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
620  sizeof...(PrimalMortarVars) == sizeof...(PrimalFields) and
621  sizeof...(PrimalMortarFluxes) == sizeof...(PrimalFluxes) and
622  sizeof...(OperatorTags) == sizeof...(PrimalFields),
623  "The number of variables must match the number of system fields.");
624  static_assert(
625  (std::is_same_v<typename PrimalVars::type,
626  typename PrimalFields::type> and
627  ...) and
628  (std::is_same_v<typename PrimalFluxesVars::type,
629  typename PrimalFluxes::type> and
630  ...) and
631  (std::is_same_v<typename PrimalMortarVars::type,
632  typename PrimalFields::type> and
633  ...) and
634  (std::is_same_v<typename PrimalMortarFluxes::type,
635  typename PrimalFluxes::type> and
636  ...) and
637  (std::is_same_v<typename OperatorTags::type,
638  typename PrimalFields::type> and
639  ...),
640  "The variables must have the same tensor types as the system fields.");
641 #ifdef SPECTRE_DEBUG
642  for (size_t d = 0; d < Dim; ++d) {
643  ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
644  mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
645  "The elliptic DG operator is currently only implemented for "
646  "Legendre-Gauss-Lobatto grids. Found basis '"
647  << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
648  << "' in dimension " << d << ".");
649  }
650 #endif // SPECTRE_DEBUG
651 
652  // This function and the one above allocate various Variables to compute
653  // intermediate quantities. It could be a performance optimization to reduce
654  // the number of these allocations and/or move some of the memory buffers
655  // into the DataBox to keep them around permanently. The latter should be
656  // informed by profiling.
657 
658  // Add boundary corrections to the auxiliary variables _before_ computing
659  // the second derivative. This is called the "flux" formulation. It is
660  // equivalent to discretizing the system in first-order form, i.e. treating
661  // the primal and auxiliary variables on the same footing, and then taking a
662  // Schur complement of the operator. The Schur complement is possible
663  // because the auxiliary equations are essentially the definition of the
664  // auxiliary variables and can therefore always be solved for the auxiliary
665  // variables by just inverting the mass matrix. This Schur-complement
666  // formulation avoids inflating the DG operator with the DOFs of the
667  // auxiliary variables. In this form it is very similar to the "primal"
668  // formulation where we get rid of the auxiliary variables through a DG
669  // theorem and thus add the auxiliary boundary corrections _after_ computing
670  // the second derivative. This involves a slightly different lifting
671  // operation with differentiation matrices, which we avoid to implement for
672  // now by using the flux-formulation.
673  auto primal_fluxes_corrected = primal_fluxes;
674  for (const auto& [mortar_id, mortar_data] : *all_mortar_data) {
675  const auto& [direction, neighbor_id] = mortar_id;
676  const bool is_internal =
677  (neighbor_id != ElementId<Dim>::external_boundary_id());
678  if constexpr (DataIsZero) {
679  if (is_internal) {
680  continue;
681  }
682  }
683  if (not directions_predicate(direction)) {
684  continue;
685  }
686  const auto face_mesh = mesh.slice_away(direction.dimension());
687  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
688  const auto& local_data = mortar_data.local_data(temporal_id);
689  const auto& remote_data = mortar_data.remote_data(temporal_id);
690  const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
691  const auto& mortar_mesh =
692  is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
693  const auto& mortar_size =
694  is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
695 
696  // This is the _strong_ auxiliary boundary correction avg(n.F_v) - n.F_v
697  auto auxiliary_boundary_corrections_on_mortar =
698  Variables<tmpl::list<PrimalMortarFluxes...>>(
699  local_data.field_data.template extract_subset<
700  tmpl::list<::Tags::NormalDotFlux<PrimalMortarFluxes>...>>());
701  const auto add_remote_contribution = [](auto& lhs,
702  const auto& rhs) noexcept {
703  for (size_t i = 0; i < lhs.size(); ++i) {
704  lhs[i] += rhs[i];
705  }
706  };
707  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_contribution(
708  get<PrimalMortarFluxes>(auxiliary_boundary_corrections_on_mortar),
710  remote_data.field_data)));
711  auxiliary_boundary_corrections_on_mortar *= -0.5;
712 
713  // Project from the mortar back down to the face if needed
714  auto auxiliary_boundary_corrections =
717  auxiliary_boundary_corrections_on_mortar, face_mesh,
719  : std::move(auxiliary_boundary_corrections_on_mortar);
720 
721  // Lift the boundary correction to the volume, but still only provide the
722  // data only on the face because it is zero everywhere else. This is the
723  // "massless" lifting operation, i.e. it involves an inverse mass matrix.
724  // The mass matrix is diagonally approximated ("mass lumping") so it
725  // reduces to a division by quadrature weights.
726  ::dg::lift_flux(make_not_null(&auxiliary_boundary_corrections),
727  mesh.extents(direction.dimension()),
728  face_normal_magnitude);
729  // The `dg::lift_flux` function contains an extra minus sign
730  auxiliary_boundary_corrections *= -1.;
731 
732  // Add the boundary corrections to the auxiliary variables
733  add_slice_to_data(make_not_null(&primal_fluxes_corrected),
734  auxiliary_boundary_corrections, mesh.extents(),
735  direction.dimension(), slice_index);
736  } // apply auxiliary boundary corrections on all mortars
737 
738  // Compute the primal equation, i.e. the actual DG operator, by taking the
739  // second derivative: -div(F_u(v)) + S_u = f(x)
740  divergence(operator_applied_to_vars, primal_fluxes_corrected, mesh,
741  inv_jacobian);
742  // This is the sign flip that makes the operator _minus_ the Laplacian for a
743  // Poisson system
744  *operator_applied_to_vars *= -1.;
745  // Using the non-boundary-corrected primal fluxes to compute sources here
746  std::apply(
747  [&operator_applied_to_vars, &primal_vars,
748  &primal_fluxes](const auto&... expanded_sources_args) noexcept {
749  SourcesComputer::apply(
750  make_not_null(&get<OperatorTags>(*operator_applied_to_vars))...,
751  expanded_sources_args..., get<PrimalVars>(primal_vars)...,
752  get<PrimalFluxesVars>(primal_fluxes)...);
753  },
754  sources_args);
755 
756  // Add boundary corrections to primal equation
757  for (auto& [mortar_id, mortar_data] : *all_mortar_data) {
758  const auto& [direction, neighbor_id] = mortar_id;
759  const bool is_internal =
760  (neighbor_id != ElementId<Dim>::external_boundary_id());
761  if constexpr (DataIsZero) {
762  if (is_internal) {
763  continue;
764  }
765  }
766  const auto face_mesh = mesh.slice_away(direction.dimension());
767  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
768  const auto [local_data, remote_data] = mortar_data.extract();
769  const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
770  const auto& mortar_mesh =
771  is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
772  const auto& mortar_size =
773  is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
774 
775  // This is the _strong_ primal boundary correction avg(n.F_u) - penalty *
776  // jump(n.F_v) - n.F_u. Note that the "internal penalty" numerical flux
777  // (as opposed to the LLF flux) uses the raw field derivatives without
778  // boundary corrections in the average, which is why we can communicate
779  // the data so early together with the auxiliary boundary data. In this
780  // case the penalty needs to include a factor N_points^2 / h (see the
781  // `penalty` function).
782  const auto penalty = elliptic::dg::penalty(
783  min(get(get<Tags::ElementSize>(local_data.field_data)),
784  get(get<Tags::ElementSize>(remote_data.field_data))),
785  std::max(get<Tags::PerpendicularNumPoints>(local_data.extra_data),
786  get<Tags::PerpendicularNumPoints>(remote_data.extra_data)),
787  penalty_parameter);
788  // Start with the penalty term
789  auto primal_boundary_corrections_on_mortar =
790  Variables<tmpl::list<PrimalMortarVars...>>(
791  local_data.field_data.template extract_subset<tmpl::list<
792  Tags::NormalDotFluxForJump<PrimalMortarVars>...>>());
793  const auto add_remote_jump_contribution = [](auto& lhs,
794  const auto& rhs) noexcept {
795  for (size_t i = 0; i < lhs.size(); ++i) {
796  lhs[i] -= rhs[i];
797  }
798  };
799  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_jump_contribution(
800  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
801  get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
802  remote_data.field_data)));
803  primal_boundary_corrections_on_mortar *= penalty;
804  const auto add_remote_avg_contribution = [](auto& lhs,
805  const auto& rhs) noexcept {
806  for (size_t i = 0; i < lhs.size(); ++i) {
807  lhs[i] += 0.5 * rhs[i];
808  }
809  };
810  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
811  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
812  get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data)));
813  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
814  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
816  remote_data.field_data)));
817  primal_boundary_corrections_on_mortar *= -1.;
818 
819  // Project from the mortar back down to the face if needed, lift and add
820  // to operator. See auxiliary boundary corrections above for details.
821  auto primal_boundary_corrections =
823  ? ::dg::project_from_mortar(primal_boundary_corrections_on_mortar,
824  face_mesh, mortar_mesh, mortar_size)
825  : std::move(primal_boundary_corrections_on_mortar);
826  ::dg::lift_flux(make_not_null(&primal_boundary_corrections),
827  mesh.extents(direction.dimension()),
828  face_normal_magnitude);
829  add_slice_to_data(operator_applied_to_vars, primal_boundary_corrections,
830  mesh.extents(), direction.dimension(), slice_index);
831  } // loop over all mortars
832 
833  // Apply DG mass matrix
834  if (massive) {
835  *operator_applied_to_vars /= get(det_inv_jacobian);
836  ::dg::apply_mass_matrix(operator_applied_to_vars, mesh);
837  }
838  }
839 
840  template <typename... FixedSourcesTags, typename ApplyBoundaryCondition,
841  typename... FluxesArgs, typename... SourcesArgs,
842  bool LocalLinearized = Linearized,
843  // This function adds nothing to the fixed sources if the operator
844  // is linearized, so it shouldn't be used in that case
847  const gsl::not_null<Variables<tmpl::list<FixedSourcesTags...>>*>
848  fixed_sources,
849  const Element<Dim>& element, const Mesh<Dim>& mesh,
850  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
851  inv_jacobian,
852  const Scalar<DataVector>& det_inv_jacobian,
853  const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
854  const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
855  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
856  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
857  const double penalty_parameter, const bool massive,
858  const ApplyBoundaryCondition& apply_boundary_condition,
859  const std::tuple<FluxesArgs...>& fluxes_args,
860  const std::tuple<SourcesArgs...>& sources_args,
862  fluxes_args_on_faces) noexcept {
863  // We just feed zero variables through the nonlinear operator to extract the
864  // constant contribution at external boundaries. Since the variables are
865  // zero the operator simplifies quite a lot. The simplification is probably
866  // not very important for performance because this function will only be
867  // called when solving a linear elliptic system and only once during
868  // initialization, but we specialize the operator for zero data nonetheless
869  // just so we can ignore internal boundaries. For internal boundaries we
870  // would unnecessarily have to copy mortar data around to emulate the
871  // communication step, so by just skipping internal boundaries we avoid
872  // that.
873  const size_t num_points = mesh.number_of_grid_points();
874  const Variables<tmpl::list<PrimalFields...>> zero_primal_vars{num_points,
875  0.};
876  Variables<tmpl::list<PrimalFluxes...>> primal_fluxes_buffer{num_points};
877  Variables<tmpl::list<AuxiliaryFields...>> unused_aux_vars_buffer{};
878  Variables<tmpl::list<AuxiliaryFluxes...>> unused_aux_fluxes_buffer{};
879  Variables<tmpl::list<FixedSourcesTags...>> operator_applied_to_zero_vars{
880  num_points};
881  // Set up data on mortars. We only need them at external boundaries.
882  ::dg::MortarMap<Dim, MortarData<size_t, tmpl::list<PrimalFields...>,
883  tmpl::list<PrimalFluxes...>>>
884  all_mortar_data{};
885  constexpr size_t temporal_id = std::numeric_limits<size_t>::max();
886  // Apply the operator to the zero variables, skipping internal boundaries
887  prepare_mortar_data<true>(
888  make_not_null(&unused_aux_vars_buffer),
889  make_not_null(&unused_aux_fluxes_buffer),
890  make_not_null(&primal_fluxes_buffer), make_not_null(&all_mortar_data),
891  zero_primal_vars, element, mesh, inv_jacobian, face_normals,
892  face_normal_magnitudes, all_mortar_meshes, all_mortar_sizes,
893  temporal_id, apply_boundary_condition, fluxes_args, sources_args,
894  fluxes_args_on_faces);
895  apply_operator<true>(make_not_null(&operator_applied_to_zero_vars),
896  make_not_null(&all_mortar_data), zero_primal_vars,
897  primal_fluxes_buffer, mesh, inv_jacobian,
898  det_inv_jacobian, face_normal_magnitudes,
899  all_mortar_meshes, all_mortar_sizes, penalty_parameter,
900  massive, temporal_id, sources_args);
901  // Impose the nonlinear (constant) boundary contribution as fixed sources on
902  // the RHS of the equations
903  *fixed_sources -= operator_applied_to_zero_vars;
904  }
905 };
906 
907 } // namespace detail
908 
909 /*!
910  * \brief Prepare data on mortars so they can be communicated to neighbors
911  *
912  * Call this function on all elements and communicate the mortar data, then call
913  * `elliptic::dg::apply_operator`.
914  */
915 template <typename System, bool Linearized, typename... Args>
916 void prepare_mortar_data(Args&&... args) noexcept {
917  detail::DgOperatorImpl<System, Linearized>::template prepare_mortar_data<
918  false>(std::forward<Args>(args)...);
919 }
920 
921 /*!
922  * \brief Apply the elliptic DG operator
923  *
924  * This function applies the elliptic DG operator on an element, assuming all
925  * data on mortars is already available. Use the
926  * `elliptic::dg::prepare_mortar_data` function to prepare mortar data on
927  * neighboring elements, then communicate the data and insert them on the
928  * "remote" side of the mortars before calling this function.
929  */
930 template <typename System, bool Linearized, typename... Args>
931 void apply_operator(Args&&... args) noexcept {
932  detail::DgOperatorImpl<System, Linearized>::template apply_operator<false>(
933  std::forward<Args>(args)...);
934 }
935 
936 /*!
937  * \brief For linear systems, impose inhomogeneous boundary conditions as
938  * contributions to the fixed sources (i.e. the RHS of the equations).
939  *
940  * This function exists because the DG operator must typically be linear, but
941  * even for linear elliptic equations we typically apply boundary conditions
942  * with a constant, and therefore nonlinear, contribution. Standard examples are
943  * inhomogeneous (i.e. non-zero) Dirichlet or Neumann boundary conditions. This
944  * nonlinear contribution can be added to the fixed sources, leaving only the
945  * linearized boundary conditions in the DG operator. For standard constant
946  * Dirichlet or Neumann boundary conditions the linearization is of course just
947  * zero.
948  *
949  * This function essentially feeds zero variables through the nonlinear operator
950  * and subtracts the result from the fixed sources: `b -= A(x=0)`.
951  */
952 template <typename System, typename... Args>
954  Args&&... args) noexcept {
955  detail::DgOperatorImpl<System, false>::
956  impose_inhomogeneous_boundary_conditions_on_source(
957  std::forward<Args>(args)...);
958 }
959 
960 } // namespace elliptic::dg
dg::SimpleBoundaryData::field_data
Variables< FieldTags > field_data
Data projected to the mortar mesh.
Definition: SimpleBoundaryData.hpp:37
Spectral::needs_projection
bool needs_projection(const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes) noexcept
Determine whether data needs to be projected between a child mesh and its parent mesh....
ElementId::external_boundary_id
static ElementId< VolumeDim > external_boundary_id() noexcept
Returns an ElementId meant for identifying data on external boundaries, which should never correspond...
EXPAND_PACK_LEFT_TO_RIGHT
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:601
utility
Divergence.hpp
index_to_slice_at
size_t index_to_slice_at(const Index< Dim > &extents, const Direction< Dim > &direction, const size_t offset=0) noexcept
Finds the index in the perpendicular dimension of an element boundary.
Definition: IndexToSliceAt.hpp:18
db::PrefixTag
Mark a struct as a prefix tag by inheriting from this.
Definition: Tag.hpp:103
elliptic::dg::Tags::PerpendicularNumPoints
Number of grid points perpendicular to an element face. Used to compute the penalty (see elliptic::dg...
Definition: DgOperator.hpp:154
elliptic::dg::prepare_mortar_data
void prepare_mortar_data(Args &&... args) noexcept
Prepare data on mortars so they can be communicated to neighbors.
Definition: DgOperator.hpp:916
LiftFlux.hpp
elliptic::dg::Tags::ElementSize
A measure of element size perpendicular to an element face. Used to compute the penalty (see elliptic...
Definition: DgOperator.hpp:160
dg::lift_flux
void lift_flux(const gsl::not_null< Variables< tmpl::list< BoundaryCorrectionTags... >> * > boundary_correction_terms, const size_t extent_perpendicular_to_boundary, Scalar< DataVector > magnitude_of_face_normal) noexcept
Lifts the flux contribution from an interface to the volume.
Definition: LiftFlux.hpp:42
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::NormalDotFlux
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:96
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
divergence
auto divergence(const Variables< FluxTags > &F, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept -> Variables< db::wrap_tags_in< Tags::div, FluxTags >>
Compute the (Euclidean) divergence of fluxes.
Spectral.hpp
Direction< Dim >
Element
Definition: Element.hpp:29
ElementId.hpp
elliptic::dg::MortarData
Boundary data on both sides of a mortar.
Definition: DgOperator.hpp:214
elliptic::dg::Tags::MortarData
Holds elliptic::dg::MortarData, i.e. boundary data on both sides of a mortar.
Definition: DgOperator.hpp:223
db::data_on_slice
Variables< tmpl::list< TagsToSlice... > > data_on_slice(const db::DataBox< TagsList > &box, const Index< VolumeDim > &element_extents, const size_t sliced_dim, const size_t fixed_index, tmpl::list< TagsToSlice... >) noexcept
Slices volume Tensors from a DataBox into a Variables
Definition: DataOnSlice.hpp:33
elliptic::dg::Tags::NormalDotFluxForJump
The quantity where is the system flux for primal variables and auxiliary variables ,...
Definition: DgOperator.hpp:169
make_array
constexpr std::array< T, Size > make_array(Args &&... args) noexcept(noexcept(MakeArray_detail::MakeArray< Size==0 >::template apply< T >(std::make_index_sequence<(Size==0 ? Size :Size - 1)>{}, std::forward< Args >(args)...)))
Create a std::array<T, Size>{{T(args...), T(args...), ...}}
Definition: MakeArray.hpp:68
dg::SimpleBoundaryData
Distinguishes between field data, which can be projected to a mortar, and extra data,...
Definition: SimpleBoundaryData.hpp:32
cstddef
Assert.hpp
add_slice_to_data
void add_slice_to_data(const gsl::not_null< Variables< tmpl::list< VolumeTags... >> * > volume_vars, const Variables< tmpl::list< SliceTags... >> &vars_on_slice, const Index< VolumeDim > &extents, const size_t sliced_dim, const size_t fixed_index) noexcept
Adds data on a codimension 1 slice to a volume quantity. The slice has a constant logical coordinate ...
Definition: SliceVariables.hpp:82
elliptic::get_sources_computer
tmpl::conditional_t< Linearized, typename detail::sources_computer_linearized< System >::type, typename System::sources_computer > get_sources_computer
The System::sources_computer or the System::sources_computer_linearized, depending on the Linearized ...
Definition: GetSourcesComputer.hpp:30
std::array
elliptic::dg::BoundaryData
::dg::SimpleBoundaryData< tmpl::append< db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFields >, db::wrap_tags_in<::Tags::NormalDotFlux, PrimalFluxes >, db::wrap_tags_in< Tags::NormalDotFluxForJump, PrimalFields >, tmpl::list< Tags::ElementSize > >, tmpl::list< Tags::PerpendicularNumPoints > > BoundaryData
Data that is projected to mortars and communicated across element boundaries.
Definition: DgOperator.hpp:183
elliptic::dg::apply_operator
void apply_operator(Args &&... args) noexcept
Apply the elliptic DG operator.
Definition: DgOperator.hpp:931
elliptic::dg
Functionality related to discontinuous Galerkin discretizations of elliptic equations.
Definition: ApplyOperator.hpp:46
Elasticity::primal_fluxes
void primal_fluxes(gsl::not_null< tnsr::II< DataVector, Dim > * > flux_for_displacement, const tnsr::ii< DataVector, Dim > &strain, const ConstitutiveRelations::ConstitutiveRelation< Dim > &constitutive_relation, const tnsr::I< DataVector, Dim > &coordinates) noexcept
Compute the fluxes for the Elasticity equation.
std::min
T min(T... args)
dg::project_from_mortar
Variables< Tags > project_from_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:85
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Element::external_boundaries
const std::unordered_set< Direction< VolumeDim > > & external_boundaries() const noexcept
The directions of the faces of the Element that are external boundaries.
Definition: Element.hpp:51
Variables.hpp
Element.hpp
dg::SimpleMortarData
Storage of boundary data on two sides of a mortar.
Definition: SimpleMortarData.hpp:22
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
DirectionMap
Definition: DirectionMap.hpp:15
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Elasticity::auxiliary_fluxes
void auxiliary_fluxes(gsl::not_null< tnsr::Ijj< DataVector, Dim > * > flux_for_strain, const tnsr::I< DataVector, Dim > &displacement) noexcept
Compute the fluxes for the auxiliary (strain) field in the first-order formulation of the Elasticity...
Tensor.hpp
elliptic::dg::zero_boundary_data_on_mortar
BoundaryData< PrimalMortarFields, PrimalMortarFluxes > zero_boundary_data_on_mortar(const Direction< Dim > &direction, const Mesh< Dim > &mesh, const Scalar< DataVector > &face_normal_magnitude, const Mesh< Dim - 1 > &mortar_mesh, const ::dg::MortarSize< Dim - 1 > &mortar_size) noexcept
Construct elliptic::dg::BoundaryData assuming the variable data on the element is zero,...
Definition: DgOperator.hpp:189
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
FixedHashMap::contains
bool contains(const key_type &key) const noexcept
Check if key is in the map.
Definition: FixedHashMap.hpp:177
elliptic::dg::penalty
DataVector penalty(const DataVector &element_size, size_t num_points, double penalty_parameter) noexcept
The penalty factor in internal penalty fluxes.
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
Direction.hpp
std::numeric_limits::max
T max(T... args)
unordered_map
elliptic::apply_boundary_condition
void apply_boundary_condition(const elliptic::BoundaryConditions::BoundaryCondition< Dim, Registrars > &boundary_condition, const db::DataBox< DbTagsList > &box, const MapKeys &map_keys_to_direction, const gsl::not_null< FieldsAndFluxes * >... fields_and_fluxes) noexcept
Apply the boundary_condition to the fields_and_fluxes with arguments from interface tags in the DataB...
Definition: ApplyBoundaryCondition.hpp:41
dg::mortar_size
MortarSize< Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, size_t dimension, const OrientationMap< Dim > &orientation) noexcept
TMPL.hpp
Element::neighbors
const Neighbors_t & neighbors() const noexcept
Information about the neighboring Elements.
Definition: Element.hpp:66
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
void impose_inhomogeneous_boundary_conditions_on_source(Args &&... args) noexcept
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: DgOperator.hpp:953
string