Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support Gauss points in elliptic DG schemes #5608

Merged
merged 7 commits into from
Nov 2, 2023
Merged

Support Gauss points in elliptic DG schemes #5608

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
5328a58
Fix static cache of max num points
nilsvu Oct 30, 2023
41430f7
Add an overload to ProjectToBoundary
nilsvu Oct 25, 2023
7d60699
Add some more docs
nilsvu Oct 30, 2023
2d3c38a
Add some comments on geometric aliasing
nilsvu Oct 25, 2023
5177e18
Add 'Quadrature' option to elliptic domains
nilsvu Oct 30, 2023
b226b13
Add Gauss quadrature to elliptic domain init
nilsvu Oct 30, 2023
9c0a5e4
Support Gauss points in elliptic DG schemes
nilsvu Oct 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -598,6 +598,9 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
db::get<domain::Tags::Faces<
Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
db::get<domain::Tags::Faces<
Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
Frame::Inertial>>>(box),
db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
db::get<elliptic::dg::Tags::PenaltyParameter>(box),
Expand Down
4 changes: 3 additions & 1 deletion src/Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "Domain/Creators/Tags/InitialExtents.hpp"
#include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
#include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
#include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
#include "Parallel/AlgorithmExecution.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/TMPL.hpp"
Expand Down Expand Up @@ -65,7 +66,8 @@ struct InitializeDomain {
using simple_tags = typename InitializeGeometry::return_tags;
using compute_tags = tmpl::list<>;
using const_global_cache_tags =
tmpl::list<domain::Tags::FunctionsOfTimeInitialize>;
tmpl::list<domain::Tags::FunctionsOfTimeInitialize,
elliptic::dg::Tags::Quadrature>;

template <typename DbTagsList, typename... InboxTags, typename Metavariables,
typename ActionList, typename ParallelComponent>
Expand Down
134 changes: 79 additions & 55 deletions src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,10 @@
#include "Elliptic/Systems/GetSourcesComputer.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleBoundaryData.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/SimpleMortarData.hpp"
#include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
Expand Down Expand Up @@ -53,9 +55,8 @@
* The DG discretization of equations in this first-order form amounts to
* projecting the equations on the set of basis functions that we also use to
* represent the fields on the computational grid. The currently implemented DG
* operator uses Lagrange interpolating polynomials w.r.t.
* Legendre-Gauss-Lobatto collocation points as basis functions. Support for
* Legendre-Gauss collocation points can be added if needed. Skipping all
* operator uses Lagrange interpolating polynomials w.r.t. Legendre-Gauss or
* Legendre-Gauss-Lobatto collocation points as basis functions. Skipping all
* further details here, the discretization results in a linear equation
* \f$A(u)=b\f$ over all grid points and primal variables. Solving the elliptic
* equations amounts to numerically inverting the DG operator \f$A\f$, typically
Expand Down Expand Up @@ -342,9 +343,10 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
#ifdef SPECTRE_DEBUG
for (size_t d = 0; d < Dim; ++d) {
ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
(mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto or
mesh.quadrature(d) == Spectral::Quadrature::Gauss),
"The elliptic DG operator is currently only implemented for "
"Legendre-Gauss-Lobatto grids. Found basis '"
"Legendre-Gauss(-Lobatto) grids. Found basis '"
<< mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
<< "' in dimension " << d << ".");
}
Expand Down Expand Up @@ -444,9 +446,8 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
const auto& face_normal = face_normals.at(direction);
const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
const auto& fluxes_args_on_face = fluxes_args_on_faces.at(direction);
const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
Variables<tmpl::list<PrimalFluxesVars..., AuxiliaryFluxesVars...>>
fluxes_on_face{};
Variables<tmpl::list<PrimalFluxesVars...>> primal_fluxes_on_face{};
Variables<tmpl::list<AuxiliaryFluxesVars...>> aux_fluxes_on_face{};
Variables<tmpl::list<::Tags::NormalDotFlux<AuxiliaryFields>...>>
n_dot_aux_fluxes{};
BoundaryData<tmpl::list<PrimalMortarVars...>,
Expand All @@ -457,16 +458,17 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
// data such as the element size will be set below.
boundary_data.field_data.initialize(face_num_points, 0.);
} else {
fluxes_on_face.initialize(face_num_points);
primal_fluxes_on_face.initialize(face_num_points);
aux_fluxes_on_face.initialize(face_num_points);
n_dot_aux_fluxes.initialize(face_num_points);
// Compute F_u(n.F_v)
fluxes_on_face.assign_subset(
data_on_slice(*auxiliary_fluxes, mesh.extents(),
direction.dimension(), slice_index));
::dg::project_contiguous_data_to_boundary(
make_not_null(&aux_fluxes_on_face), *auxiliary_fluxes, mesh,
direction);
EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
make_not_null(
&get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
face_normal, get<AuxiliaryFluxesVars>(aux_fluxes_on_face)));
std::apply(
[&boundary_data,
&n_dot_aux_fluxes](const auto&... expanded_fluxes_args_on_face) {
Expand All @@ -485,13 +487,13 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
// to the auxiliary variables. The reason is essentially that the
// internal penalty flux uses average(grad(u)) in place of average(v),
// i.e. the raw primal field derivatives without boundary corrections.
fluxes_on_face.assign_subset(
data_on_slice(*primal_fluxes, mesh.extents(), direction.dimension(),
slice_index));
::dg::project_contiguous_data_to_boundary(
make_not_null(&primal_fluxes_on_face), *primal_fluxes, mesh,
direction);
EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
boundary_data.field_data)),
face_normal, get<PrimalFluxesVars>(fluxes_on_face)));
face_normal, get<PrimalFluxesVars>(primal_fluxes_on_face)));
// Compute n.F_u(n.F_v) for jump term, re-using the memory buffer from
// above
EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
Expand Down Expand Up @@ -562,8 +564,8 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
if constexpr (AllDataIsZero) {
dirichlet_vars.initialize(face_num_points, 0.);
} else {
data_on_slice(make_not_null(&dirichlet_vars), primal_vars,
mesh.extents(), direction.dimension(), slice_index);
::dg::project_contiguous_data_to_boundary(
make_not_null(&dirichlet_vars), primal_vars, mesh, direction);
}
apply_boundary_condition(
direction, make_not_null(&get<PrimalVars>(dirichlet_vars))...,
Expand All @@ -573,21 +575,21 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
// The n.F_u (Neumann-type conditions) are done, but we have to compute
// the n.F_v (Dirichlet-type conditions) from the Dirichlet fields. We
// re-use the memory buffer from above.
fluxes_on_face.initialize(face_num_points);
aux_fluxes_on_face.initialize(face_num_points);
n_dot_aux_fluxes.initialize(face_num_points);
std::apply(
[&fluxes_on_face,
[&aux_fluxes_on_face,
&dirichlet_vars](const auto&... expanded_fluxes_args_on_face) {
FluxesComputer::apply(
make_not_null(&get<AuxiliaryFluxesVars>(fluxes_on_face))...,
expanded_fluxes_args_on_face...,
get<PrimalVars>(dirichlet_vars)...);
FluxesComputer::apply(make_not_null(&get<AuxiliaryFluxesVars>(
aux_fluxes_on_face))...,
expanded_fluxes_args_on_face...,
get<PrimalVars>(dirichlet_vars)...);
},
fluxes_args_on_face);
EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
make_not_null(
&get<::Tags::NormalDotFlux<AuxiliaryFields>>(n_dot_aux_fluxes)),
face_normal, get<AuxiliaryFluxesVars>(fluxes_on_face)));
face_normal, get<AuxiliaryFluxesVars>(aux_fluxes_on_face)));
std::apply(
[&boundary_data,
&n_dot_aux_fluxes](const auto&... expanded_fluxes_args_on_face) {
Expand Down Expand Up @@ -712,9 +714,10 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
#ifdef SPECTRE_DEBUG
for (size_t d = 0; d < Dim; ++d) {
ASSERT(mesh.basis(d) == Spectral::Basis::Legendre and
mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto,
(mesh.quadrature(d) == Spectral::Quadrature::GaussLobatto or
mesh.quadrature(d) == Spectral::Quadrature::Gauss),
"The elliptic DG operator is currently only implemented for "
"Legendre-Gauss-Lobatto grids. Found basis '"
"Legendre-Gauss(-Lobatto) grids. Found basis '"
<< mesh.basis(d) << "' and quadrature '" << mesh.quadrature(d)
<< "' in dimension " << d << ".");
}
Expand Down Expand Up @@ -772,7 +775,6 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
has_any_boundary_corrections = true;

const auto face_mesh = mesh.slice_away(direction.dimension());
const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
const auto& local_data = mortar_data.local_data(temporal_id);
const auto& remote_data = mortar_data.remote_data(temporal_id);
const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
Expand Down Expand Up @@ -805,22 +807,35 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
mortar_mesh, mortar_size, mortar_jacobians.at(mortar_id),
face_mesh, face_jacobians.at(direction))
: std::move(auxiliary_boundary_corrections_on_mortar);

// Lift the boundary correction to the volume, but still only provide the
// data only on the face because it is zero everywhere else. This is the
// "massless" lifting operation, i.e. it involves an inverse mass matrix.
// The mass matrix is diagonally approximated ("mass lumping") so it
// reduces to a division by quadrature weights.
::dg::lift_flux(make_not_null(&auxiliary_boundary_corrections),
mesh.extents(direction.dimension()),
face_normal_magnitude);
// The `dg::lift_flux` function contains an extra minus sign
// The `dg::lift_flux` function contains an extra minus sign (for
// evolution systems)
auxiliary_boundary_corrections *= -1.;

// Add the boundary corrections to the auxiliary variables
add_slice_to_data(make_not_null(&primal_fluxes_corrected),
auxiliary_boundary_corrections, mesh.extents(),
direction.dimension(), slice_index);
if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
// Lift the boundary correction to the volume, but still only provide
// the data only on the face because it is zero everywhere else. This is
// the "massless" lifting operation, i.e. it involves an inverse mass
// matrix. The mass matrix is diagonally approximated ("mass lumping")
// so it reduces to a division by quadrature weights.
::dg::lift_flux(make_not_null(&auxiliary_boundary_corrections),
mesh.extents(direction.dimension()),
face_normal_magnitude);

// Add the boundary corrections to the auxiliary variables
add_slice_to_data(make_not_null(&primal_fluxes_corrected),
auxiliary_boundary_corrections, mesh.extents(),
direction.dimension(),
index_to_slice_at(mesh.extents(), direction));
} else {
// We already have the `face_jacobians = det(J) * magnitude(n)` here, so
// just pass a constant 1 for `magnitude(n)`. This could be optimized to
// avoid allocating the vector of ones.
::dg::lift_boundary_terms_gauss_points(
make_not_null(&primal_fluxes_corrected), det_inv_jacobian, mesh,
direction, auxiliary_boundary_corrections,
Scalar<DataVector>{face_mesh.number_of_grid_points(), 1.},
face_jacobians.at(direction));
}
} // apply auxiliary boundary corrections on all mortars

