11 #include "DataStructures/DataBox/Tag.hpp"
12 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
18 #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
19 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
20 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
24 #include "Utilities/ProtocolHelpers.hpp"
145 template <
size_t Dim,
typename FluxesComputerTag,
typename FieldTagsList,
146 typename AuxiliaryFieldTagsList>
149 template <
size_t Dim,
typename FluxesComputerTag,
typename... FieldTags,
150 typename... AuxiliaryFieldTags>
152 tmpl::list<FieldTags...>,
153 tmpl::list<AuxiliaryFieldTags...>>
156 using fluxes_computer_tag = FluxesComputerTag;
157 using FluxesComputer =
typename fluxes_computer_tag::type;
167 template <
typename Tag>
173 template <
typename Tag>
180 struct PenaltyParameter {
183 "The prefactor to the penalty term of the flux. Values closer to one "
184 "lead to better-conditioned problems, but on curved meshes the penalty "
185 "parameter may need to be increased to keep the problem well-defined."};
186 static type lower_bound() noexcept {
return 1.; }
188 using options = tmpl::list<PenaltyParameter>;
190 "The internal penalty flux for elliptic systems."};
191 static std::string name() noexcept {
return "InternalPenalty"; }
193 FirstOrderInternalPenalty() =
default;
194 explicit FirstOrderInternalPenalty(
const double penalty_parameter) noexcept
195 : penalty_parameter_(penalty_parameter) {}
198 void pup(PUP::er& p) noexcept { p | penalty_parameter_; }
200 using variables_tags = tmpl::list<FieldTags..., AuxiliaryFieldTags...>;
210 fluxes_computer_tag>;
213 fluxes_computer_tag>;
215 using package_field_tags =
216 tmpl::list<::Tags::NormalDotFlux<AuxiliaryFieldTags>...,
217 NormalDotDivFlux<AuxiliaryFieldTags>...,
218 NormalDotNormalDotFlux<AuxiliaryFieldTags>..., ElementSize>;
219 using package_extra_tags = tmpl::list<PerpendicularNumPoints>;
221 template <
typename... FluxesArgs>
224 AuxiliaryFieldTags>::type*>... packaged_n_dot_aux_fluxes,
226 AuxiliaryFieldTags>::type*>... packaged_n_dot_div_aux_fluxes,
228 typename NormalDotNormalDotFlux<AuxiliaryFieldTags>::
229 type*>... packaged_n_dot_principal_n_dot_aux_fluxes,
234 const typename ::Tags::NormalDotFlux<
235 AuxiliaryFieldTags>::type&... normal_dot_auxiliary_field_fluxes,
236 const typename ::Tags::div<
237 ::
Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
239 const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
240 const FluxesComputer& fluxes_computer,
241 const FluxesArgs&... fluxes_args)
const noexcept {
242 *perpendicular_num_points = volume_mesh.extents(direction.dimension());
243 get(*element_size) = 2. /
get(face_normal_magnitude);
244 auto principal_div_aux_field_fluxes =
make_with_value<Variables<tmpl::list<
247 auto principal_ndot_aux_field_fluxes =
make_with_value<Variables<tmpl::list<
250 fluxes_computer.apply(
253 principal_div_aux_field_fluxes))...,
254 fluxes_args..., div_auxiliary_field_fluxes...);
255 fluxes_computer.apply(
258 principal_ndot_aux_field_fluxes))...,
259 fluxes_args..., normal_dot_auxiliary_field_fluxes...);
260 const auto apply_normal_dot =
261 [&principal_div_aux_field_fluxes, &principal_ndot_aux_field_fluxes,
262 &interface_unit_normal](
263 auto field_tag_v,
auto packaged_normal_dot_auxiliary_field_flux,
264 auto packaged_normal_dot_div_auxiliary_field_flux,
265 auto packaged_normal_dot_principal_ndot_aux_field_flux,
266 const auto& normal_dot_auxiliary_field_flux) noexcept {
267 using field_tag = decltype(field_tag_v);
269 *packaged_normal_dot_auxiliary_field_flux =
270 normal_dot_auxiliary_field_flux;
273 packaged_normal_dot_div_auxiliary_field_flux,
274 interface_unit_normal,
276 principal_div_aux_field_fluxes));
278 packaged_normal_dot_principal_ndot_aux_field_flux,
279 interface_unit_normal,
281 principal_ndot_aux_field_fluxes));
284 FieldTags{}, packaged_n_dot_aux_fluxes, packaged_n_dot_div_aux_fluxes,
285 packaged_n_dot_principal_n_dot_aux_fluxes,
286 normal_dot_auxiliary_field_fluxes));
291 FieldTags>::type*>... numerical_flux_for_fields,
293 AuxiliaryFieldTags>::type*>... numerical_flux_for_auxiliary_fields,
294 const typename ::Tags::NormalDotFlux<
295 AuxiliaryFieldTags>::type&... normal_dot_auxiliary_flux_interiors,
296 const typename NormalDotDivFlux<
297 AuxiliaryFieldTags>::type&... normal_dot_div_auxiliary_flux_interiors,
298 const typename NormalDotNormalDotFlux<
299 AuxiliaryFieldTags>::type&... ndot_ndot_aux_flux_interiors,
301 const size_t perpendicular_num_points_interior,
302 const typename ::Tags::NormalDotFlux<AuxiliaryFieldTags>::
303 type&... minus_normal_dot_auxiliary_flux_exteriors,
304 const typename NormalDotDivFlux<
305 AuxiliaryFieldTags>::type&... minus_normal_dot_div_aux_flux_exteriors,
306 const typename NormalDotDivFlux<
307 AuxiliaryFieldTags>::type&... ndot_ndot_aux_flux_exteriors,
309 const size_t perpendicular_num_points_exterior)
const noexcept {
310 const auto penalty = ::elliptic::dg::penalty(
311 min(
get(element_size_interior),
get(element_size_exterior)),
312 std::max(perpendicular_num_points_interior,
313 perpendicular_num_points_exterior),
316 const auto assemble_numerical_fluxes = [&
penalty](
317 const auto numerical_flux_for_field,
318 const auto numerical_flux_for_auxiliary_field,
319 const auto& normal_dot_auxiliary_flux_interior,
320 const auto& normal_dot_div_auxiliary_flux_interior,
321 const auto& ndot_ndot_aux_flux_interior,
322 const auto& minus_normal_dot_auxiliary_flux_exterior,
323 const auto& minus_normal_dot_div_aux_flux_exterior,
324 const auto& ndot_ndot_aux_flux_exterior) noexcept {
325 for (
auto it = numerical_flux_for_auxiliary_field->begin();
326 it != numerical_flux_for_auxiliary_field->end(); it++) {
328 numerical_flux_for_auxiliary_field->get_tensor_index(it);
334 *it = 0.5 * (normal_dot_auxiliary_flux_interior.get(index) -
335 minus_normal_dot_auxiliary_flux_exterior.get(index));
337 for (
auto it = numerical_flux_for_field->begin();
338 it != numerical_flux_for_field->end(); it++) {
339 const auto index = numerical_flux_for_field->get_tensor_index(it);
340 *it = 0.5 * (normal_dot_div_auxiliary_flux_interior.get(index) -
341 minus_normal_dot_div_aux_flux_exterior.get(index)) -
342 penalty * (ndot_ndot_aux_flux_interior.get(index) -
343 ndot_ndot_aux_flux_exterior.get(index));
347 numerical_flux_for_fields, numerical_flux_for_auxiliary_fields,
348 normal_dot_auxiliary_flux_interiors,
349 normal_dot_div_auxiliary_flux_interiors, ndot_ndot_aux_flux_interiors,
350 minus_normal_dot_auxiliary_flux_exteriors,
351 minus_normal_dot_div_aux_flux_exteriors, ndot_ndot_aux_flux_exteriors));
362 template <
typename... FluxesArgs>
363 void compute_dirichlet_boundary(
365 FieldTags>::type*>... numerical_flux_for_fields,
367 AuxiliaryFieldTags>::type*>... numerical_flux_for_auxiliary_fields,
368 const typename FieldTags::type&... dirichlet_fields,
370 const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
372 const FluxesComputer& fluxes_computer,
373 const FluxesArgs&... fluxes_args)
const noexcept {
374 const auto penalty = ::elliptic::dg::penalty(
375 2. /
get(face_normal_magnitude),
376 volume_mesh.extents(direction.dimension()), penalty_parameter_);
379 auto dirichlet_auxiliary_field_fluxes =
382 interface_unit_normal,
384 fluxes_computer.apply(
387 dirichlet_auxiliary_field_fluxes))...,
388 fluxes_args..., dirichlet_fields...);
389 const auto apply_normal_dot_aux =
390 [&interface_unit_normal, &dirichlet_auxiliary_field_fluxes ](
391 auto auxiliary_field_tag_v,
392 const auto numerical_flux_for_auxiliary_field) noexcept {
393 using auxiliary_field_tag = decltype(auxiliary_field_tag_v);
395 numerical_flux_for_auxiliary_field, interface_unit_normal,
400 AuxiliaryFieldTags{}, numerical_flux_for_auxiliary_fields));
403 auto principal_dirichlet_auxiliary_field_fluxes =
406 interface_unit_normal,
408 fluxes_computer.apply(
411 principal_dirichlet_auxiliary_field_fluxes))...,
412 fluxes_args..., *numerical_flux_for_auxiliary_fields...);
413 const auto assemble_dirichlet_penalty = [
414 &interface_unit_normal, &
penalty,
415 &principal_dirichlet_auxiliary_field_fluxes
416 ](
auto field_tag_v,
const auto numerical_flux_for_field) noexcept {
417 using field_tag = decltype(field_tag_v);
419 numerical_flux_for_field, interface_unit_normal,
421 principal_dirichlet_auxiliary_field_fluxes));
422 for (
auto it = numerical_flux_for_field->begin();
423 it != numerical_flux_for_field->end(); it++) {
428 assemble_dirichlet_penalty(FieldTags{}, numerical_flux_for_fields));