diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp index d5b0af6eb574..a832daa2e194 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp @@ -28,6 +28,7 @@ #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp" #include "Elliptic/DiscontinuousGalerkin/Tags.hpp" #include "Elliptic/Systems/GetFluxesComputer.hpp" +#include "Elliptic/Systems/GetModifyBoundaryData.hpp" #include "Elliptic/Systems/GetSourcesComputer.hpp" #include "Elliptic/Utilities/ApplyAt.hpp" #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp" @@ -498,6 +499,9 @@ template struct DgOperator { + using system = System; + using temporal_id_tag = TemporalIdTag; + private: static constexpr size_t Dim = System::volume_dim; @@ -530,15 +534,17 @@ template , typename SourcesArgsTags = - elliptic::get_sources_argument_tags> + elliptic::get_sources_argument_tags, + typename ModifyBoundaryDataArgsTags = + elliptic::get_modify_boundary_data_args_tags> struct ImposeInhomogeneousBoundaryConditionsOnSource; /// \cond template + typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags> struct ImposeInhomogeneousBoundaryConditionsOnSource< System, FixedSourcesTag, tmpl::list, - tmpl::list> { + tmpl::list, tmpl::list> { private: static constexpr size_t Dim = System::volume_dim; using BoundaryConditionsBase = typename System::boundary_conditions_base; @@ -636,12 +642,16 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource< Frame::Inertial>>>(box), db::get<::Tags::Mortars, Dim>>(box), db::get<::Tags::Mortars<::Tags::MortarSize, Dim>>(box), + db::get<::Tags::Mortars, + Dim>>(box), db::get<::Tags::Mortars>(box), db::get(box), db::get(box), apply_boundary_condition, std::forward_as_tuple(db::get(box)...), std::forward_as_tuple(db::get(box)...), - fluxes_args_on_faces); + fluxes_args_on_faces, + std::forward_as_tuple(db::get(box)...)); // Move the mutated data back into the DataBox db::mutate( diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp index edcd681eebd0..0c66e6a343cb 100644 --- a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp @@ -517,11 +517,11 @@ struct DgOperatorImpl, // --- This is essentially a break to communicate the mortar data --- - template static void apply_operator( const gsl::not_null>*> @@ -686,11 +686,6 @@ struct DgOperatorImpl, const auto& neighbor_id = mortar_id.id(); const bool is_internal = (neighbor_id != ElementId::external_boundary_id()); - if constexpr (AllDataIsZero) { - if (is_internal) { - continue; - } - } if (not directions_predicate(direction)) { continue; } @@ -992,6 +987,7 @@ struct DgOperatorImpl, template , face_jacobian_times_inv_jacobians, const ::dg::MortarMap>& all_mortar_meshes, const ::dg::MortarMap>& all_mortar_sizes, + const ::dg::MortarMap>& mortar_jacobians, const ::dg::MortarMap>& penalty_factors, const bool massive, const ::dg::Formulation formulation, const ApplyBoundaryCondition& apply_boundary_condition, const std::tuple& fluxes_args, const std::tuple& sources_args, - const DirectionMap>& - fluxes_args_on_faces) { + const DirectionMap>& fluxes_args_on_faces, + const std::tuple& modify_boundary_data_args) { // We just feed zero variables through the nonlinear operator to extract the // constant contribution at external boundaries. Since the variables are // zero the operator simplifies quite a lot. The simplification is probably // not very important for performance because this function will only be // called when solving a linear elliptic system and only once during // initialization, but we specialize the operator for zero data nonetheless - // just so we can ignore internal boundaries. For internal boundaries we - // would unnecessarily have to copy mortar data around to emulate the - // communication step, so by just skipping internal boundaries we avoid - // that. + // just so we can ignore internal boundaries. Only when the system modifies + // boundary data do we need to handle internal boundaries (see below), + // otherwise we can skip them. const size_t num_points = mesh.number_of_grid_points(); const Variables> zero_primal_vars{num_points, 0.}; @@ -1042,7 +1038,7 @@ struct DgOperatorImpl, unused_deriv_vars_buffer{}; Variables> operator_applied_to_zero_vars{ num_points}; - // Set up data on mortars. We only need them at external boundaries. + // Set up data on mortars ::dg::MortarMap, tmpl::list>> all_mortar_data{}; @@ -1054,15 +1050,53 @@ struct DgOperatorImpl, element, mesh, inv_jacobian, face_normals, all_mortar_meshes, all_mortar_sizes, temporal_id, apply_boundary_condition, fluxes_args); - apply_operator( - make_not_null(&operator_applied_to_zero_vars), - make_not_null(&all_mortar_data), zero_primal_vars, primal_fluxes_buffer, - element, mesh, inv_jacobian, det_inv_jacobian, det_jacobian, - det_times_inv_jacobian, face_normals, face_normal_vectors, - face_normal_magnitudes, face_jacobians, - face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes, - {}, penalty_factors, massive, formulation, temporal_id, - fluxes_args_on_faces, sources_args); + // Modify internal boundary data if needed, e.g. to transform from one + // variable to another when crossing element boundaries. This is a nonlinear + // operation in the sense that feeding zero through the operator is nonzero, + // so we evaluate it here and add the contribution to the fixed sources, + // just like inhomogeneous external boundary conditions. + if constexpr (not std::is_same_v) { + for (const auto& [direction, neighbors] : element.neighbors()) { + for (const auto& neighbor_id : neighbors) { + const ::dg::MortarId mortar_id{direction, neighbor_id}; + ASSERT(not all_mortar_data.contains(mortar_id), + "Mortar data with ID " << mortar_id << " already exists."); + auto& mortar_data = all_mortar_data[mortar_id]; + // Manufacture zero mortar data, store as local data, then treat as + // received from neighbor and apply modifications + BoundaryData, tmpl::list> + remote_boundary_data_on_mortar{}; + remote_boundary_data_on_mortar.field_data.initialize( + all_mortar_meshes.at(mortar_id).number_of_grid_points(), 0.); + mortar_data.local_insert(temporal_id, remote_boundary_data_on_mortar); + // Modify "received" mortar data + std::apply( + [&remote_boundary_data_on_mortar, + &mortar_id](const auto&... args) { + System::modify_boundary_data::apply( + make_not_null(&get( + remote_boundary_data_on_mortar.field_data))..., + make_not_null(&get<::Tags::NormalDotFlux>( + remote_boundary_data_on_mortar.field_data))..., + mortar_id, args...); + }, + modify_boundary_data_args); + // Insert as remote mortar data + mortar_data.remote_insert(temporal_id, + std::move(remote_boundary_data_on_mortar)); + } + } + } + apply_operator(make_not_null(&operator_applied_to_zero_vars), + make_not_null(&all_mortar_data), zero_primal_vars, + primal_fluxes_buffer, element, mesh, inv_jacobian, + det_inv_jacobian, det_jacobian, det_times_inv_jacobian, + face_normals, face_normal_vectors, face_normal_magnitudes, + face_jacobians, face_jacobian_times_inv_jacobians, + all_mortar_meshes, all_mortar_sizes, mortar_jacobians, + penalty_factors, massive, formulation, temporal_id, + fluxes_args_on_faces, sources_args); // Impose the nonlinear (constant) boundary contribution as fixed sources on // the RHS of the equations *fixed_sources -= operator_applied_to_zero_vars; @@ -1094,7 +1128,7 @@ void prepare_mortar_data(Args&&... args) { */ template void apply_operator(Args&&... args) { - detail::DgOperatorImpl::template apply_operator( + detail::DgOperatorImpl::apply_operator( std::forward(args)...); } diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp index 12b3a960e2b2..ef792a7720d3 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp @@ -274,12 +274,14 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { // external boundaries. ::Tags::deriv, tmpl::size_t, Frame::Inertial>>>, - tmpl::list<::Tags::Mortars, Dim>, - ::Tags::Mortars<::Tags::MortarSize, Dim>, - ::Tags::Mortars, - Dim>, - ::Tags::Mortars>>; + tmpl::list< + ::Tags::Mortars, Dim>, + ::Tags::Mortars<::Tags::MortarSize, Dim>, + ::Tags::Mortars, Dim>, + ::Tags::Mortars, + Dim>, + ::Tags::Mortars>>; using argument_tags = tmpl::append< tmpl::list, domain::Tags::Element, domain::Tags::NeighborMesh, domain::Tags::ElementMap, @@ -318,6 +320,8 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { const gsl::not_null<::dg::MortarMap>*> mortar_meshes, const gsl::not_null<::dg::MortarMap>*> mortar_sizes, + const gsl::not_null<::dg::MortarMap>*> + all_mortar_inertial_coords, const gsl::not_null<::dg::MortarMap>*> mortar_jacobians, const gsl::not_null<::dg::MortarMap>*> @@ -333,9 +337,10 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { apply(face_directions, faces_inertial_coords, face_normals, face_normal_vectors, face_normal_magnitudes, face_jacobians, face_jacobian_times_inv_jacobian, deriv_unnormalized_face_normals, - mortar_meshes, mortar_sizes, mortar_jacobians, penalty_factors, mesh, - element, neighbor_meshes, element_map, inv_jacobian, domain, - functions_of_time, penalty_parameter, nullptr, nullptr, amr_data...); + mortar_meshes, mortar_sizes, all_mortar_inertial_coords, + mortar_jacobians, penalty_factors, mesh, element, neighbor_meshes, + element_map, inv_jacobian, domain, functions_of_time, + penalty_parameter, nullptr, nullptr, amr_data...); } template static void apply( @@ -359,6 +364,8 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { const gsl::not_null<::dg::MortarMap>*> mortar_meshes, const gsl::not_null<::dg::MortarMap>*> mortar_sizes, + const gsl::not_null<::dg::MortarMap>*> + all_mortar_inertial_coords, const gsl::not_null<::dg::MortarMap>*> mortar_jacobians, const gsl::not_null<::dg::MortarMap>*> @@ -450,6 +457,7 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { // Mortars (internal directions) mortar_meshes->clear(); mortar_sizes->clear(); + all_mortar_inertial_coords->clear(); mortar_jacobians->clear(); penalty_factors->clear(); const auto& element_id = element.id(); @@ -471,8 +479,11 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { const auto& mortar_size = mortar_sizes->at(mortar_id); const auto mortar_logical_coords = detail::mortar_logical_coordinates( mortar_mesh, mortar_size, direction); - const auto mortar_inertial_coords = - element_map(mortar_logical_coords, 0., functions_of_time); + const auto& mortar_inertial_coords = + all_mortar_inertial_coords + ->emplace(mortar_id, element_map(mortar_logical_coords, 0., + functions_of_time)) + .first->second; Scalar perpendicular_element_size{}; if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) { auto& mortar_jacobian = (*mortar_jacobians)[mortar_id]; diff --git a/src/Elliptic/Protocols/FirstOrderSystem.hpp b/src/Elliptic/Protocols/FirstOrderSystem.hpp index 8355d9b58493..2cb218e1b460 100644 --- a/src/Elliptic/Protocols/FirstOrderSystem.hpp +++ b/src/Elliptic/Protocols/FirstOrderSystem.hpp @@ -154,7 +154,7 @@ struct test_fields_and_fluxes, * The function is expected to _add_ the sources \f$S_\alpha\f$ to the * output buffers. It must also have the following alias: * - * - `const_global_cache_tags`: the subset of `argument_tags` that can be + * - `const_global_cache_tags`: the subset of `argument_tags` that can be * retrieved from _any_ element's DataBox, because they are stored in the * global cache. * @@ -164,6 +164,35 @@ struct test_fields_and_fluxes, * boundary conditions. Boundary conditions can be factory-created from this * base class. Currently this should be a specialization of * `elliptic::BoundaryConditions::BoundaryCondition`. + * + * - `modify_boundary_data` (optional, can be `void`): A class that can modify + * boundary data received from a neighboring element, e.g. to transform from + * one variable to another across element boundaries. This can be used to + * solve for different variables in some regions of the domain to absorb + * singularities. For example, when solving $-\Delta u = f$ we could define + * $u=u_R + u_P$ in some region of the domain with a known (often singular) + * field $u_P$ and solve only for the regular field $u_R$ in this region. So, + * when receiving boundary data from outside this region, we subtract $u_P$, + * and when receiving boundary data from inside this region, we add $u_P$. We + * also add $\Delta u_P$ to the volume fixed sources $f$ inside the + * regularized region. + * + * The `modify_boundary_data` must have an `argument_tags` type alias and an + * `apply` function that takes these arguments in this order: + * + * 1. The primal fields and the normal-dot-fluxes on the mortar as not-null + * pointer. These hold the received data from the neighbor and can be + * modified. Note: the normal-dot-fluxes have been computed with the + * neighbor's normal, so from the perspective of the receiving element they + * are $-n_i F^i$ where $n_i$ is the face normal of the receiving element. + * 2. The `const DirectionalId& mortar_id` identifying the mortar. + * 3. The `argument_tags`. + * + * Currently, modification made by this function must not depend on the + * variables, meaning that the modification can only be adding or subtracting + * a precomputed field. This is a simplification so the linearized operator is + * not modified at all and can be relaxed if needed (then we also need + * `modify_boundary_data_linearized`). */ struct FirstOrderSystem { template @@ -191,6 +220,7 @@ struct FirstOrderSystem { static_assert( std::is_base_of_v); + using modify_boundary_data = typename ConformingType::modify_boundary_data; }; }; diff --git a/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp b/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp index 99da70e77018..f7988b9095d1 100644 --- a/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp @@ -94,5 +94,6 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition<3>; + using modify_boundary_data = void; }; } // namespace BnsInitialData diff --git a/src/Elliptic/Systems/Elasticity/FirstOrderSystem.hpp b/src/Elliptic/Systems/Elasticity/FirstOrderSystem.hpp index 7064278a563d..9b309a994268 100644 --- a/src/Elliptic/Systems/Elasticity/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/Elasticity/FirstOrderSystem.hpp @@ -54,5 +54,6 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition; + using modify_boundary_data = void; }; } // namespace Elasticity diff --git a/src/Elliptic/Systems/GetModifyBoundaryData.hpp b/src/Elliptic/Systems/GetModifyBoundaryData.hpp new file mode 100644 index 000000000000..b49ee26cee18 --- /dev/null +++ b/src/Elliptic/Systems/GetModifyBoundaryData.hpp @@ -0,0 +1,24 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "Utilities/TMPL.hpp" + +namespace elliptic { +namespace detail { +struct NoModifyBoundaryData { + using argument_tags = tmpl::list<>; +}; +} // namespace detail + +/// The `argument_tags` of the `System::modify_boundary_data`, or an empty list +/// if `System::modify_boundary_data` is `void`. +template +using get_modify_boundary_data_args_tags = typename tmpl::conditional_t< + std::is_same_v, + detail::NoModifyBoundaryData, + typename System::modify_boundary_data>::argument_tags; +} // namespace elliptic diff --git a/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp b/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp index 43282d66da88..6b90cf226f49 100644 --- a/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/Poisson/FirstOrderSystem.hpp @@ -94,5 +94,6 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition; + using modify_boundary_data = void; }; } // namespace Poisson diff --git a/src/Elliptic/Systems/Punctures/FirstOrderSystem.hpp b/src/Elliptic/Systems/Punctures/FirstOrderSystem.hpp index 1d8260ca8145..61bf3637416c 100644 --- a/src/Elliptic/Systems/Punctures/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/Punctures/FirstOrderSystem.hpp @@ -42,6 +42,7 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition<3>; + using modify_boundary_data = void; }; } // namespace Punctures diff --git a/src/Elliptic/Systems/ScalarGaussBonnet/FirstOrderSystem.hpp b/src/Elliptic/Systems/ScalarGaussBonnet/FirstOrderSystem.hpp index 0298db100f6a..09d799a17401 100644 --- a/src/Elliptic/Systems/ScalarGaussBonnet/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/ScalarGaussBonnet/FirstOrderSystem.hpp @@ -47,6 +47,7 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition<3>; + using modify_boundary_data = void; }; } // namespace sgb diff --git a/src/Elliptic/Systems/Xcts/FirstOrderSystem.hpp b/src/Elliptic/Systems/Xcts/FirstOrderSystem.hpp index d7e37c04f750..78cedbca3292 100644 --- a/src/Elliptic/Systems/Xcts/FirstOrderSystem.hpp +++ b/src/Elliptic/Systems/Xcts/FirstOrderSystem.hpp @@ -230,6 +230,7 @@ struct FirstOrderSystem using boundary_conditions_base = elliptic::BoundaryConditions::BoundaryCondition<3>; + using modify_boundary_data = void; }; } // namespace Xcts diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp index 5e2d863c6633..811ada859772 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp @@ -251,6 +251,91 @@ struct Metavariables { void pup(PUP::er& /*p*/) {} }; +template +struct ModifiedPoissonSystem + : Poisson::FirstOrderSystem { + struct modify_boundary_data { + using argument_tags = tmpl::list< + domain::Tags::Element, + ::Tags::Mortars, Dim>>; + static void apply( + const gsl::not_null*> field, + const gsl::not_null*> normal_dot_flux, + const DirectionalId& mortar_id, const Element& element, + const DirectionalIdMap>& + all_mortar_coords) { + // Assuming that in the region x > 0.5 (block 1) we decompose the field as + // u = u_R + u_P with u_P = 2x + 3y + const auto& element_id = element.id(); + const bool element_is_modified = element_id.block_id() == 1; + const bool neighbor_is_modified = mortar_id.id().block_id() == 1; + if (element_is_modified == neighbor_is_modified) { + return; + } + const auto& x = all_mortar_coords.at(mortar_id); + const DataVector singular_field = 2. * get<0>(x) + 3. * get<1>(x); + const auto& direction = mortar_id.direction(); + // Assuming rectilinear grid to simplify the normal vector + const double singular_normal_dot_flux = + (direction.dimension() == 0 ? 2. : 3.) * direction.sign(); + // In modified elements, the field is u_R = u - u_P + const double sign = element_is_modified ? -1. : 1.; + get(*field) += sign * singular_field; + // Minus sign because we are modifying the _received_ fluxes, which were + // computed using the neighbor's face normal + get(*normal_dot_flux) -= sign * singular_normal_dot_flux; + } + }; +}; + +template +struct ModifiedPoissonSolution : Poisson::Solutions::ProductOfSinusoids { + using Poisson::Solutions::ProductOfSinusoids::ProductOfSinusoids; + std::unique_ptr get_clone() + const override { + return std::make_unique(*this); + } + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(ModifiedPoissonSolution); // NOLINT + template + + tuples::TaggedTuple variables( + const tnsr::I& x, + tmpl::list /*meta*/) const { + auto vars = Poisson::Solutions::ProductOfSinusoids::variables( + x, tmpl::list{}); + // Using coordinates to determine in which block we are. This only works for + // Gauss points, because GaussLobatto points are duplicate on boundaries. + const bool element_is_modified = get<0>(x)[0] > 0.5; + if (not element_is_modified) { + return vars; + } + const DataVector singular_field = 2. * get<0>(x) + 3. * get<1>(x); + if constexpr (tmpl::list_contains_v, + Poisson::Tags::Field>) { + get(get>(vars)) -= singular_field; + } + using deriv_field = ::Tags::deriv, + tmpl::size_t, Frame::Inertial>; + if constexpr (tmpl::list_contains_v, + deriv_field>) { + get<0>(get(vars)) -= 2.; + get<1>(get(vars)) -= 3.; + } + using flux_tag = ::Tags::Flux, + tmpl::size_t, Frame::Inertial>; + if constexpr (tmpl::list_contains_v, + flux_tag>) { + get<0>(get(vars)) -= 2.; + get<1>(get(vars)) -= 3.; + } + return vars; + } +}; + +template +PUP::able::PUP_ID ModifiedPoissonSolution::my_PUP_ID = 0; // NOLINT + template < typename System, bool Linearized, typename AnalyticSolution, size_t Dim = System::volume_dim, @@ -1417,4 +1502,39 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") { analytic_solution_operator_approx, {}, true); } } + { + INFO("2D with modified boundary data"); + using system = ModifiedPoissonSystem<2>; + const ModifiedPoissonSolution<2> analytic_solution{{{M_PI, M_PI}}}; + const auto dirichlet_bc = + elliptic::BoundaryConditions::AnalyticSolution{ + analytic_solution.get_clone(), + elliptic::BoundaryConditionType::Dirichlet}; + // In block 1 (x > 0.5) we decompose u = U_R + u_P with u_P = 2x + 3y and + // solve for u_R. + const domain::creators::AlignedLattice<2> domain_creator{ + {{{0., 0.5, 1.}, {0., 1.}}}, + {{1, 1}}, + {{12, 12}}, + {}, + {}, + {}, + {{{{dirichlet_bc.get_clone(), dirichlet_bc.get_clone()}}, + {{dirichlet_bc.get_clone(), dirichlet_bc.get_clone()}}}}}; + Approx analytic_solution_aux_approx = + Approx::custom().epsilon(1.e-4).scale(1.); + Approx analytic_solution_operator_approx = + Approx::custom().epsilon(1.e-4).scale(1.); + for (const auto& [dg_formulation, massive] : + cartesian_product(make_array(::dg::Formulation::StrongInertial, + ::dg::Formulation::StrongLogical, + ::dg::Formulation::WeakInertial), + make_array(true, false))) { + test_dg_operator( + domain_creator, penalty_parameter, massive, + Spectral::Quadrature::Gauss, dg_formulation, analytic_solution, + analytic_solution_aux_approx, analytic_solution_operator_approx, {}, + true); + } + } }