DgOperator.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/range/join.hpp>
7 #include <cstddef>
8 #include <string>
9 #include <tuple>
10 #include <unordered_map>
11 #include <utility>
12 
13 #include "DataStructures/DataVector.hpp"
14 #include "DataStructures/SliceVariables.hpp"
18 #include "Domain/Structure/DirectionMap.hpp"
21 #include "Domain/Structure/IndexToSliceAt.hpp"
22 #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
23 #include "Elliptic/Systems/GetSourcesComputer.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 std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
278  internal_face_normals,
279  const std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
280  external_face_normals,
282  internal_face_normal_magnitudes,
284  external_face_normal_magnitudes,
285  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
286  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
287  const TemporalId& temporal_id,
288  const ApplyBoundaryCondition& apply_boundary_condition,
289  const std::tuple<FluxesArgs...>& fluxes_args,
290  const std::tuple<SourcesArgs...>& sources_args,
292  fluxes_args_on_internal_faces,
294  fluxes_args_on_external_faces,
295  const DirectionsPredicate& directions_predicate =
296  AllDirections{}) noexcept {
297  static_assert(
298  sizeof...(PrimalVars) == sizeof...(PrimalFields) and
299  sizeof...(AuxiliaryVars) == sizeof...(AuxiliaryFields) and
300  sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
301  sizeof...(AuxiliaryFluxesVars) == sizeof...(AuxiliaryFluxes),
302  "The number of variables must match the number of system fields.");
303  static_assert(
304  (std::is_same_v<typename PrimalVars::type,
305  typename PrimalFields::type> and
306  ...) and
307  (std::is_same_v<typename AuxiliaryVars::type,
308  typename AuxiliaryFields::type> and
309  ...) and
310  (std::is_same_v<typename PrimalFluxesVars::type,
311  typename PrimalFluxes::type> and
312  ...) and
313  (std::is_same_v<typename AuxiliaryFluxesVars::type,
314  typename AuxiliaryFluxes::type> and
315  ...),
316  "The variables must have the same tensor types as the system fields.");
317 #ifdef SPECTRE_DEBUG
318  for (size_t d = 0; d < Dim; ++d) {
319  ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
320  mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
321  "The elliptic DG operator is currently only implemented for "
322  "Legendre-Gauss-Lobatto grids. Found basis '"
323  << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
324  << "' in dimension " << d << ".");
325  }
326 #endif // SPECTRE_DEBUG
327  const size_t num_points = mesh.number_of_grid_points();
328 
329  // This function and the one below allocate various Variables to compute
330  // intermediate quantities. It could be a performance optimization to reduce
331  // the number of these allocations and/or move some of the memory buffers
332  // into the DataBox to keep them around permanently. The latter should be
333  // informed by profiling.
334 
335  // Compute the auxiliary variables, and from those the primal fluxes. The
336  // auxiliary variables are the variables `v` in the auxiliary equations
337  // `-div(F_v(u)) + S_v + v = 0` where `F_v` and `S_v` are the auxiliary
338  // fluxes and sources. From the auxiliary variables we compute the primal
339  // fluxes `F_u(v)` for the primal equation `-div(F_u(v)) + S_u = f(x)`. Note
340  // that before taking the second derivative, boundary corrections from the
341  // first derivative have to be added to `F_u(v)`. Therefore the second
342  // derivative is taken after a communication break in the `apply_operator`
343  // function below.
344  if constexpr (DataIsZero) {
345  (void)auxiliary_vars;
346  (void)auxiliary_fluxes;
347  primal_fluxes->initialize(num_points, 0.);
348  } else {
349  // Compute the auxiliary fluxes
350  auxiliary_fluxes->initialize(num_points);
351  std::apply(
352  [&auxiliary_fluxes,
353  &primal_vars](const auto&... expanded_fluxes_args) noexcept {
354  FluxesComputer::apply(
355  make_not_null(&get<AuxiliaryFluxesVars>(*auxiliary_fluxes))...,
356  expanded_fluxes_args..., get<PrimalVars>(primal_vars)...);
357  },
358  fluxes_args);
359  divergence(auxiliary_vars, *auxiliary_fluxes, mesh, inv_jacobian);
360  // Subtract the sources. Invert signs because the sources computers _add_
361  // the sources, i.e. they compute the `+ S_v` term defined above.
362  *auxiliary_vars *= -1.;
363  std::apply(
364  [&auxiliary_vars,
365  &primal_vars](const auto&... expanded_sources_args) noexcept {
366  SourcesComputer::apply(
367  make_not_null(&get<AuxiliaryVars>(*auxiliary_vars))...,
368  expanded_sources_args..., get<PrimalVars>(primal_vars)...);
369  },
370  sources_args);
371  *auxiliary_vars *= -1.;
372  // Compute the primal fluxes
373  primal_fluxes->initialize(num_points);
374  std::apply(
375  [&primal_fluxes,
376  &auxiliary_vars](const auto&... expanded_fluxes_args) noexcept {
377  FluxesComputer::apply(
378  make_not_null(&get<PrimalFluxesVars>(*primal_fluxes))...,
379  expanded_fluxes_args...,
380  get<AuxiliaryVars>(*auxiliary_vars)...);
381  },
382  fluxes_args);
383  }
384 
385  // Populate the mortar data on this element's side of the boundary so it's
386  // ready to be sent to neighbors.
387  for (const auto& direction : [&element]() noexcept -> decltype(auto) {
388  if constexpr (DataIsZero) {
389  // Skipping internal boundaries for zero data because they won't
390  // contribute boundary corrections anyway.
391  return element.external_boundaries();
392  } else {
393  return boost::join(element.internal_boundaries(),
394  element.external_boundaries());
395  };
396  }()) {
397  if (not directions_predicate(direction)) {
398  continue;
399  }
400  const bool is_internal = element.neighbors().contains(direction);
401  const auto face_mesh = mesh.slice_away(direction.dimension());
402  const size_t face_num_points = face_mesh.number_of_grid_points();
403  const auto& face_normal = is_internal
404  ? internal_face_normals.at(direction)
405  : external_face_normals.at(direction);
406  const auto& face_normal_magnitude =
407  is_internal ? internal_face_normal_magnitudes.at(direction)
408  : external_face_normal_magnitudes.at(direction);
409  const auto& fluxes_args_on_face =
410  is_internal ? fluxes_args_on_internal_faces.at(direction)
411  : fluxes_args_on_external_faces.at(direction);
412  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
413  Variables<tmpl::list<PrimalFluxesVars..., AuxiliaryFluxesVars...>>
414  fluxes_on_face{face_num_points};
415  Variables<tmpl::list<::Tags::NormalDotFlux<AuxiliaryFields>...>>
416  n_dot_aux_fluxes{face_num_points};
417  BoundaryData<tmpl::list<PrimalMortarVars...>,
418  tmpl::list<PrimalMortarFluxes...>>
419  boundary_data{face_num_points};
420  if constexpr (DataIsZero) {
421  // Just setting all boundary field data to zero. Variable-independent
422  // data such as the element size will be set below.
423  boundary_data.field_data.initialize(face_num_points, 0.);
424  } else {
425  // Compute F_u(n.F_v)
426  fluxes_on_face.assign_subset(
427  data_on_slice(*auxiliary_fluxes, mesh.extents(),
428  direction.dimension(), slice_index));
429  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
431  &get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
432  face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
433  std::apply(
434  [&boundary_data, &n_dot_aux_fluxes](
435  const auto&... expanded_fluxes_args_on_face) noexcept {
436  FluxesComputer::apply(
438  boundary_data.field_data))...,
439  expanded_fluxes_args_on_face...,
441  n_dot_aux_fluxes)...);
442  },
443  fluxes_args_on_face);
444  // Compute n.F_u
445  //
446  // For the internal penalty flux we can already slice the n.F_u to the
447  // faces at this point, before the boundary corrections have been added
448  // to the auxiliary variables. The reason is essentially that the
449  // internal penalty flux uses average(grad(u)) in place of average(v),
450  // i.e. the raw primal field derivatives without boundary corrections.
451  fluxes_on_face.assign_subset(
452  data_on_slice(*primal_fluxes, mesh.extents(), direction.dimension(),
453  slice_index));
454  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
456  boundary_data.field_data)),
457  face_normal, get<PrimalFluxesVars>(fluxes_on_face)));
458  // Compute n.F_u(n.F_v) for jump term, re-using the memory buffer from
459  // above
460  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
461  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
462  boundary_data.field_data)),
463  face_normal,
465  boundary_data.field_data)));
466  }
467 
468  // Collect the remaining data that's needed on both sides of the boundary
469  // These are actually constant throughout the solve, so a performance
470  // optimization could be to store them in the DataBox.
471  get<Tags::PerpendicularNumPoints>(boundary_data.extra_data) =
472  mesh.extents(direction.dimension());
473  get(get<Tags::ElementSize>(boundary_data.field_data)) =
474  2. / get(face_normal_magnitude);
475 
476  if (is_internal) {
477  if constexpr (not DataIsZero) {
478  // Project boundary data on internal faces to mortars
479  for (const auto& neighbor_id : element.neighbors().at(direction)) {
480  const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
481  const auto& mortar_mesh = all_mortar_meshes.at(mortar_id);
482  const auto& mortar_size = all_mortar_sizes.at(mortar_id);
483  auto projected_boundary_data =
485  ? boundary_data.project_to_mortar(face_mesh, mortar_mesh,
486  mortar_size)
487  : std::move(boundary_data);
488  (*all_mortar_data)[mortar_id].local_insert(
489  temporal_id, std::move(projected_boundary_data));
490  }
491  }
492  } else {
493  // No need to do projections on external boundaries
494  const ::dg::MortarId<Dim> mortar_id{
496  (*all_mortar_data)[mortar_id].local_insert(temporal_id, boundary_data);
497 
498  // -------------------------
499  // Apply boundary conditions
500  // -------------------------
501  //
502  // To apply boundary conditions we fill the boundary data with
503  // "exterior" or "ghost" data and set it as remote mortar data, so
504  // external boundaries behave just like internal boundaries when
505  // applying boundary corrections.
506  //
507  // The `apply_boundary_conditions` invocable is expected to impose
508  // boundary conditions by modifying the fields and fluxes that are
509  // passed by reference. Dirichlet-type boundary conditions are imposed
510  // by modifying the fields, and Neumann-type boundary conditions are
511  // imposed by modifying the interior n.F_u. Note that all data passed to
512  // the boundary conditions is taken from the "interior" side of the
513  // boundary, i.e. with a normal vector that points _out_ of the
514  // computational domain.
515  auto dirichlet_vars = data_on_slice(primal_vars, mesh.extents(),
516  direction.dimension(), slice_index);
518  direction, make_not_null(&get<PrimalVars>(dirichlet_vars))...,
520  boundary_data.field_data))...);
521 
522  // The n.F_u (Neumann-type conditions) are done, but we have to compute
523  // the n.F_v (Dirichlet-type conditions) from the Dirichlet fields. We
524  // re-use the memory buffer from above.
525  std::apply(
526  [&fluxes_on_face, &dirichlet_vars](
527  const auto&... expanded_fluxes_args_on_face) noexcept {
528  FluxesComputer::apply(
529  make_not_null(&get<AuxiliaryFluxesVars>(fluxes_on_face))...,
530  expanded_fluxes_args_on_face...,
531  get<PrimalVars>(dirichlet_vars)...);
532  },
533  fluxes_args_on_face);
534  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
536  &get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
537  face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
538  std::apply(
539  [&boundary_data, &n_dot_aux_fluxes](
540  const auto&... expanded_fluxes_args_on_face) noexcept {
541  FluxesComputer::apply(
543  boundary_data.field_data))...,
544  expanded_fluxes_args_on_face...,
546  n_dot_aux_fluxes)...);
547  },
548  fluxes_args_on_face);
549 
550  // Invert the sign of the fluxes to account for the inverted normal on
551  // exterior faces. Also multiply by 2 and add the interior fluxes to
552  // impose the boundary conditions on the _average_ instead of just
553  // setting the fields on the exterior.
554  const auto invert_sign_and_impose_on_average =
555  [](const auto exterior_n_dot_flux,
556  const auto& interior_n_dot_flux) {
557  for (size_t i = 0; i < interior_n_dot_flux.size(); ++i) {
558  (*exterior_n_dot_flux)[i] *= -2.;
559  (*exterior_n_dot_flux)[i] += interior_n_dot_flux[i];
560  }
561  };
562  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
564  boundary_data.field_data)),
566  all_mortar_data->at(mortar_id)
567  .local_data(temporal_id)
568  .field_data)));
569  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign_and_impose_on_average(
571  boundary_data.field_data)),
573  all_mortar_data->at(mortar_id)
574  .local_data(temporal_id)
575  .field_data)));
576 
577  // Compute n.F_u(n.F_v) for jump term
578  EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
579  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
580  boundary_data.field_data)),
581  face_normal,
583  boundary_data.field_data)));
584  const auto invert_sign = [](const auto exterior_n_dot_flux) {
585  for (size_t i = 0; i < exterior_n_dot_flux->size(); ++i) {
586  (*exterior_n_dot_flux)[i] *= -1.;
587  }
588  };
589  EXPAND_PACK_LEFT_TO_RIGHT(invert_sign(
590  make_not_null(&get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
591  boundary_data.field_data))));
592 
593  // Store the exterior boundary data on the mortar
594  all_mortar_data->at(mortar_id).remote_insert(temporal_id,
595  std::move(boundary_data));
596  } // if (is_internal)
597  } // loop directions
598  }
599 
600  // --- This is essentially a break to communicate the mortar data ---
601 
602  template <bool DataIsZero, typename... OperatorTags, typename... PrimalVars,
603  typename... PrimalFluxesVars, typename... PrimalMortarVars,
604  typename... PrimalMortarFluxes, typename TemporalId,
605  typename... FluxesArgs, typename... SourcesArgs,
606  typename DirectionsPredicate = AllDirections>
607  static void apply_operator(
608  const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
609  operator_applied_to_vars,
611  Dim, MortarData<TemporalId, tmpl::list<PrimalMortarVars...>,
612  tmpl::list<PrimalMortarFluxes...>>>*>
613  all_mortar_data,
614  const Variables<tmpl::list<PrimalVars...>>& primal_vars,
615  // Taking the primal fluxes computed in the `prepare_mortar_data` function
616  // by const-ref here because other code might use them and so we don't
617  // want to modify them by adding boundary corrections. E.g. linearized
618  // sources use the nonlinear fields and fluxes as background fields.
619  const Variables<tmpl::list<PrimalFluxesVars...>>& primal_fluxes,
620  const Mesh<Dim>& mesh,
621  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
622  inv_jacobian,
624  internal_face_normal_magnitudes,
626  external_face_normal_magnitudes,
627  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
628  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
629  const double penalty_parameter, const TemporalId& temporal_id,
630  const std::tuple<SourcesArgs...>& sources_args,
631  const DirectionsPredicate& directions_predicate =
632  AllDirections{}) noexcept {
633  static_assert(
634  sizeof...(PrimalVars) == sizeof...(PrimalFields) and
635  sizeof...(PrimalFluxesVars) == sizeof...(PrimalFluxes) and
636  sizeof...(PrimalMortarVars) == sizeof...(PrimalFields) and
637  sizeof...(PrimalMortarFluxes) == sizeof...(PrimalFluxes) and
638  sizeof...(OperatorTags) == sizeof...(PrimalFields),
639  "The number of variables must match the number of system fields.");
640  static_assert(
641  (std::is_same_v<typename PrimalVars::type,
642  typename PrimalFields::type> and
643  ...) and
644  (std::is_same_v<typename PrimalFluxesVars::type,
645  typename PrimalFluxes::type> and
646  ...) and
647  (std::is_same_v<typename PrimalMortarVars::type,
648  typename PrimalFields::type> and
649  ...) and
650  (std::is_same_v<typename PrimalMortarFluxes::type,
651  typename PrimalFluxes::type> and
652  ...) and
653  (std::is_same_v<typename OperatorTags::type,
654  typename PrimalFields::type> and
655  ...),
656  "The variables must have the same tensor types as the system fields.");
657 #ifdef SPECTRE_DEBUG
658  for (size_t d = 0; d < Dim; ++d) {
659  ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
660  mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
661  "The elliptic DG operator is currently only implemented for "
662  "Legendre-Gauss-Lobatto grids. Found basis '"
663  << mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
664  << "' in dimension " << d << ".");
665  }
666 #endif // SPECTRE_DEBUG
667 
668  // This function and the one above allocate various Variables to compute
669  // intermediate quantities. It could be a performance optimization to reduce
670  // the number of these allocations and/or move some of the memory buffers
671  // into the DataBox to keep them around permanently. The latter should be
672  // informed by profiling.
673 
674  // Add boundary corrections to the auxiliary variables _before_ computing
675  // the second derivative. This is called the "flux" formulation. It is
676  // equivalent to discretizing the system in first-order form, i.e. treating
677  // the primal and auxiliary variables on the same footing, and then taking a
678  // Schur complement of the operator. The Schur complement is possible
679  // because the auxiliary equations are essentially the definition of the
680  // auxiliary variables and can therefore always be solved for the auxiliary
681  // variables by just inverting the mass matrix. This Schur-complement
682  // formulation avoids inflating the DG operator with the DOFs of the
683  // auxiliary variables. In this form it is very similar to the "primal"
684  // formulation where we get rid of the auxiliary variables through a DG
685  // theorem and thus add the auxiliary boundary corrections _after_ computing
686  // the second derivative. This involves a slightly different lifting
687  // operation with differentiation matrices, which we avoid to implement for
688  // now by using the flux-formulation.
689  auto primal_fluxes_corrected = primal_fluxes;
690  for (const auto& [mortar_id, mortar_data] : *all_mortar_data) {
691  const auto& [direction, neighbor_id] = mortar_id;
692  const bool is_internal =
693  (neighbor_id != ElementId<Dim>::external_boundary_id());
694  if constexpr (DataIsZero) {
695  if (is_internal) {
696  continue;
697  }
698  }
699  if (not directions_predicate(direction)) {
700  continue;
701  }
702  const auto face_mesh = mesh.slice_away(direction.dimension());
703  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
704  const auto& local_data = mortar_data.local_data(temporal_id);
705  const auto& remote_data = mortar_data.remote_data(temporal_id);
706  const auto& face_normal_magnitude =
707  is_internal ? internal_face_normal_magnitudes.at(direction)
708  : external_face_normal_magnitudes.at(direction);
709  const auto& mortar_mesh =
710  is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
711  const auto& mortar_size =
712  is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
713 
714  // This is the _strong_ auxiliary boundary correction avg(n.F_v) - n.F_v
715  auto auxiliary_boundary_corrections_on_mortar =
716  Variables<tmpl::list<PrimalMortarFluxes...>>(
717  local_data.field_data.template extract_subset<
718  tmpl::list<::Tags::NormalDotFlux<PrimalMortarFluxes>...>>());
719  const auto add_remote_contribution = [](auto& lhs,
720  const auto& rhs) noexcept {
721  for (size_t i = 0; i < lhs.size(); ++i) {
722  lhs[i] += rhs[i];
723  }
724  };
725  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_contribution(
726  get<PrimalMortarFluxes>(auxiliary_boundary_corrections_on_mortar),
728  remote_data.field_data)));
729  auxiliary_boundary_corrections_on_mortar *= -0.5;
730 
731  // Project from the mortar back down to the face if needed
732  auto auxiliary_boundary_corrections =
735  auxiliary_boundary_corrections_on_mortar, face_mesh,
737  : std::move(auxiliary_boundary_corrections_on_mortar);
738 
739  // Lift the boundary correction to the volume, but still only provide the
740  // data only on the face because it is zero everywhere else. This is the
741  // "massless" lifting operation, i.e. it involves an inverse mass matrix.
742  // The mass matrix is diagonally approximated ("mass lumping") so it
743  // reduces to a division by quadrature weights.
744  ::dg::lift_flux(make_not_null(&auxiliary_boundary_corrections),
745  mesh.extents(direction.dimension()),
746  face_normal_magnitude);
747  // The `dg::lift_flux` function contains an extra minus sign
748  auxiliary_boundary_corrections *= -1.;
749 
750  // Add the boundary corrections to the auxiliary variables
751  add_slice_to_data(make_not_null(&primal_fluxes_corrected),
752  auxiliary_boundary_corrections, mesh.extents(),
753  direction.dimension(), slice_index);
754  } // apply auxiliary boundary corrections on all mortars
755 
756  // Compute the primal equation, i.e. the actual DG operator, by taking the
757  // second derivative: -div(F_u(v)) + S_u = f(x)
758  divergence(operator_applied_to_vars, primal_fluxes_corrected, mesh,
759  inv_jacobian);
760  // This is the sign flip that makes the operator _minus_ the Laplacian for a
761  // Poisson system
762  *operator_applied_to_vars *= -1.;
763  // Using the non-boundary-corrected primal fluxes to compute sources here
764  std::apply(
765  [&operator_applied_to_vars, &primal_vars,
766  &primal_fluxes](const auto&... expanded_sources_args) noexcept {
767  SourcesComputer::apply(
768  make_not_null(&get<OperatorTags>(*operator_applied_to_vars))...,
769  expanded_sources_args..., get<PrimalVars>(primal_vars)...,
770  get<PrimalFluxesVars>(primal_fluxes)...);
771  },
772  sources_args);
773 
774  // Add boundary corrections to primal equation
775  for (auto& [mortar_id, mortar_data] : *all_mortar_data) {
776  const auto& [direction, neighbor_id] = mortar_id;
777  const bool is_internal =
778  (neighbor_id != ElementId<Dim>::external_boundary_id());
779  if constexpr (DataIsZero) {
780  if (is_internal) {
781  continue;
782  }
783  }
784  const auto face_mesh = mesh.slice_away(direction.dimension());
785  const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
786  const auto [local_data, remote_data] = mortar_data.extract();
787  const auto& face_normal_magnitude =
788  is_internal ? internal_face_normal_magnitudes.at(direction)
789  : external_face_normal_magnitudes.at(direction);
790  const auto& mortar_mesh =
791  is_internal ? all_mortar_meshes.at(mortar_id) : face_mesh;
792  const auto& mortar_size =
793  is_internal ? all_mortar_sizes.at(mortar_id) : full_mortar_size;
794 
795  // This is the _strong_ primal boundary correction avg(n.F_u) - penalty *
796  // jump(n.F_v) - n.F_u. Note that the "internal penalty" numerical flux
797  // (as opposed to the LLF flux) uses the raw field derivatives without
798  // boundary corrections in the average, which is why we can communicate
799  // the data so early together with the auxiliary boundary data. In this
800  // case the penalty needs to include a factor N_points^2 / h (see the
801  // `penalty` function).
802  const auto penalty = elliptic::dg::penalty(
803  min(get(get<Tags::ElementSize>(local_data.field_data)),
804  get(get<Tags::ElementSize>(remote_data.field_data))),
805  std::max(get<Tags::PerpendicularNumPoints>(local_data.extra_data),
806  get<Tags::PerpendicularNumPoints>(remote_data.extra_data)),
807  penalty_parameter);
808  // Start with the penalty term
809  auto primal_boundary_corrections_on_mortar =
810  Variables<tmpl::list<PrimalMortarVars...>>(
811  local_data.field_data.template extract_subset<tmpl::list<
812  Tags::NormalDotFluxForJump<PrimalMortarVars>...>>());
813  const auto add_remote_jump_contribution = [](auto& lhs,
814  const auto& rhs) noexcept {
815  for (size_t i = 0; i < lhs.size(); ++i) {
816  lhs[i] -= rhs[i];
817  }
818  };
819  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_jump_contribution(
820  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
821  get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
822  remote_data.field_data)));
823  primal_boundary_corrections_on_mortar *= penalty;
824  const auto add_remote_avg_contribution = [](auto& lhs,
825  const auto& rhs) noexcept {
826  for (size_t i = 0; i < lhs.size(); ++i) {
827  lhs[i] += 0.5 * rhs[i];
828  }
829  };
830  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
831  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
833  local_data.field_data)));
834  EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
835  get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
837  remote_data.field_data)));
838  primal_boundary_corrections_on_mortar *= -1.;
839 
840  // Project from the mortar back down to the face if needed, lift and add
841  // to operator. See auxiliary boundary corrections above for details.
842  auto primal_boundary_corrections =
844  ? ::dg::project_from_mortar(primal_boundary_corrections_on_mortar,
845  face_mesh, mortar_mesh, mortar_size)
846  : std::move(primal_boundary_corrections_on_mortar);
847  ::dg::lift_flux(make_not_null(&primal_boundary_corrections),
848  mesh.extents(direction.dimension()),
849  face_normal_magnitude);
850  add_slice_to_data(operator_applied_to_vars, primal_boundary_corrections,
851  mesh.extents(), direction.dimension(), slice_index);
852  } // loop over all mortars
853  }
854 
855  template <typename... FixedSourcesTags, typename ApplyBoundaryCondition,
856  typename... FluxesArgs, typename... SourcesArgs,
857  bool LocalLinearized = Linearized,
858  // This function adds nothing to the fixed sources if the operator
859  // is linearized, so it shouldn't be used in that case
862  const gsl::not_null<Variables<tmpl::list<FixedSourcesTags...>>*>
863  fixed_sources,
864  const Element<Dim>& element, const Mesh<Dim>& mesh,
865  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
866  inv_jacobian,
867  const std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
868  external_face_normals,
870  external_face_normal_magnitudes,
871  const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
872  const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
873  const double penalty_parameter,
874  const ApplyBoundaryCondition& apply_boundary_condition,
875  const std::tuple<FluxesArgs...>& fluxes_args,
876  const std::tuple<SourcesArgs...>& sources_args,
878  fluxes_args_on_internal_faces,
880  fluxes_args_on_external_faces) noexcept {
881  // We just feed zero variables through the nonlinear operator to extract the
882  // constant contribution at external boundaries. Since the variables are
883  // zero the operator simplifies quite a lot. The simplification is probably
884  // not very important for performance because this function will only be
885  // called when solving a linear elliptic system and only once during
886  // initialization, but we specialize the operator for zero data nonetheless
887  // just so we can ignore internal boundaries. For internal boundaries we
888  // would unnecessarily have to copy mortar data around to emulate the
889  // communication step, so by just skipping internal boundaries we avoid
890  // that.
891  const size_t num_points = mesh.number_of_grid_points();
892  const Variables<tmpl::list<PrimalFields...>> zero_primal_vars{num_points,
893  0.};
894  Variables<tmpl::list<PrimalFluxes...>> primal_fluxes_buffer{num_points};
895  Variables<tmpl::list<AuxiliaryFields...>> unused_aux_vars_buffer{};
896  Variables<tmpl::list<AuxiliaryFluxes...>> unused_aux_fluxes_buffer{};
897  Variables<tmpl::list<FixedSourcesTags...>> operator_applied_to_zero_vars{
898  num_points};
899  // Set up data on mortars. We only need them at external boundaries.
900  ::dg::MortarMap<Dim, MortarData<size_t, tmpl::list<PrimalFields...>,
901  tmpl::list<PrimalFluxes...>>>
902  all_mortar_data{};
903  constexpr size_t temporal_id = std::numeric_limits<size_t>::max();
904  // Apply the operator to the zero variables, skipping internal boundaries
905  prepare_mortar_data<true>(
906  make_not_null(&unused_aux_vars_buffer),
907  make_not_null(&unused_aux_fluxes_buffer),
908  make_not_null(&primal_fluxes_buffer),
909  make_not_null(&all_mortar_data), zero_primal_vars, element, mesh,
910  inv_jacobian, {}, external_face_normals, {},
911  external_face_normal_magnitudes, all_mortar_meshes, all_mortar_sizes,
912  temporal_id, apply_boundary_condition, fluxes_args, sources_args,
913  fluxes_args_on_internal_faces, fluxes_args_on_external_faces);
914  apply_operator<true>(make_not_null(&operator_applied_to_zero_vars),
915  make_not_null(&all_mortar_data), zero_primal_vars,
916  primal_fluxes_buffer, mesh, inv_jacobian, {},
917  external_face_normal_magnitudes, all_mortar_meshes,
918  all_mortar_sizes, penalty_parameter, temporal_id,
919  sources_args);
920  // Impose the nonlinear (constant) boundary contribution as fixed sources on
921  // the RHS of the equations
922  *fixed_sources -= operator_applied_to_zero_vars;
923  }
924 };
925 
926 } // namespace detail
927 
928 /*!
929  * \brief Prepare data on mortars so they can be communicated to neighbors
930  *
931  * Call this function on all elements and communicate the mortar data, then call
932  * `elliptic::dg::apply_operator`.
933  */
934 template <typename System, bool Linearized, typename... Args>
935 void prepare_mortar_data(Args&&... args) noexcept {
936  detail::DgOperatorImpl<System, Linearized>::template prepare_mortar_data<
937  false>(std::forward<Args>(args)...);
938 }
939 
940 /*!
941  * \brief Apply the elliptic DG operator
942  *
943  * This function applies the elliptic DG operator on an element, assuming all
944  * data on mortars is already available. Use the
945  * `elliptic::dg::prepare_mortar_data` function to prepare mortar data on
946  * neighboring elements, then communicate the data and insert them on the
947  * "remote" side of the mortars before calling this function.
948  */
949 template <typename System, bool Linearized, typename... Args>
950 void apply_operator(Args&&... args) noexcept {
951  detail::DgOperatorImpl<System, Linearized>::template apply_operator<false>(
952  std::forward<Args>(args)...);
953 }
954 
955 /*!
956  * \brief For linear systems, impose inhomogeneous boundary conditions as
957  * contributions to the fixed sources (i.e. the RHS of the equations).
958  *
959  * This function exists because the DG operator must typically be linear, but
960  * even for linear elliptic equations we typically apply boundary conditions
961  * with a constant, and therefore nonlinear, contribution. Standard examples are
962  * inhomogeneous (i.e. non-zero) Dirichlet or Neumann boundary conditions. This
963  * nonlinear contribution can be added to the fixed sources, leaving only the
964  * linearized boundary conditions in the DG operator. For standard constant
965  * Dirichlet or Neumann boundary conditions the linearization is of course just
966  * zero.
967  *
968  * This function essentially feeds zero variables through the nonlinear operator
969  * and subtracts the result from the fixed sources: `b -= A(x=0)`.
970  */
971 template <typename System, typename... Args>
973  Args&&... args) noexcept {
974  detail::DgOperatorImpl<System, false>::
975  impose_inhomogeneous_boundary_conditions_on_source(
976  std::forward<Args>(args)...);
977 }
978 
979 } // namespace elliptic::dg
dg::SimpleBoundaryData::field_data
Variables< FieldTags > field_data
Data projected to the mortar mesh.
Definition: SimpleBoundaryData.hpp:37
elliptic::apply_boundary_condition
void apply_boundary_condition(const elliptic::BoundaryConditions::BoundaryCondition< Dim, Registrars > &boundary_condition, const db::DataBox< DbTagsList > &box, const Direction< Dim > &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:36
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:565
utility
Divergence.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
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:935
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::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::internal_boundaries
const std::unordered_set< Direction< VolumeDim > > & internal_boundaries() const noexcept
The directions of the faces of the Element that are internal boundaries.
Definition: Element.hpp:57
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:950
elliptic::dg
Functionality related to discontinuous Galerkin discretizations of elliptic equations.
Definition: ApplyOperator.hpp:45
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:48
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
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: ReadSpecThirdOrderPiecewisePolynomial.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:972
string