From 2d30004964345b0e2a0b3b314bcc0a2300bd09e2 Mon Sep 17 00:00:00 2001 From: Michael Murphy Date: Mon, 21 Jul 2025 11:50:40 -0400 Subject: [PATCH 1/2] Allow Cartoon bases to have one grid point in a mesh --- .../Spectral/GetSpectralQuantityForMesh.hpp | 17 +++++++ .../Spectral/LogicalCoordinates.cpp | 9 ++++ .../Spectral/LogicalCoordinates.hpp | 3 +- .../Spectral/Test_LogicalCoordinates.cpp | 48 +++++++++++++++++++ 4 files changed, 76 insertions(+), 1 deletion(-) diff --git a/src/NumericalAlgorithms/Spectral/GetSpectralQuantityForMesh.hpp b/src/NumericalAlgorithms/Spectral/GetSpectralQuantityForMesh.hpp index ff1f87e05550..e70ff30a972c 100644 --- a/src/NumericalAlgorithms/Spectral/GetSpectralQuantityForMesh.hpp +++ b/src/NumericalAlgorithms/Spectral/GetSpectralQuantityForMesh.hpp @@ -48,6 +48,23 @@ decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) { default: ERROR("Missing quadrature case for spectral quantity"); } + case Basis::Cartoon: + switch (mesh.quadrature(0)) { + case Quadrature::AxialSymmetry: + return f( + std::integral_constant{}, + std::integral_constant{}, + num_points); + case Quadrature::SphericalSymmetry: + return f(std::integral_constant{}, + std::integral_constant{}, + num_points); + default: + ERROR( + "Only Axial and Spherical Symmetry quadratures are allowed for " + "a Cartoon basis."); + } case Basis::FiniteDifference: switch (mesh.quadrature(0)) { case Quadrature::CellCentered: diff --git a/src/NumericalAlgorithms/Spectral/LogicalCoordinates.cpp b/src/NumericalAlgorithms/Spectral/LogicalCoordinates.cpp index 0cf6e7dbb9aa..7c1a234c8124 100644 --- a/src/NumericalAlgorithms/Spectral/LogicalCoordinates.cpp +++ b/src/NumericalAlgorithms/Spectral/LogicalCoordinates.cpp @@ -65,6 +65,15 @@ void logical_coordinates( } break; } + case Spectral::Basis::Cartoon: { + if (mesh.extents(d) != 1) { + ERROR("Only 1 grid point is allowed in a Cartoon basis."); + } + for (IndexIterator index(mesh.extents()); index; ++index) { + logical_coords->get(d)[index.collapsed_index()] = 0.0; + } + break; + } // NOLINTNEXTLINE(bugprone-branch-clone) case Spectral::Basis::Chebyshev: [[fallthrough]]; diff --git a/src/NumericalAlgorithms/Spectral/LogicalCoordinates.hpp b/src/NumericalAlgorithms/Spectral/LogicalCoordinates.hpp index c8950ce81b91..056fd95ef537 100644 --- a/src/NumericalAlgorithms/Spectral/LogicalCoordinates.hpp +++ b/src/NumericalAlgorithms/Spectral/LogicalCoordinates.hpp @@ -57,7 +57,8 @@ struct not_null; * corresponding to the Legendre Gauss points of \f$\cos \theta\f$ (and thus * have no points at the poles), while the azimuth has Equiangular quadrature * which are distributed uniformly including the left endpoint, but not the - * right. + * right. For the Cartoon basis, only one logical coordinate is allowed and it + * is never used for computations. * * \example * \snippet Test_LogicalCoordinates.cpp logical_coordinates_example diff --git a/tests/Unit/NumericalAlgorithms/Spectral/Test_LogicalCoordinates.cpp b/tests/Unit/NumericalAlgorithms/Spectral/Test_LogicalCoordinates.cpp index d5716d2d129f..4678c079bec3 100644 --- a/tests/Unit/NumericalAlgorithms/Spectral/Test_LogicalCoordinates.cpp +++ b/tests/Unit/NumericalAlgorithms/Spectral/Test_LogicalCoordinates.cpp @@ -13,6 +13,7 @@ #include "Domain/CoordinateMaps/Affine.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" #include "Domain/CoordinateMaps/CoordinateMap.tpp" +#include "Domain/CoordinateMaps/Identity.hpp" #include "Domain/CoordinateMaps/ProductMaps.hpp" #include "Domain/CoordinateMaps/ProductMaps.tpp" #include "Domain/Structure/Direction.hpp" @@ -106,6 +107,53 @@ SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Spectral.LogicalCoordinates", CHECK(x_3d[2][0] == -32.0); CHECK(x_3d[2][15] == 74.0); + const Mesh<3> mesh_spherical{ + {3, 1, 1}, + {Spectral::Basis::Legendre, Spectral::Basis::Cartoon, + Spectral::Basis::Cartoon}, + {Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::SphericalSymmetry, + Spectral::Quadrature::SphericalSymmetry}}; + using Affine = domain::CoordinateMaps::Affine; + using Identity = domain::CoordinateMaps::Identity<1>; + const Identity identity_cartoon_map; + const Affine affine_x_map(-1.0, 1.0, -1.0, 1.0); + const domain::CoordinateMap< + Frame::ElementLogical, Frame::Grid, + domain::CoordinateMaps::ProductOf3Maps> + map_spherical{{affine_x_map, identity_cartoon_map, identity_cartoon_map}}; + const auto x_spherical = map_spherical(logical_coordinates(mesh_spherical)); + CHECK(x_spherical[1][0] == 0.0); + CHECK(x_spherical[1][1] == 0.0); + CHECK(x_spherical[1][2] == 0.0); + CHECK(x_spherical[2][0] == 0.0); + CHECK(x_spherical[2][1] == 0.0); + CHECK(x_spherical[2][2] == 0.0); + + const Mesh<3> mesh_axial{ + {2, 3, 1}, + {Spectral::Basis::Legendre, Spectral::Basis::Legendre, + Spectral::Basis::Cartoon}, + {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::AxialSymmetry}}; + const domain::CoordinateMap< + Frame::ElementLogical, Frame::Grid, + domain::CoordinateMaps::ProductOf3Maps> + map_axial{{affine_x_map, affine_x_map, identity_cartoon_map}}; + const auto x_axial = map_axial(logical_coordinates(mesh_axial)); + CHECK(x_axial[2][0] == 0.0); + CHECK(x_axial[2][1] == 0.0); + CHECK(x_axial[2][2] == 0.0); + CHECK(x_axial[2][3] == 0.0); + CHECK(x_axial[2][4] == 0.0); + CHECK(x_axial[2][5] == 0.0); + + const Mesh<1> mesh_cartoon_error{2, Spectral::Basis::Cartoon, + Spectral::Quadrature::AxialSymmetry}; + CHECK_THROWS_WITH((logical_coordinates(mesh_cartoon_error)), + Catch::Matchers::ContainsSubstring( + "Only 1 grid point is allowed in a Cartoon basis.")); + TestHelpers::db::test_compute_tag>( "ElementLogicalCoordinates"); TestHelpers::db::test_compute_tag>( From c6556d3edf556cb417d89d665e4c141953617c3f Mon Sep 17 00:00:00 2001 From: Michael Murphy Date: Mon, 21 Jul 2025 11:59:37 -0400 Subject: [PATCH 2/2] Implement DG `cartoon_partial_derivatives` --- src/DataStructures/Tags/TempTensor.hpp | 28 ++ .../LinearOperators/PartialDerivatives.hpp | 73 +++ .../LinearOperators/PartialDerivatives.tpp | 429 +++++++++++++++++- .../Test_PartialDerivatives.cpp | 363 +++++++++++++++ 4 files changed, 891 insertions(+), 2 deletions(-) diff --git a/src/DataStructures/Tags/TempTensor.hpp b/src/DataStructures/Tags/TempTensor.hpp index 0b66beb67714..1003d43ff677 100644 --- a/src/DataStructures/Tags/TempTensor.hpp +++ b/src/DataStructures/Tags/TempTensor.hpp @@ -16,6 +16,9 @@ struct Inertial; } // namespace Frame /// \endcond +/*! + * \brief Contains objects related to Tags + */ namespace Tags { template struct TempTensor : db::SimpleTag { @@ -25,6 +28,31 @@ struct TempTensor : db::SimpleTag { } }; +/*! + * \brief A struct that constructs a `TempTensor` given a label \p N and type + * \p TensorType. Call with `make_temp_tensor::apply` + * + * \see convert_to_temp_tensors + */ +template +struct make_temp_tensor { + template + struct apply { + using type = ::Tags::TempTensor; + }; +}; + +/*! + * \brief Takes a `tmpl::list` of `Tensor` types and a label \p N which is + * applied to all `TempTensor`s in the resulting list. + * + * \see make_temp_tensor + */ +template +using convert_to_temp_tensors = + tmpl::transform::template apply>; + /// @{ /// \ingroup PeoGroup /// Variables Tags for temporary tensors inside a function. diff --git a/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp b/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp index 68dc3bd16d78..194da1b9501d 100644 --- a/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp +++ b/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp @@ -111,6 +111,79 @@ auto logical_partial_derivative(const Tensor& u, Frame::ElementLogical>; /// @} +/// @{ +/*! + * \ingroup NumericalAlgorithmsGroup + * \brief Compute the partial derivatives of each variable, using the Cartoon + * method for directions not in the computational domain. + * + * For a spacetime with some Killing vector \f$\xi^a\f$, not only does the + * metric respect the isometry via \f$\mathcal{L}_\xi g_{ab} = 0\f$, but the + * fields must as well. + * + * In the case of axial symmetry about the \f$y\f$-axis, the Killing vector is + * \f$\xi^a = (0, -z, 0, x)\f$, so we have that for a vector \f$V^a\f$, + * + * \f{align*}{ + * \mathcal{L}_\xi V^a &= \xi^b \partial_b V^a - V^b \partial_b \xi^a + * = -z \partial_x V^a + x \partial_z V^a - V^b \partial_b \xi^a = 0. + * \f} + * + * Choosing the computational domain to be the \f$z=0\f$ plane, get the result: + * + * \f{align*}{ + * \partial_z V^a &= \frac{V^b \partial_b \xi^a}{x}. + * \f} + * + * In the case that $x=0$ is in the computational domain, L'Hôpital's + * rule is used. For a general tensor, + * \f$T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m}\f$ we have + * + * \f{align*}{ + * \partial_z T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m} &= + * \left\{ + * \begin{array}{ll} + * \partial_x(\sum_k + * T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k} + * \xi^{a_k} & \\ + * \;\;\;\;\;-\sum_k + * T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k} + * \xi^{c_k} + * ) & \mathrm{if} \; x=0 \\ + * (\sum_k + * T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k} + * \xi^{a_k} & \\ + * \;-\sum_k + * T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k} + * \xi^{c_k})\frac{1}{x} & \mathrm{otherwise} + * \end{array}\right. + * \f} + * + * In the case of spherical symmetry, another Killing vector + * \f$\xi'^a = (0, -y, x, 0) \f$ is used as well, with the above equation + * easily generalizing. In that case, the computational domain is only the + * \f$x\f$-axis. + * + * This function computes the partial derivatives of _all_ variables in + * `VariableTags`. The additional parameter `inertial_coords` is used for + * division by the \f$x\f$ coordinates. If \f$x=0\f$ is included in the domain, + * it is assumed to be present only at the first index and is handled by + * L'Hôpital's rule. + * + * The mesh is required to have the Cartoon basis in the last and potentially + * second-to-last coordinates and the inverse Jacobian is accordingly used only + * in the first and potentially second dimensions. + */ +template = nullptr> +void cartoon_partial_derivatives( + gsl::not_null*> du, const Variables& u, + const Mesh& mesh, + const InverseJacobian& inverse_jacobian_3d, + const tnsr::I& inertial_coords); +/// @} + /// @{ /// \ingroup NumericalAlgorithmsGroup /// \brief Compute the partial derivatives of each variable with respect to diff --git a/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.tpp b/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.tpp index 1ebc603bc7f3..5a4bf9abdacd 100644 --- a/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.tpp +++ b/src/NumericalAlgorithms/LinearOperators/PartialDerivatives.tpp @@ -9,15 +9,21 @@ #include "DataStructures/DataBox/Prefixes.hpp" #include "DataStructures/DataVector.hpp" #include "DataStructures/Matrix.hpp" +#include "DataStructures/Tags/TempTensor.hpp" +#include "DataStructures/Tensor/IndexType.hpp" +#include "DataStructures/Tensor/Metafunctions.hpp" #include "DataStructures/Transpose.hpp" #include "DataStructures/Variables.hpp" +#include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/DifferentiationMatrix.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Quadrature.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" #include "NumericalAlgorithms/SphericalHarmonics/SpherepackCache.hpp" #include "Utilities/Algorithm.hpp" #include "Utilities/Blas.hpp" #include "Utilities/ContainerHelpers.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeArray.hpp" #include "Utilities/MemoryHelpers.hpp" @@ -190,6 +196,425 @@ void partial_derivatives( inverse_jacobian); } +constexpr tnsr::iab Killing_vector_derivatives() { + // holds derivative of both (0,-y,x,0) for \partial_y and (0,-z,0,x) for + // \partial_z + tnsr::iab da_killing_vectors{}; + for (size_t i = 0; i < 3; ++i) { + for (size_t KV_deriv_index = 0; KV_deriv_index < 4; ++KV_deriv_index) { + for (size_t Killing_index = 0; Killing_index < 4; ++Killing_index) { + if (i == 1 and KV_deriv_index == 2 and Killing_index == 1) { + // y derivative + da_killing_vectors.get(i, KV_deriv_index, Killing_index) = -1.0; + } else if (i == 1 and KV_deriv_index == 1 and Killing_index == 2) { + // y derivative + da_killing_vectors.get(i, KV_deriv_index, Killing_index) = 1.0; + } else if (i == 2 and KV_deriv_index == 3 and Killing_index == 1) { + // z derivative + da_killing_vectors.get(i, KV_deriv_index, Killing_index) = -1.0; + } else if (i == 2 and KV_deriv_index == 1 and Killing_index == 3) { + // z derivative + da_killing_vectors.get(i, KV_deriv_index, Killing_index) = 1.0; + } else { + da_killing_vectors.get(i, KV_deriv_index, Killing_index) = 0.0; + } + } + } + } + return da_killing_vectors; +} + +template +void cartoon_contraction( + const gsl::not_null*> result_tensor, + const Tensor& tensor, const size_t deriv_index, + const Spectral::Quadrature quad_type) { + // This function does the contractions in Eqn. (218) of the SXS book's + // Numerical Method's chapter, specifically the thing on the RHS that you take + // the derivative of or divide by x. + // Note that the result_tensor is written over and will have + // (*result_tensor)[0].size() == tensor[0].size() + auto& contract_tensor = *result_tensor; + ASSERT(quad_type == Spectral::Quadrature::AxialSymmetry or + quad_type == Spectral::Quadrature::SphericalSymmetry, + "Must pass a valid Cartoon quadrature"); + + constexpr tnsr::iab da_killing_vectors = + Killing_vector_derivatives(); + const std::array::rank()> valences = + tensor.index_valences(); + const auto type_index = tensor.index_types(); + + const size_t valence_size = valences.size(); + auto input_tensor_array = make_array(0_st); + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array Killing_indices; + Killing_indices[0] = deriv_index; + + for (size_t i = 0; i < tensor.size(); ++i) { + const auto tensor_index = tensor.get_tensor_index(i); + contract_tensor.get(tensor_index) = 0.0 * tensor.get(tensor_index); + for (size_t rank = 0; rank < valence_size; ++rank) { + const double sign = (valences[rank] == UpLo::Up) ? 1.0 : -1.0; + const size_t max_dummy = type_index[rank] == IndexType::Spacetime ? 4 : 3; + const size_t shift_index = + type_index[rank] == IndexType::Spacetime ? 0 : 1; + + // loop over dimension of rank + for (size_t dummy = 0; dummy < max_dummy; ++dummy) { + if (valences[rank] == UpLo::Lo) { + // covariant index + Killing_indices[1] = tensor_index[rank] + shift_index; + Killing_indices[2] = dummy + shift_index; + } else { + // contravariant + Killing_indices[1] = dummy + shift_index; + Killing_indices[2] = tensor_index[rank] + shift_index; + } + if (da_killing_vectors.get(Killing_indices) != 0) { + for (size_t j = 0; j < tensor_index.size(); ++j) { + input_tensor_array[j] = tensor_index[j]; + } + input_tensor_array[rank] = dummy; + + contract_tensor.get(tensor_index) += + sign * tensor.get(input_tensor_array) * + da_killing_vectors.get(Killing_indices); + } + } + } + } +} + +template = nullptr> +void cartoon_derivative( + TensorMetafunctions::prepend_spatial_index& d_tensor, + const InputTensorDataType& tensor, const DataVector& safe_x_coords, + const Spectral::Quadrature quad_type) { + // This function assumes `safe_x_coords` does not contain zero (safe in the + // sense that having x in the denominator of an expression will not cause + // an FPE). + // Therefore, if the domain actually does have zero, the coordinates must be + // made safe by replacing all occurances of zero with a non-zero value, and + // the output from this function is then not valid at those coords, one must + // then manually use L'Hopital's rule + ASSERT(quad_type == Spectral::Quadrature::AxialSymmetry or + quad_type == Spectral::Quadrature::SphericalSymmetry, + "Must pass a valid Cartoon quadrature"); + ASSERT( + not equal_within_roundoff(0.0, safe_x_coords[0], + std::numeric_limits::epsilon() * 100.0, + max(safe_x_coords)), + "Invalid coordinate: safe_x_coords[0] is too close to zero (" + << safe_x_coords[0] + << "). Division by this value may cause an FPE if the value is zero. " + "Maximum coordinate: " + << max(safe_x_coords)); + + const bool spherical_sym = + quad_type == Spectral::Quadrature::SphericalSymmetry; + const size_t start_deriv_index = spherical_sym ? 1 : 2; + + InputTensorDataType tensor_holder; + + for (size_t deriv_index = start_deriv_index; deriv_index < 3; ++deriv_index) { + cartoon_contraction(make_not_null(&tensor_holder), tensor, deriv_index, + quad_type); + for (size_t component_index = 0; component_index < tensor_holder.size(); + ++component_index) { + const auto input_index = tensor_holder.get_tensor_index(component_index); + const auto output_index = prepend(input_index, deriv_index); + d_tensor.get(output_index) = + tensor_holder.get(input_index) / safe_x_coords; + } + } +} + +template +void cartoon_partial_derivatives_apply( + const gsl::not_null*> du, + const Variables& u, const Mesh& mesh, + const InverseJacobian& inverse_jacobian_3d, + const tnsr::I& inertial_coords) { + using DerivativeTags = + tmpl::front>>; + static_assert(Dim == 3); + static_assert( + std::is_same_v< + tmpl::transform>, + tmpl::transform, DerivativeFrame>, + tmpl::bind>>); + ASSERT((Comp_dim == 2 and + mesh.quadrature(2) == Spectral::Quadrature::AxialSymmetry) or + (Comp_dim == 1 and + (mesh.quadrature(1) == Spectral::Quadrature::SphericalSymmetry and + mesh.quadrature(2) == Spectral::Quadrature::SphericalSymmetry)), + "Invalid Quadrature combinations: axial symmetry requires 2 " + "non-Cartoon dimensions, spherical symmetry requires 1 non-Cartoon " + "dimension. Got: " + << mesh.quadrature()); + + using ValueType = typename Variables::value_type; + auto& partial_derivatives_of_u = *du; + + if (UNLIKELY(partial_derivatives_of_u.number_of_grid_points() != + mesh.number_of_grid_points())) { + partial_derivatives_of_u.initialize(mesh.number_of_grid_points()); + } + + const size_t vars_size = + u.number_of_grid_points() * + Variables::number_of_independent_components; + const auto logical_derivs_data = + cpp20::make_unique_for_overwrite( + (Comp_dim > 1 ? (Comp_dim + 1) : Comp_dim) * vars_size); + + std::array logical_derivs{}; + for (size_t i = 0; i < Comp_dim; ++i) { + gsl::at(logical_derivs, i) = &(logical_derivs_data[i * vars_size]); + } + Variables temp{}; + if constexpr (Comp_dim > 1) { + temp.set_data_ref(&logical_derivs_data[Comp_dim * vars_size], vars_size); + } + + if constexpr (Comp_dim == 1) { + partial_derivatives_detail:: + LogicalImpl::apply( + make_not_null(&logical_derivs), &partial_derivatives_of_u, &temp, u, + mesh.slice_through(0)); + } else { + partial_derivatives_detail:: + LogicalImpl::apply( + make_not_null(&logical_derivs), &partial_derivatives_of_u, &temp, u, + mesh.slice_through(0, 1)); + } + + std::array const_logical_derivs{}; + for (size_t i = 0; i < Comp_dim; ++i) { + gsl::at(const_logical_derivs, i) = gsl::at(logical_derivs, i); + } + + const Spectral::Quadrature quad_type = mesh.quadrature(2); + const auto numerical_deriv_in_this_direction = [](const size_t deriv_num) { + if constexpr (Comp_dim == 2) { + return deriv_num == 0 or deriv_num == 1; + } else { + return deriv_num == 0; + } + }; + // only the 1st (possibly 2nd) dimensions of the Jacobian are relevant here + const InverseJacobian + inverse_jacobian{inverse_jacobian_3d.begin()->size()}; + for (size_t i = 0; i < Comp_dim; ++i) { + for (size_t j = 0; j < Comp_dim; ++j) { + make_const_view(make_not_null(&inverse_jacobian.get(i, j)), + inverse_jacobian_3d.get(i, j), 0, + inverse_jacobian_3d.begin()->size()); + } + } + // pdu points to du + double* pdu = du->data(); + const size_t num_grid_points = du->number_of_grid_points(); + DataVector lhs{}; + DataVector logical_du{}; + DataVector safe_x_coords = get<0>(inertial_coords); + const bool element_contains_zero_of_symmetry_axis = equal_within_roundoff( + 0.0, safe_x_coords[0], std::numeric_limits::epsilon() * 100.0, + max(safe_x_coords)); + // Having x=0 requires both not dividing by zero (done here by setting to + // arbitrary non-zero value) and remembering to then do L'Hopital + if (element_contains_zero_of_symmetry_axis) { + safe_x_coords[0] = 1.0; + if constexpr (Comp_dim == 2) { + const size_t x_extents = mesh.extents(0); + const size_t y_extents = mesh.extents(1); + for (size_t i = 1; i < y_extents; ++i) { + ASSERT( + equal_within_roundoff( + 0.0, safe_x_coords[i * x_extents], + std::numeric_limits::epsilon() * 100.0, + max(safe_x_coords)), + "The passed inertial coordinates do not follow the required format " + "of a rectangular mesh with x=0 being located at the first index " + "of each constant-y subdomain of the x DataVector. x=" + << safe_x_coords); + safe_x_coords[i * x_extents] = 1.0; + } + } + } + // If doing L'Hopital's rule, we need to store temporary forms of the data + // in two different tensors (hence the 0 & 1 labels) for each tensor type + // They only have to hold the x=0 positions, so of size y_extents + using VariableTypes = tmpl::remove_duplicates< + tmpl::transform>>; + + using TempTags0 = ::Tags::convert_to_temp_tensors; + using TempTags1 = ::Tags::convert_to_temp_tensors; + using VarTags = tmpl::append; + + Variables temp_vars{}; + if (element_contains_zero_of_symmetry_axis) { + const size_t y_extents = mesh.extents(1); + temp_vars.initialize(y_extents); + } + + std::array, Comp_dim> indices{}; + for (size_t deriv_index = 0; deriv_index < Comp_dim; ++deriv_index) { + for (size_t d = 0; d < Comp_dim; ++d) { + gsl::at(gsl::at(indices, d), deriv_index) = + InverseJacobian::get_storage_index(d, deriv_index); + } + } + + // the non-cartoon derivatives need to know where they are within `u` to + // contract with the inverse Jacobian to compute the inertial derivatives + size_t global_component_index = 0; + tmpl::for_each( + [&u, &du, &lhs, &pdu, &num_grid_points, &logical_du, + &const_logical_derivs, &inverse_jacobian, &indices, &temp_vars, + &global_component_index, &safe_x_coords, + &numerical_deriv_in_this_direction, &quad_type, &mesh, + &element_contains_zero_of_symmetry_axis]( + tmpl::type_ /*meta*/) { + auto& tensor = get(u); + + TensorMetafunctions::prepend_spatial_index< + typename TensorTag::type, Dim, UpLo::Lo, Frame::Inertial> + cart_deriv_tensor; + cartoon_derivative( + cart_deriv_tensor, tensor, safe_x_coords, quad_type); + + for (size_t component_index = 0; + component_index < TensorTag::type::size(); + ++component_index, ++global_component_index) { + for (size_t deriv_index = 0; deriv_index < Dim; ++deriv_index) { + // lhs points to pdu (shifts by num grid points below) + lhs.set_data_ref(pdu, num_grid_points); + + if (numerical_deriv_in_this_direction(deriv_index)) { + // do normal derivative for x (and possibly y) derivative + // clang-tidy: const cast is fine since we won't modify the data + // and we need it to easily hook into the expression templates. + logical_du.set_data_ref( + const_cast( // NOLINT + gsl::at(const_logical_derivs, 0)) + // NOLINT + global_component_index * num_grid_points, + num_grid_points); + lhs = (*(inverse_jacobian.begin() + + gsl::at(indices[0], deriv_index))) * + logical_du; + if constexpr (Comp_dim == 2) { + // clang-tidy: const cast is fine since we won't modify the data + // and we need it to easily hook into the expression templates. + logical_du.set_data_ref( + const_cast( // NOLINT + gsl::at(const_logical_derivs, 1)) + // NOLINT + global_component_index * num_grid_points, + num_grid_points); + lhs += (*(inverse_jacobian.begin() + + gsl::at(gsl::at(indices, 1), deriv_index))) * + logical_du; + } + } else { + // do cartoon dervative + const auto input_tensor_index = + tensor.get_tensor_index(component_index); + const auto output_tensor_index = + prepend(input_tensor_index, size_t{deriv_index}); + lhs = cart_deriv_tensor.get(output_tensor_index); + } + // clang-tidy: no pointer arithmetic + pdu += num_grid_points; // NOLINT + } + } + if (element_contains_zero_of_symmetry_axis) { + // need to use L'Hopital + using dtensor_type = + tmpl::at>; + auto& dtensor = get(*du); + + // if spherical symmetry,the zero is only at the first index + // if axial symmetry, the zero is repeated every x_extents times + const size_t x_extents = mesh.extents(0); + // for spherical symmetry, y_extents = 1, so we could use a double + // type rather than a size 1 DataVector, but it requires a lot of code + // reuse + const size_t y_extents = mesh.extents(1); + + auto& dx_tensor = + get<::Tags::TempTensor<0, typename TensorTag::type>>(temp_vars); + auto& contracted_dx_tensor = + get<::Tags::TempTensor<1, typename TensorTag::type>>(temp_vars); + + for (size_t component_index = 0; component_index < dx_tensor.size(); + ++component_index) { + const auto output_index = + dx_tensor.get_tensor_index(component_index); + // the 0 is because we want the x derivative + const size_t input_index = + dtensor.get_storage_index(prepend(output_index, size_t{0})); + for (size_t i = 0; i < y_extents; ++i) { + dx_tensor[component_index][i] = + dtensor[input_index][i * x_extents]; + } + } + + const size_t start_index = + quad_type == Spectral::Quadrature::SphericalSymmetry ? 1 : 2; + for (size_t deriv_index = start_index; deriv_index < 3; + ++deriv_index) { + cartoon_contraction(make_not_null(&contracted_dx_tensor), dx_tensor, + deriv_index, quad_type); + for (size_t component_index = 0; + component_index < contracted_dx_tensor.size(); + ++component_index) { + const auto input_index = + contracted_dx_tensor.get_tensor_index(component_index); + const size_t output_index = + dtensor.get_storage_index(prepend(input_index, deriv_index)); + for (size_t i = 0; i < y_extents; ++i) { + dtensor[output_index][i * x_extents] = + contracted_dx_tensor[component_index][i]; + } + } + } + } + }); +} + +template > +void cartoon_partial_derivatives( + const gsl::not_null*> du, + const Variables& u, const Mesh& mesh, + const InverseJacobian& inverse_jacobian_3d, + const tnsr::I& inertial_coords) { + if (mesh.basis(0) != Spectral::Basis::Cartoon and + mesh.basis(2) == Spectral::Basis::Cartoon) { + // computational dimension needs to be a constexpr + if (mesh.basis(1) == Spectral::Basis::Cartoon) { + cartoon_partial_derivatives_apply<1, ResultTags, VariableTags, Dim, + DerivativeFrame>( + du, u, mesh, inverse_jacobian_3d, inertial_coords); + } else { + cartoon_partial_derivatives_apply<2, ResultTags, VariableTags, Dim, + DerivativeFrame>( + du, u, mesh, inverse_jacobian_3d, inertial_coords); + } + } else { + ERROR("Bases do not match valid Cartoon pattern."); + } +} + template void partial_derivatives( @@ -253,8 +678,8 @@ partial_derivatives( Variables, DerivativeFrame>> partial_derivatives_of_u(u.number_of_grid_points()); - partial_derivatives(make_not_null(&partial_derivatives_of_u), - u, mesh, inverse_jacobian); + partial_derivatives(make_not_null(&partial_derivatives_of_u), u, mesh, + inverse_jacobian); return partial_derivatives_of_u; } diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PartialDerivatives.cpp b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PartialDerivatives.cpp index 33dcd7fa89b7..b1f6fc690484 100644 --- a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PartialDerivatives.cpp +++ b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PartialDerivatives.cpp @@ -21,6 +21,8 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Index.hpp" #include "DataStructures/IndexIterator.hpp" +#include "DataStructures/Tags/TempTensor.hpp" +#include "DataStructures/Tensor/Metafunctions.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" #include "DataStructures/VariablesTag.hpp" @@ -648,6 +650,360 @@ void test_partial_derivatives_spherical_shell() { } } } + +template +DataVector cartoon_func(const tnsr::I& coords) { + if constexpr (Spherical) { + // a radial function, f(r) = f(x) because the computational domain is the x + // axis + return 0.01 * pow(get<0>(coords), 4) + 0.3 * pow(get<0>(coords), 3) - + 0.1 * pow(get<0>(coords), 2) - 2.0 * get<0>(coords) - 1.5; + } else { + // an axially symmetric function about the y axis, + // f(\sqrt{x^2 + z^2}, y) = f(x, y) because the computational domain is the + // x-y plane + return square(get<1>(coords)) + square(get<0>(coords)) * get<1>(coords); + } +} + +template +DataVector cartoon_dfunc(size_t deriv_index, + const tnsr::I& coords) { + if constexpr (Spherical) { + if (deriv_index == 0) { + return 0.04 * pow(get<0>(coords), 3) + 0.9 * pow(get<0>(coords), 2) - + 0.2 * get<0>(coords) - 2.0; + } else { + return {get<0>(coords).size(), 0.0}; + } + } else { + if (deriv_index == 0) { + return 2.0 * get<0>(coords) * get<1>(coords); + } else if (deriv_index == 1) { + return 2.0 * get<1>(coords) + square(get<0>(coords)); + } else { + return {get<0>(coords).size(), 0.0}; + } + } +} + +template +void test_cartoon(const double x_start) { + Mesh<3> mesh; + tnsr::I coords; + InverseJacobian + inv_jacobian; + + using Identity1D = domain::CoordinateMaps::Identity<1>; + const Identity1D identity_cartoon_map; + + if constexpr (Spherical) { + // spherical symmetry + const size_t num_grid_pts = 8; + const double x_end = 4.0; + + mesh = Mesh<3>{{{num_grid_pts, 1, 1}}, + {{Spectral::Basis::Legendre, Spectral::Basis::Cartoon, + Spectral::Basis::Cartoon}}, + {{Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::SphericalSymmetry, + Spectral::Quadrature::SphericalSymmetry}}}; + + const Affine affine_x_map(-1.0, 1.0, x_start, x_end); + + using Cartoon_map_combination = + domain::CoordinateMaps::ProductOf3Maps; + const domain::CoordinateMap + map{{affine_x_map, identity_cartoon_map, identity_cartoon_map}}; + inv_jacobian = map.inv_jacobian(logical_coordinates(mesh)); + coords = map(logical_coordinates(mesh)); + } else { + // axial symmetry + const size_t num_x_grid_pts = 6; + const double x_end = 3.25; + const size_t num_y_grid_pts = 7; + const double y_start = -2.5; + const double y_end = 4.0; + + mesh = Mesh<3>{{{num_x_grid_pts, num_y_grid_pts, 1}}, + {{Spectral::Basis::Legendre, Spectral::Basis::Legendre, + Spectral::Basis::Cartoon}}, + {{Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::GaussLobatto, + Spectral::Quadrature::AxialSymmetry}}}; + + const Affine affine_x_map(-1.0, 1.0, x_start, x_end); + const Affine affine_y_map(-1.0, 1.0, y_start, y_end); + + using Cartoon_map_combination = + domain::CoordinateMaps::ProductOf3Maps; + const domain::CoordinateMap + map{{affine_x_map, affine_y_map, identity_cartoon_map}}; + inv_jacobian = map.inv_jacobian(logical_coordinates(mesh)); + coords = map(logical_coordinates(mesh)); + } + + using Tempabc = + ::Tags::TempTensor<0, tnsr::abc>; + + using VarTags = + tmpl::list<::Tags::TempScalar<0>, ::Tags::Tempa<0, 3, Frame::Inertial>, + ::Tags::TempA<0, 3, Frame::Inertial>, + ::Tags::TempAb<0, 3, Frame::Inertial>, + ::Tags::Tempab<0, 3, Frame::Inertial>, Tempabc>; + + Variables vars{mesh.number_of_grid_points()}; + + using d_VarTags = + tmpl::transform, Frame::Grid>>; + + Variables expected_d_vars{mesh.number_of_grid_points()}; + + // Here we create "prefactors", which serve to fill out our tensors by being + // multiplied by our cartoon_func(), which themselves make our tensors have + // some nontrivial spatial derivative + // The point of these prefactors is to ensure the tensors respect the + // symmetry of the spacetime: it is not sufficient that each component + // follows the symmetry, rather the entire tensor must satisfy + // \mathcal{L}_\xi (tensor) = 0. Each rank has it's own form of prefactor + Variables prefactor_vars{mesh.number_of_grid_points()}; + + using TempiA = + ::Tags::TempTensor<0, tnsr::iA>; + using TempiAb = + ::Tags::TempTensor<0, tnsr::iAb>; + using Tempiab = + ::Tags::TempTensor<0, tnsr::iab>; + using Tempiabc = + ::Tags::TempTensor<0, TensorMetafunctions::prepend_spatial_index< + tnsr::abc, 3, + UpLo::Lo, Frame::Inertial>>; + + // prepending a spatial index to VarTags (= PrefactorVarTags if it existed) + using d_PrefactorVarTags = tmpl::list<::Tags::Tempi<0, 3, Frame::Inertial>, + ::Tags::Tempia<0, 3, Frame::Inertial>, + TempiA, TempiAb, Tempiab, Tempiabc>; + + Variables d_prefactor_vars{mesh.number_of_grid_points()}; + + const size_t dv_size = get<0>(coords).size(); + + // in the following, the time component of the indices is arbirarily taken + // to be 1 + auto& scalar = get<::Tags::TempScalar<0>>(prefactor_vars); + auto& d_scalar = get<::Tags::Tempi<0, 3>>(d_prefactor_vars); + // scalar, doing nothing + scalar.get() = DataVector(dv_size, 1.0); + // \partial_i of 1 is zero + for (size_t i = 0; i < index_dim<0>(d_scalar); ++i) { + d_scalar.get(i) = DataVector(dv_size, 0.0); + } + + auto& one_form = get<::Tags::Tempa<0, 3>>(prefactor_vars); + auto& d_one_form = get<::Tags::Tempia<0, 3>>(d_prefactor_vars); + auto& vector = get<::Tags::TempA<0, 3>>(prefactor_vars); + auto& d_vector = get(d_prefactor_vars); + if constexpr (Spherical) { + // spherical case, vector/one form, using x^a + get<0>(vector) = get<0>(one_form) = DataVector(dv_size, 1.0); + for (size_t i = 1; i < index_dim<0>(vector); ++i) { + vector.get(i) = one_form.get(i) = coords.get(i - 1); + } + // partial_i of x_a is \delta_ia + for (size_t i = 0; i < index_dim<0>(d_vector); ++i) { + for (size_t a = 0; a < index_dim<1>(d_vector); ++a) { + if ((i + 1) == a and a != 0) { + d_vector.get(i, a) = d_one_form.get(i, a) = DataVector(dv_size, 1.0); + } else { + d_vector.get(i, a) = d_one_form.get(i, a) = DataVector(dv_size, 0.0); + } + } + } + } else { + // axial case, vector/one form, using x^a = (0, -z, 0, x) (pure rotation) + get<0>(vector) = get<0>(one_form) = DataVector(dv_size, 0.0); + get<1>(vector) = get<1>(one_form) = -1.0 * coords.get(2); + get<2>(vector) = get<2>(one_form) = DataVector(dv_size, 0.0); + get<3>(vector) = get<3>(one_form) = get<0>(coords); + + // \partial_i of x_a is (0, \delta_i2, 0, \delta_i0) + for (size_t i = 0; i < index_dim<0>(d_vector); ++i) { + for (size_t a = 0; a < index_dim<1>(d_vector); ++a) { + if (i == 2 and a == 1) { + d_vector.get(i, a) = d_one_form.get(i, a) = DataVector(dv_size, -1.0); + } else if (i == 0 and a == 3) { + d_vector.get(i, a) = d_one_form.get(i, a) = DataVector(dv_size, 1.0); + } else { + d_vector.get(i, a) = d_one_form.get(i, a) = DataVector(dv_size, 0.0); + } + } + } + } + + auto& mixed = get<::Tags::TempAb<0, 3>>(prefactor_vars); + auto& d_mixed = get(d_prefactor_vars); + auto& rank2 = get<::Tags::Tempab<0, 3>>(prefactor_vars); + auto& d_rank2 = get(d_prefactor_vars); + // filling space with (essenially) projector to tangent space of sphere + // P_ab = \delta_ab + x_a x_b + // can have arbitrary function in from of each term, not doing here + // (the real projector is \delta_ij - x_i x_j / r^2) + for (size_t a = 0; a < index_dim<0>(rank2); ++a) { + for (size_t b = 0; b < index_dim<1>(rank2); ++b) { + if (a == 0 and b == 0) { + mixed.get(a, b) = rank2.get(a, b) = DataVector(dv_size, 1.0); + } else if (a == 0) { + mixed.get(a, b) = rank2.get(a, b) = coords.get(b - 1); + } else if (b == 0) { + mixed.get(a, b) = rank2.get(a, b) = coords.get(a - 1); + } else { + if (a == b) { + mixed.get(a, b) = rank2.get(a, b) = + DataVector(dv_size, 1.0) + coords.get(a - 1) * coords.get(b - 1); + } else { + mixed.get(a, b) = rank2.get(a, b) = + coords.get(a - 1) * coords.get(b - 1); + } + } + } + } + // \partial_i of (\delta_ab + x_a x_b) is (x_a \delta_ib + x_b \delta_ia) + for (size_t i = 0; i < index_dim<0>(d_rank2); ++i) { + for (size_t a = 0; a < index_dim<1>(d_rank2); ++a) { + for (size_t b = 0; b < index_dim<2>(d_rank2); ++b) { + d_mixed.get(i, a, b) = d_rank2.get(i, a, b) = DataVector(dv_size, 0.0); + if ((i + 1) == b) { + if (a == 0) { + d_mixed.get(i, a, b) += 1.0; + d_rank2.get(i, a, b) += 1.0; + } else { + d_mixed.get(i, a, b) += coords.get(a - 1); + d_rank2.get(i, a, b) += coords.get(a - 1); + } + } + if ((i + 1) == a) { + if (b == 0) { + d_mixed.get(i, a, b) += 1.0; + d_rank2.get(i, a, b) += 1.0; + } else { + d_mixed.get(i, a, b) += coords.get(b - 1); + d_rank2.get(i, a, b) += coords.get(b - 1); + } + } + } + } + } + + auto& rank3 = get(prefactor_vars); + auto& d_rank3 = get(d_prefactor_vars); + // Rank 3 is important because of \Phi_iab, the spatial derivative of the + // metric, which though not a tensor mathematically still must obey + // \mathcal{L}_\xi \Phi_{iab} = 0 + // Using the form + // x_a x_b x_c + \delta_ab x_c \delta_ac x_b \delta_bc x_a + for (size_t a = 0; a < index_dim<0>(rank3); ++a) { + for (size_t b = 0; b < index_dim<1>(rank3); ++b) { + for (size_t c = 0; c < index_dim<2>(rank3); ++c) { + rank3.get(a, b, c) = + (a == 0 ? DataVector(dv_size, 1.0) : coords.get(a - 1)) * + (b == 0 ? DataVector(dv_size, 1.0) : coords.get(b - 1)) * + (c == 0 ? DataVector(dv_size, 1.0) : coords.get(c - 1)); + if (a == b) { + rank3.get(a, b, c) += + (c == 0 ? DataVector(dv_size, 1.0) : coords.get(c - 1)); + } + if (a == c) { + rank3.get(a, b, c) += + (b == 0 ? DataVector(dv_size, 1.0) : coords.get(b - 1)); + } + if (b == c) { + rank3.get(a, b, c) += + (a == 0 ? DataVector(dv_size, 1.0) : coords.get(a - 1)); + } + } + } + } + // \partial_i of (x_a x_b x_c + \delta_ab x_c \delta_ac x_b \delta_bc x_a) is + // \delta_ia x_b x_c + \delta_ib x_a x_c + \delta_ic x_a x_b + + // \delta_ab \delta_ic + \delta_ac \delta_ib + \delta_bc \delta_ia + for (size_t i = 0; i < index_dim<0>(d_rank3); ++i) { + for (size_t a = 0; a < index_dim<1>(d_rank3); ++a) { + for (size_t b = 0; b < index_dim<2>(d_rank3); ++b) { + for (size_t c = 0; c < index_dim<3>(d_rank3); ++c) { + d_rank3.get(i, a, b, c) = DataVector(dv_size, 0.0); + if ((i + 1) == a) { + d_rank3.get(i, a, b, c) += + (b == 0 ? DataVector(dv_size, 1.0) : coords.get(b - 1)) * + (c == 0 ? DataVector(dv_size, 1.0) : coords.get(c - 1)); + } + if ((i + 1) == b) { + d_rank3.get(i, a, b, c) += + (a == 0 ? DataVector(dv_size, 1.0) : coords.get(a - 1)) * + (c == 0 ? DataVector(dv_size, 1.0) : coords.get(c - 1)); + } + if ((i + 1) == c) { + d_rank3.get(i, a, b, c) += + (a == 0 ? DataVector(dv_size, 1.0) : coords.get(a - 1)) * + (b == 0 ? DataVector(dv_size, 1.0) : coords.get(b - 1)); + } + if (a == b and (i + 1) == c) { + d_rank3.get(i, a, b, c) += 1.0; + } + if (a == c and (i + 1) == b) { + d_rank3.get(i, a, b, c) += 1.0; + } + if (b == c and (i + 1) == a) { + d_rank3.get(i, a, b, c) += 1.0; + } + } + } + } + } + + tmpl::for_each([&vars, &prefactor_vars, &expected_d_vars, + &d_prefactor_vars, &coords]( + tmpl::type_ /*meta*/) { + auto& tensor = get(vars); + const auto& prefactor_tensor = get(prefactor_vars); + using deriv_tag = tmpl::apply< + tmpl::bind<::Tags::deriv, tmpl::_1, tmpl::size_t<3>, Frame::Grid>, + tensor_tag>; + auto& d_tensor = get(expected_d_vars); + const auto& d_prefactor_tensor = + get>>( + d_prefactor_vars); + + for (size_t storage_index = 0; storage_index < tensor.size(); + ++storage_index) { + tensor[storage_index] = + prefactor_tensor[storage_index] * cartoon_func(coords); + + const auto input_index = tensor.get_tensor_index(storage_index); + for (size_t deriv_index = 0; deriv_index < 3; ++deriv_index) { + const auto output_index = prepend(input_index, size_t{deriv_index}); + + d_tensor.get(output_index) = + prefactor_tensor.get(input_index) * + cartoon_dfunc(deriv_index, coords) + + cartoon_func(coords) * + d_prefactor_tensor.get(output_index); + } + } + }); + + const auto to_inertial = + domain::CoordinateMap>{}; + Variables d_vars{mesh.number_of_grid_points()}; + cartoon_partial_derivatives(make_not_null(&d_vars), vars, mesh, inv_jacobian, + to_inertial(coords)); + // `d_vars` holds our partial derivatives we want to test against the truth + const Approx local_approx = Approx::custom().epsilon(1e-10).scale(1.0); + CHECK_VARIABLES_CUSTOM_APPROX(d_vars, expected_d_vars, local_approx); +} } // namespace // [[Timeout, 60]] @@ -772,6 +1128,13 @@ SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.PartialDerivs", test_partial_derivatives_3d, one_var>(mesh_3d); + // testing with L'Hopital + test_cartoon(0.0); + test_cartoon(0.0); + // testing without L'Hopital + test_cartoon(1.0); + test_cartoon(1.0); + TestHelpers::db::test_prefix_tag< Tags::deriv, tmpl::size_t<3>, Frame::Grid>>( "deriv(Var1)");