From 99e3b5eeab620f447c8346ec3e27982d6ade9007 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 23 Aug 2024 14:04:11 -0700 Subject: [PATCH] Add `logical_divergence` function --- .../LinearOperators/Divergence.cpp | 46 +++++++++++++++++++ .../LinearOperators/Divergence.hpp | 23 ++++++++++ .../LinearOperators/Divergence.tpp | 33 +++++++++++++ .../LinearOperators/Test_Divergence.cpp | 13 +++++- 4 files changed, 113 insertions(+), 2 deletions(-) diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.cpp b/src/NumericalAlgorithms/LinearOperators/Divergence.cpp index b5cba1a54c47..09185380fcfa 100644 --- a/src/NumericalAlgorithms/LinearOperators/Divergence.cpp +++ b/src/NumericalAlgorithms/LinearOperators/Divergence.cpp @@ -6,6 +6,7 @@ #include #include +#include "DataStructures/ApplyMatrices.hpp" #include "DataStructures/ComplexDataVector.hpp" #include "DataStructures/DataBox/Tag.hpp" #include "DataStructures/DataVector.hpp" @@ -61,6 +62,26 @@ void divergence(const gsl::not_null*> div_input, } } +template +void logical_divergence(const gsl::not_null div_flux, + const FluxTensor& flux, const Mesh& mesh) { + // Note: This function hasn't been optimized much at all. Feel free to + // optimize if needed! + static const Matrix identity_matrix{}; + for (size_t d = 0; d < Dim; ++d) { + auto matrices = make_array(std::cref(identity_matrix)); + gsl::at(matrices, d) = + Spectral::differentiation_matrix(mesh.slice_through(d)); + for (size_t storage_index = 0; storage_index < div_flux->size(); + ++storage_index) { + const auto div_flux_index = div_flux->get_tensor_index(storage_index); + const auto flux_index = prepend(div_flux_index, d); + div_flux->get(div_flux_index) += + apply_matrices(matrices, flux.get(flux_index), mesh.extents()); + } + } +} + #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) #define DIM(data) BOOST_PP_TUPLE_ELEM(1, data) #define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data) @@ -85,3 +106,28 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector, ComplexDataVector), (1, 2, 3), #undef DIM #undef FRAME #undef INSTANTIATE + +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) +#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data) +#define TENSOR(data) BOOST_PP_TUPLE_ELEM(2, data) + +#define INSTANTIATION_SCALAR(r, data) \ + template void logical_divergence( \ + const gsl::not_null*> div_flux, \ + const tnsr::I& flux, \ + const Mesh& mesh); +#define INSTANTIATION_TENSOR(r, data) \ + template void logical_divergence( \ + const gsl::not_null* > div_flux, \ + const TensorMetafunctions::prepend_spatial_index< \ + tnsr::TENSOR(data) < DTYPE(data), DIM(data), Frame::Inertial>, \ + DIM(data), UpLo::Up, Frame::ElementLogical > &flux, \ + const Mesh& mesh); + +GENERATE_INSTANTIATIONS(INSTANTIATION_SCALAR, (DataVector, ComplexDataVector), + (1, 2, 3)) +GENERATE_INSTANTIATIONS(INSTANTIATION_TENSOR, (DataVector, ComplexDataVector), + (1, 2, 3), (i, I)) + +#undef INSTANTIATION diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp index 88a6f759cd1e..ef10009e2764 100644 --- a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp +++ b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp @@ -89,6 +89,29 @@ void divergence(gsl::not_null*> div_input, DerivativeFrame>& inverse_jacobian); /// @} +/// @{ +/*! + * \brief Compute the divergence of fluxes in logical coordinates + * + * Applies the logical differentiation matrix to the fluxes in each dimension + * and sums over dimensions. + * + * \see divergence + */ +template +void logical_divergence(gsl::not_null div_flux, + const FluxTensor& flux, const Mesh& mesh); + +template +auto logical_divergence(const Variables& flux, const Mesh& mesh) + -> Variables>; + +template +void logical_divergence( + gsl::not_null>*> div_flux, + const Variables>& flux, const Mesh& mesh); +/// @} + namespace Tags { /*! * \ingroup DataBoxTagsGroup diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp index 295a861996a7..14f94ee87f42 100644 --- a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp +++ b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp @@ -96,3 +96,36 @@ void divergence( }; EXPAND_PACK_LEFT_TO_RIGHT(apply_div(FluxTags{}, DivTags{})); } + +template +Variables> logical_divergence( + const Variables& flux, const Mesh& mesh) { + Variables> div_flux( + flux.number_of_grid_points()); + logical_divergence(make_not_null(&div_flux), flux, mesh); + return div_flux; +} + +template +void logical_divergence( + const gsl::not_null>*> div_flux, + const Variables>& flux, const Mesh& mesh) { + static_assert( + (... and + std::is_same_v, + SpatialIndex>), + "The first index of each flux must be an upper spatial index in " + "element-logical coordinates."); + static_assert( + (... and + std::is_same_v< + typename ResultTags::type, + TensorMetafunctions::remove_first_index>), + "The result tensors must have the same type as the flux tensors with " + "their first index removed."); + div_flux->initialize(mesh.number_of_grid_points(), 0.); + // Note: This function hasn't been optimized much at all. Feel free to + // optimize if needed! + EXPAND_PACK_LEFT_TO_RIGHT(logical_divergence( + make_not_null(&get(*div_flux)), get(flux), mesh)); +} diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp index f0bb4165a833..f59f932d42e1 100644 --- a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp +++ b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp @@ -17,8 +17,10 @@ #include "DataStructures/DataBox/Prefixes.hpp" #include "DataStructures/DataBox/Tag.hpp" #include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/Determinant.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" +#include "DataStructures/Variables/FrameTransform.hpp" #include "DataStructures/VariablesTag.hpp" #include "Domain/CoordinateMaps/Affine.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" @@ -143,6 +145,7 @@ void test_divergence_impl( const auto xi = logical_coordinates(mesh); const auto x = coordinate_map(xi); const auto inv_jacobian = coordinate_map.inv_jacobian(xi); + const auto det_jacobian = determinant(coordinate_map.jacobian(xi)); MathFunctions::TensorProduct f(1.0, std::move(functions)); using flux_tags = two_fluxes; Variables fluxes(num_grid_points); @@ -154,7 +157,7 @@ void test_divergence_impl( using DivFluxTag = Tags::div; get(expected_div_fluxes) = FluxTag::divergence_of_flux(f, x); }); - const auto div_fluxes = divergence(fluxes, mesh, inv_jacobian); + auto div_fluxes = divergence(fluxes, mesh, inv_jacobian); CHECK(div_fluxes.size() == expected_div_fluxes.size()); CHECK(Dim * div_fluxes.size() == fluxes.size()); Approx local_approx = Approx::custom().epsilon(1.e-11).scale(1.); @@ -166,6 +169,12 @@ void test_divergence_impl( const auto& expected = get>>(div_fluxes); CHECK_ITERABLE_CUSTOM_APPROX(expected, div_vector, local_approx); + + // Test logical divergence + auto logical_fluxes = + transform::first_index_to_different_frame(fluxes, inv_jacobian); + logical_divergence(make_not_null(&div_fluxes), logical_fluxes, mesh); + CHECK_VARIABLES_CUSTOM_APPROX(div_fluxes, expected_div_fluxes, local_approx); } template @@ -314,7 +323,7 @@ void test_divergence_compute() { } } // namespace -// [[Timeout, 10]] +// [[Timeout, 20]] SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.Divergence", "[NumericalAlgorithms][LinearOperators][Unit]") { test_divergence();