// Compute the primal equation, i.e. the actual DG operator, by taking the
Expand Down Expand Up @@ -868,7 +883,6 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
}

const auto face_mesh = mesh.slice_away(direction.dimension());
const size_t slice_index = index_to_slice_at(mesh.extents(), direction);
const auto [local_data, remote_data] = mortar_data.extract();
const auto& face_normal_magnitude = face_normal_magnitudes.at(direction);
const auto& mortar_mesh =
Expand Down Expand Up @@ -927,11 +941,20 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
mortar_mesh, mortar_size, mortar_jacobians.at(mortar_id),
face_mesh, face_jacobians.at(direction))
: std::move(primal_boundary_corrections_on_mortar);
::dg::lift_flux(make_not_null(&primal_boundary_corrections),
mesh.extents(direction.dimension()),
face_normal_magnitude);
add_slice_to_data(operator_applied_to_vars, primal_boundary_corrections,
mesh.extents(), direction.dimension(), slice_index);
if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
::dg::lift_flux(make_not_null(&primal_boundary_corrections),
mesh.extents(direction.dimension()),
face_normal_magnitude);
add_slice_to_data(operator_applied_to_vars, primal_boundary_corrections,
mesh.extents(), direction.dimension(),
index_to_slice_at(mesh.extents(), direction));
} else {
::dg::lift_boundary_terms_gauss_points(
operator_applied_to_vars, det_inv_jacobian, mesh, direction,
primal_boundary_corrections,
Scalar<DataVector>{face_mesh.number_of_grid_points(), 1.},
face_jacobians.at(direction));
}
} // loop over all mortars

