diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp index 5d9ad8b0071f..b0b535844e27 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp @@ -403,6 +403,11 @@ struct ReceiveMortarDataAndApplyOperator< Frame::Inertial>>(box), db::get>(box), + db::get< + domain::Tags::DetJacobian>( + box), + db::get>(box), db::get>>(box), db::get>(box), db::get>(box), + db::get< + domain::Tags::DetJacobian>( + box), + db::get>(box), db::get>>(box), db::get>>( box), diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp index 5a7182a60983..2f30a92117b8 100644 --- a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp @@ -596,6 +596,9 @@ struct DgOperatorImpl, const InverseJacobian& inv_jacobian, const Scalar& det_inv_jacobian, + const Scalar& det_jacobian, + const InverseJacobian& det_times_inv_jacobian, const DirectionMap>& face_normal_magnitudes, const DirectionMap>& face_jacobians, const DirectionMap, if (local_data_is_zero) { operator_applied_to_vars->initialize(num_points, 0.); } else { - ASSERT(formulation == ::dg::Formulation::StrongInertial, - "Only the strong formulation is currently implemented."); - divergence(operator_applied_to_vars, primal_fluxes, mesh, inv_jacobian); - // This is the sign flip that makes the operator _minus_ the Laplacian for - // a Poisson system - *operator_applied_to_vars *= -1.; + // "Massive" operators retain the factors from the volume integral: + // \int_volume div(F) \phi_p = w_p det(J)_p div(F)_p + // Here, `w` are the quadrature weights (the diagonal logical mass matrix + // with mass-lumping) and det(J) is the Jacobian determinant. The + // quantities are evaluated at the grid point `p`. + if (formulation == ::dg::Formulation::StrongInertial) { + // Compute strong divergence: + // div(F) = (J^\hat{i}_i)_p \sum_q (D_\hat{i})_pq (F^i)_q. + divergence(operator_applied_to_vars, primal_fluxes, mesh, + massive ? det_times_inv_jacobian : inv_jacobian); + // This is the sign flip that makes the operator _minus_ the Laplacian + // for a Poisson system + *operator_applied_to_vars *= -1.; + } else { + // Compute weak divergence: + // F^i \partial_i \phi = 1/w_p \sum_q + // (D^T_\hat{i})_pq (w det(J) J^\hat{i}_i F^i)_q + weak_divergence(operator_applied_to_vars, primal_fluxes, mesh, + det_times_inv_jacobian); + if (not massive) { + *operator_applied_to_vars *= get(det_inv_jacobian); + } + } if constexpr (not std::is_same_v) { + Variables> sources{num_points, 0.}; std::apply( - [&operator_applied_to_vars, &primal_vars, + [&sources, &primal_vars, &primal_fluxes](const auto&... expanded_sources_args) { SourcesComputer::apply( - make_not_null( - &get(*operator_applied_to_vars))..., + make_not_null(&get(sources))..., expanded_sources_args..., get(primal_vars)..., get(primal_fluxes)...); }, sources_args); + if (massive) { + sources *= get(det_jacobian); + } + *operator_applied_to_vars += sources; } - if (massive) { - *operator_applied_to_vars /= get(det_inv_jacobian); - ::dg::apply_mass_matrix(operator_applied_to_vars, mesh); - } + } + if (massive) { + ::dg::apply_mass_matrix(operator_applied_to_vars, mesh); } // Add boundary corrections @@ -802,18 +825,20 @@ struct DgOperatorImpl, get>( remote_data.field_data))); primal_boundary_corrections_on_mortar *= penalty; - const auto add_remote_avg_contribution = [](auto& lhs, const auto& rhs) { + const auto add_avg_contribution = [](auto& lhs, const auto& rhs, + const double prefactor) { for (size_t i = 0; i < lhs.size(); ++i) { - lhs[i] += 0.5 * rhs[i]; + lhs[i] += prefactor * rhs[i]; } }; - EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution( + EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution( get(primal_boundary_corrections_on_mortar), - get<::Tags::NormalDotFlux>(local_data.field_data))); - EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution( + get<::Tags::NormalDotFlux>(local_data.field_data), + formulation == ::dg::Formulation::StrongInertial ? 0.5 : -0.5)); + EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution( get(primal_boundary_corrections_on_mortar), - get<::Tags::NormalDotFlux>( - remote_data.field_data))); + get<::Tags::NormalDotFlux>(remote_data.field_data), + 0.5)); // Project from the mortar back down to the face if needed, lift and add // to operator. See auxiliary boundary corrections above for details. @@ -913,6 +938,9 @@ struct DgOperatorImpl, const InverseJacobian& inv_jacobian, const Scalar& det_inv_jacobian, + const Scalar& det_jacobian, + const InverseJacobian& det_times_inv_jacobian, const DirectionMap>& face_normals, const DirectionMap>& face_normal_vectors, const DirectionMap>& face_normal_magnitudes, @@ -965,10 +993,10 @@ struct DgOperatorImpl, 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, face_normal_magnitudes, - face_jacobians, face_jacobian_times_inv_jacobians, all_mortar_meshes, - all_mortar_sizes, {}, penalty_parameter, massive, formulation, - temporal_id, sources_args); + element, mesh, inv_jacobian, det_inv_jacobian, det_jacobian, + det_times_inv_jacobian, face_normal_magnitudes, face_jacobians, + face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes, + {}, penalty_parameter, massive, formulation, temporal_id, sources_args); // Impose the nonlinear (constant) boundary contribution as fixed sources on // the RHS of the equations *fixed_sources -= operator_applied_to_zero_vars; diff --git a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp index 38b856ffea6b..4a6eea012460 100644 --- a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp @@ -158,6 +158,9 @@ struct SubdomainOperator domain::Tags::InverseJacobian, domain::Tags::DetInvJacobian, + domain::Tags::DetJacobian, + domain::Tags::DetTimesInvJacobian, domain::Tags::Faces>, domain::Tags::Faces; - static type create_from_options(const type value) { - if (not(value == ::dg::Formulation::StrongInertial)) { - ERROR_NO_TRACE("Only the strong formulation is currently supported."); - } - return value; - } + static type create_from_options(const type value) { return value; } }; } // namespace Tags diff --git a/support/Pipelines/Bbh/InitialData.yaml b/support/Pipelines/Bbh/InitialData.yaml index afce7c4ce3b1..3afd3eb95265 100644 --- a/support/Pipelines/Bbh/InitialData.yaml +++ b/support/Pipelines/Bbh/InitialData.yaml @@ -85,7 +85,7 @@ Discretization: PenaltyParameter: 1. Massive: True Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "BbhVolume" diff --git a/tests/InputFiles/Poisson/Lorentzian.yaml b/tests/InputFiles/Poisson/Lorentzian.yaml new file mode 100644 index 000000000000..24104621cc17 --- /dev/null +++ b/tests/InputFiles/Poisson/Lorentzian.yaml @@ -0,0 +1,117 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +Executable: SolvePoisson3D +Testing: + Check: parse;execute_check_output + Timeout: 40 + Priority: High +ExpectedOutput: + - LorentzianReductions.h5 + - LorentzianVolume0.h5 +OutputFileChecks: + - Label: Discretization error + Subfile: /ErrorNorms.dat + FileGlob: LorentzianReductions.h5 + SkipColumns: [0, 1, 2] + AbsoluteTolerance: 0.05 + +--- + +Background: &solution Lorentzian + +InitialGuess: Zero + +RandomizeInitialGuess: None + +DomainCreator: + Sphere: + InnerRadius: 10 + OuterRadius: &outer_radius 1e9 + Interior: + FillWithSphericity: 0. + InitialRefinement: 0 + InitialGridPoints: 6 + UseEquiangularMap: True + EquatorialCompression: None + RadialPartitioning: [&outer_shell_inner_radius 20.] + RadialDistribution: [Linear, &outer_shell_distribution Inverse] + WhichWedges: All + TimeDependentMaps: None + OuterBoundaryCondition: + AnalyticSolution: + Solution: *solution + Field: Dirichlet + +Discretization: + DiscontinuousGalerkin: + PenaltyParameter: 1. + Massive: True + Quadrature: GaussLobatto + Formulation: WeakInertial + +Observers: + VolumeFileName: "LorentzianVolume" + ReductionFileName: "LorentzianReductions" + +LinearSolver: + ConvergenceCriteria: + MaxIterations: 10 + RelativeResidual: 0. + AbsoluteResidual: 1.e-3 + Verbosity: Quiet + + Multigrid: + Iterations: 1 + MaxLevels: Auto + PreSmoothing: True + PostSmoothingAtBottom: False + Verbosity: Silent + OutputVolumeData: True + + SchwarzSmoother: + Iterations: 3 + MaxOverlap: 2 + Verbosity: Silent + SubdomainSolver: + ExplicitInverse: + WriteMatrixToFile: None + ObservePerCoreReductions: False + +RadiallyCompressedCoordinates: + InnerRadius: *outer_shell_inner_radius + OuterRadius: *outer_radius + Compression: *outer_shell_distribution + +EventsAndTriggers: + - Trigger: HasConverged + Events: + - ObserveNorms: + SubfileName: ErrorNorms + TensorsToObserve: + - Name: Error(Field) + NormType: L2Norm + Components: Sum + - ObserveFields: + SubfileName: VolumeData + VariablesToObserve: + - Field + - Error(Field) + - RadiallyCompressedCoordinates + - FixedSource(Field) + InterpolateToMesh: None + CoordinatesFloatingPointType: Double + FloatingPointTypes: [Double] + +PhaseChangeAndTriggers: + +BuildMatrix: + MatrixSubfileName: Matrix + Verbosity: Verbose + +Parallelization: + ElementDistribution: NumGridPoints + +ResourceInfo: + AvoidGlobalProc0: false + Singletons: Auto diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml index 63f200b14964..117fdc9681e3 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml +++ b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml @@ -59,7 +59,7 @@ Discretization: PenaltyParameter: 1. Massive: False Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "PoissonProductOfSinusoids3DVolume" diff --git a/tests/InputFiles/Xcts/BinaryBlackHole.yaml b/tests/InputFiles/Xcts/BinaryBlackHole.yaml index 8cd288e72168..9304a80cf900 100644 --- a/tests/InputFiles/Xcts/BinaryBlackHole.yaml +++ b/tests/InputFiles/Xcts/BinaryBlackHole.yaml @@ -117,7 +117,7 @@ Discretization: PenaltyParameter: 1. Massive: True Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "BbhVolume" diff --git a/tests/InputFiles/Xcts/HeadOnBns.yaml b/tests/InputFiles/Xcts/HeadOnBns.yaml index 520ed737f50e..b42e2661589e 100644 --- a/tests/InputFiles/Xcts/HeadOnBns.yaml +++ b/tests/InputFiles/Xcts/HeadOnBns.yaml @@ -108,7 +108,7 @@ Discretization: PenaltyParameter: 1. Massive: True Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "BnsVolume" diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml index 174c6121d1ab..18254ff5b843 100644 --- a/tests/InputFiles/Xcts/KerrSchild.yaml +++ b/tests/InputFiles/Xcts/KerrSchild.yaml @@ -72,7 +72,7 @@ Discretization: PenaltyParameter: 1. Massive: True Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "KerrSchildVolume" diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml index 9c5fcaeef299..59549741a3f6 100644 --- a/tests/InputFiles/Xcts/TovStar.yaml +++ b/tests/InputFiles/Xcts/TovStar.yaml @@ -62,7 +62,7 @@ Discretization: PenaltyParameter: 1. Massive: True Quadrature: GaussLobatto - Formulation: StrongInertial + Formulation: WeakInertial Observers: VolumeFileName: "TovStarVolume"