// Apply DG mass matrix
Expand All @@ -956,6 +979,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
const Scalar<DataVector>& det_inv_jacobian,
const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
const double penalty_parameter, const bool massive,
Expand Down Expand Up @@ -996,12 +1020,12 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
face_normal_magnitudes, all_mortar_meshes, all_mortar_sizes,
temporal_id, apply_boundary_condition, fluxes_args, sources_args,
fluxes_args_on_faces);
apply_operator<true>(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, {},
all_mortar_meshes, all_mortar_sizes, {},
penalty_parameter, massive, temporal_id, sources_args);
apply_operator<true>(
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, all_mortar_meshes, all_mortar_sizes, {},
penalty_parameter, massive, 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;
Expand Down
19 changes: 12 additions & 7 deletions src/Elliptic/DiscontinuousGalerkin/Initialization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include "Domain/Structure/DirectionMap.hpp"
#include "Domain/Structure/Element.hpp"
#include "Domain/Structure/IndexToSliceAt.hpp"
#include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
#include "NumericalAlgorithms/Spectral/Mesh.hpp"
Expand Down Expand Up @@ -52,9 +53,14 @@ void InitializeGeometry<Dim>::operator()(
const std::vector<std::array<size_t, Dim>>& initial_refinement,
const Domain<Dim>& domain,
const domain::FunctionsOfTimeMap& functions_of_time,
const Spectral::Quadrature quadrature,
const ElementId<Dim>& element_id) const {
// Mesh
const auto quadrature = Spectral::Quadrature::GaussLobatto;
ASSERT(quadrature == Spectral::Quadrature::GaussLobatto or
quadrature == Spectral::Quadrature::Gauss,
"The elliptic DG scheme supports Gauss and Gauss-Lobatto "
"grids, but the chosen quadrature is: "
<< quadrature);
*mesh = domain::Initialization::create_initial_mesh(initial_extents,
element_id, quadrature);
// Element
Expand Down Expand Up @@ -84,16 +90,15 @@ void deriv_unnormalized_face_normals_impl(
if (element.external_boundaries().empty()) {
return;
}
ASSERT(mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto,
"Slicing the Hessian to the boundary currently supports only "
"Gauss-Lobatto grids. Add support to "
"'elliptic::dg::InitializeFacesAndMortars'.");
// If the accuracy of this derivative is insufficient we could also compute
// it on a higher-order grid and then project it down.
// On Gauss grids we could compute the derivative on a Gauss-Lobatto grid and
// slice it.
const auto deriv_inv_jac =
partial_derivative(inv_jacobian, mesh, inv_jacobian);
for (const auto& direction : element.external_boundaries()) {
const auto deriv_inv_jac_on_face =
data_on_slice(deriv_inv_jac, mesh.extents(), direction.dimension(),
index_to_slice_at(mesh.extents(), direction));
::dg::project_tensor_to_boundary(deriv_inv_jac, mesh, direction);
auto& deriv_unnormalized_face_normal =
(*deriv_unnormalized_face_normals)[direction];
for (size_t i = 0; i < Dim; ++i) {
Expand Down
Loading
Loading