From c3caa14650778fd8c94f31dde04ed7bcb120d1e6 Mon Sep 17 00:00:00 2001 From: Abdelhadi KARA Date: Tue, 3 Dec 2024 13:17:23 +0100 Subject: [PATCH] Split mappings so they represent one equation Closes #409 See merge request gysela-developpers/gyselalibxx!779 -------------------------------------------- Co-authored-by: Emily Bourne --- .../geometryRTheta/diocotron/diocotron.cpp | 32 ++- .../vortex_merger/vortex_merger.cpp | 31 ++- .../advection/advection_domain.hpp | 48 ++-- .../advection/spline_foot_finder.hpp | 2 +- .../connectivity/onion_patch_locator.hpp | 171 +++++++------ .../connectivity/utils_patch_locators.hpp | 13 +- .../spline/multipatch_spline_evaluator_2d.hpp | 2 + .../advection_2d_rp/advection_all_tests.cpp | 63 +++-- .../advection_selected_test.cpp | 61 +++-- .../advection_simulation_utils.hpp | 129 ++++++---- .../advection_2d_rp/test_cases.hpp | 16 +- .../advection_field_gtest.cpp | 30 +-- .../test_cases_adv_field.hpp | 5 +- .../polar_poisson/polarpoissonfemsolver.cpp | 4 +- .../polar_poisson/test_cases.cpp | 24 +- .../polar_poisson/test_cases.hpp | 1 + .../quadrature/tests_L1_and_L2_norms.cpp | 2 +- .../patch_locator_onion_shape_2patches.cpp | 94 ++++--- .../spline/multipatch_spline_evaluator.cpp | 35 ++- .../spline/null_extrapolation_rules.cpp | 1 + .../spline/spline_testing_tools.hpp | 10 +- .../sll/mapping/cartesian_to_circular.hpp | 233 ++++++++++++++++++ .../sll/mapping/cartesian_to_czarny.hpp | 153 ++++++++++++ .../sll/mapping/circular_to_cartesian.hpp | 31 +-- .../sll/mapping/czarny_to_cartesian.hpp | 43 +--- .../sll/mapping/discrete_to_cartesian.hpp | 2 +- .../sll/mapping/inverse_jacobian_matrix.hpp | 2 +- .../tests/mapping_execution_space_access.cpp | 88 ++++--- vendor/sll/tests/mapping_jacobian.cpp | 12 +- .../tests/mapping_jacobian_matrix_coef.cpp | 6 +- vendor/sll/tests/metric_tensor.cpp | 8 +- vendor/sll/tests/polar_bsplines.cpp | 2 +- vendor/sll/tests/polar_splines.cpp | 4 +- vendor/sll/tests/pseudo_cartesian_matrix.cpp | 4 +- vendor/sll/tests/refined_discrete_mapping.cpp | 4 +- 35 files changed, 952 insertions(+), 414 deletions(-) create mode 100644 vendor/sll/include/sll/mapping/cartesian_to_circular.hpp create mode 100644 vendor/sll/include/sll/mapping/cartesian_to_czarny.hpp diff --git a/simulations/geometryRTheta/diocotron/diocotron.cpp b/simulations/geometryRTheta/diocotron/diocotron.cpp index 3034ff7b2..1bfabb6dd 100644 --- a/simulations/geometryRTheta/diocotron/diocotron.cpp +++ b/simulations/geometryRTheta/diocotron/diocotron.cpp @@ -6,8 +6,8 @@ #include +#include #include -#include #include #include @@ -51,7 +51,8 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver< SplineRThetaEvaluatorNullBound>; using DiscreteMappingBuilder = DiscreteToCartesianBuilder; -using Mapping = CircularToCartesian; +using LogicalToPhysicalMapping = CircularToCartesian; +using PhysicalToLogicalMapping = CartesianToCircular; namespace fs = std::filesystem; @@ -111,10 +112,11 @@ int main(int argc, char** argv) ddc::PeriodicExtrapolationRule(), ddc::PeriodicExtrapolationRule()); - const Mapping mapping; + const LogicalToPhysicalMapping to_physical_mapping; + const PhysicalToLogicalMapping to_logical_mapping; DiscreteMappingBuilder const discrete_mapping_builder( Kokkos::DefaultHostExecutionSpace(), - mapping, + to_physical_mapping, builder, spline_evaluator_extrapol); DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); @@ -161,12 +163,16 @@ int main(int argc, char** argv) PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator); - AdvectionPhysicalDomain advection_domain(mapping); + AdvectionPhysicalDomain advection_domain(to_physical_mapping, to_logical_mapping); - SplineFootFinder - find_feet(time_stepper, advection_domain, mapping, builder, spline_evaluator_extrapol); + SplineFootFinder find_feet( + time_stepper, + advection_domain, + to_physical_mapping, + builder, + spline_evaluator_extrapol); - BslAdvectionRTheta advection_operator(interpolator, find_feet, mapping); + BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping); @@ -195,7 +201,7 @@ int main(int argc, char** argv) // --- Predictor corrector operator --------------------------------------------------------------- #if defined(PREDCORR) BslPredCorrRTheta predcorr_operator( - mapping, + to_physical_mapping, advection_operator, builder, spline_evaluator, @@ -203,7 +209,7 @@ int main(int argc, char** argv) #elif defined(EXPLICIT_PREDCORR) BslExplicitPredCorrRTheta predcorr_operator( advection_domain, - mapping, + to_physical_mapping, advection_operator, mesh_rp, builder, @@ -213,7 +219,7 @@ int main(int argc, char** argv) #elif defined(IMPLICIT_PREDCORR) BslImplicitPredCorrRTheta predcorr_operator( advection_domain, - mapping, + to_physical_mapping, advection_operator, mesh_rp, builder, @@ -267,10 +273,10 @@ int main(int argc, char** argv) host_t> coords_y(mesh_rp); host_t jacobian(mesh_rp); ddc::for_each(mesh_rp, [&](IdxRTheta const irp) { - CoordXY coords_xy = mapping(ddc::coordinate(irp)); + CoordXY coords_xy = to_physical_mapping(ddc::coordinate(irp)); coords_x(irp) = ddc::select(coords_xy); coords_y(irp) = ddc::select(coords_xy); - jacobian(irp) = mapping.jacobian(ddc::coordinate(irp)); + jacobian(irp) = to_physical_mapping.jacobian(ddc::coordinate(irp)); }); diff --git a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp index 5006a3eed..19e32b0f7 100644 --- a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp +++ b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp @@ -51,7 +51,8 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver< SplineRThetaEvaluatorNullBound>; using DiscreteMappingBuilder = DiscreteToCartesianBuilder; -using CircularMapping = CircularToCartesian; +using LogicalToPhysicalMapping = CircularToCartesian; +using PhysicalToLogicalMapping = CartesianToCircular; } // end namespace @@ -109,10 +110,11 @@ int main(int argc, char** argv) ddc::PeriodicExtrapolationRule(), ddc::PeriodicExtrapolationRule()); - const CircularMapping mapping; + const LogicalToPhysicalMapping to_physical_mapping; + const PhysicalToLogicalMapping to_logical_mapping; DiscreteMappingBuilder const discrete_mapping_builder( Kokkos::DefaultHostExecutionSpace(), - mapping, + to_physical_mapping, builder, spline_evaluator_extrapol); DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); @@ -139,12 +141,16 @@ int main(int argc, char** argv) PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator); - AdvectionPhysicalDomain advection_domain(mapping); + AdvectionPhysicalDomain advection_domain(to_physical_mapping, to_logical_mapping); - SplineFootFinder - find_feet(time_stepper, advection_domain, mapping, builder, spline_evaluator_extrapol); + SplineFootFinder find_feet( + time_stepper, + advection_domain, + to_physical_mapping, + builder, + spline_evaluator_extrapol); - BslAdvectionRTheta advection_operator(interpolator, find_feet, mapping); + BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping); @@ -173,7 +179,7 @@ int main(int argc, char** argv) // --- Predictor corrector operator --------------------------------------------------------------- BslImplicitPredCorrRTheta predcorr_operator( advection_domain, - mapping, + to_physical_mapping, advection_operator, grid, builder, @@ -225,23 +231,24 @@ int main(int argc, char** argv) host_t> coords_y(grid); host_t jacobian(grid); ddc::for_each(grid, [&](IdxRTheta const irp) { - CoordXY coords_xy = mapping(ddc::coordinate(irp)); + CoordXY coords_xy = to_physical_mapping(ddc::coordinate(irp)); coords_x(irp) = ddc::select(coords_xy); coords_y(irp) = ddc::select(coords_xy); - jacobian(irp) = mapping.jacobian(ddc::coordinate(irp)); + jacobian(irp) = to_physical_mapping.jacobian(ddc::coordinate(irp)); }); double const tau(1e-10); double const phi_max(1.); - VortexMergerEquilibria equilibrium(mapping, grid, builder, spline_evaluator, poisson_solver); + VortexMergerEquilibria + equilibrium(to_physical_mapping, grid, builder, spline_evaluator, poisson_solver); std::function const function = [&](double const x) { return x * x; }; host_t rho_eq(grid); equilibrium.set_equilibrium(rho_eq, function, phi_max, tau); - VortexMergerDensitySolution solution(mapping); + VortexMergerDensitySolution solution(to_physical_mapping); host_t rho(grid); solution.set_initialisation(rho, rho_eq, eps, sigma, x_star_1, y_star_1, x_star_2, y_star_2); diff --git a/src/geometryRTheta/advection/advection_domain.hpp b/src/geometryRTheta/advection/advection_domain.hpp index 8fcc840f1..7a45234f1 100644 --- a/src/geometryRTheta/advection/advection_domain.hpp +++ b/src/geometryRTheta/advection/advection_domain.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -35,7 +36,7 @@ * @see BslAdvectionRTheta * @see IFootFinder */ -template +template class AdvectionDomain { public: @@ -62,8 +63,8 @@ class AdvectionDomain * * */ -template -class AdvectionPhysicalDomain : public AdvectionDomain +template +class AdvectionPhysicalDomain : public AdvectionDomain { public: /** @@ -80,16 +81,23 @@ class AdvectionPhysicalDomain : public AdvectionDomain using CoordXY_adv = Coord; private: - Mapping const& m_mapping; + LogicalToPhysicalMapping const& m_to_cartesian_mapping; + PhysicalToLogicalMapping const& m_to_curvilinear_mapping; public: /** * @brief Instantiate a AdvectionPhysicalDomain advection domain. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping from the logical domain to the physical domain. + * @param[in] to_logical_mapping + * The mapping from the physical domain to the logical domain. */ - AdvectionPhysicalDomain(Mapping const& mapping) : m_mapping(mapping) {}; + AdvectionPhysicalDomain( + LogicalToPhysicalMapping const& to_physical_mapping, + PhysicalToLogicalMapping const& to_logical_mapping) + : m_to_cartesian_mapping(to_physical_mapping) + , m_to_curvilinear_mapping(to_logical_mapping) {}; ~AdvectionPhysicalDomain() {}; /** @@ -136,18 +144,18 @@ class AdvectionPhysicalDomain : public AdvectionDomain using namespace ddc; IdxRangeRTheta const idx_range_rp = get_idx_range(feet_coords_rp); - CoordXY coord_center(m_mapping(CoordRTheta(0, 0))); + CoordXY coord_center(m_to_cartesian_mapping(CoordRTheta(0, 0))); ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(feet_coords_rp(irp)); - CoordXY const coord_xy = m_mapping(coord_rp); + CoordXY const coord_xy = m_to_cartesian_mapping(coord_rp); CoordXY const feet_xy = coord_xy - dt * advection_field(irp); if (norm_inf(feet_xy - coord_center) < 1e-15) { feet_coords_rp(irp) = CoordRTheta(0, 0); } else { - feet_coords_rp(irp) = m_mapping(feet_xy); + feet_coords_rp(irp) = m_to_curvilinear_mapping(feet_xy); ddc::select(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range( ddc::select(feet_coords_rp(irp)), IdxRangeTheta(idx_range_rp)); @@ -203,8 +211,8 @@ class AdvectionPhysicalDomain : public AdvectionDomain * @see Curvilinear2DToCartesian * @see DiscreteToCartesian */ -template -class AdvectionPseudoCartesianDomain : public AdvectionDomain +template +class AdvectionPseudoCartesianDomain : public AdvectionDomain { public: /** @@ -227,21 +235,23 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain private: - Mapping const& m_mapping; + LogicalToPhysicalMapping const& m_to_cartesian_mapping; double m_epsilon; public: /** * @brief Instantiate an AdvectionPseudoCartesianDomain advection domain. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping from the logical domain to the physical domain. * @param[in] epsilon * @f$ \varepsilon @f$ parameter used for the linearization of the * advection field around the central point. */ - AdvectionPseudoCartesianDomain(Mapping const& mapping, double epsilon = 1e-12) - : m_mapping(mapping) + AdvectionPseudoCartesianDomain( + LogicalToPhysicalMapping const& to_physical_mapping, + double epsilon = 1e-12) + : m_to_cartesian_mapping(to_physical_mapping) , m_epsilon(epsilon) {}; ~AdvectionPseudoCartesianDomain() {}; @@ -286,10 +296,11 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain host_t> const& advection_field, double const dt) const { - static_assert(!std::is_same_v>); + static_assert( + !std::is_same_v>); IdxRangeRTheta const idx_range_rp = get_idx_range(advection_field); - CircularToCartesian const pseudo_Cartesian_mapping; + CircularToCartesian const pseudo_Cartesian_mapping; CoordXY_adv const center_xy_pseudo_cart = CoordXY_adv(pseudo_Cartesian_mapping(CoordRTheta(0., 0.))); @@ -302,7 +313,8 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain if (norm_inf(feet_xy_pseudo_cart - center_xy_pseudo_cart) < 1e-15) { feet_coords_rp(irp) = CoordRTheta(0, 0); } else { - feet_coords_rp(irp) = pseudo_Cartesian_mapping(feet_xy_pseudo_cart); + CartesianToCircular const inv_pseudo_Cartesian_mapping; + feet_coords_rp(irp) = inv_pseudo_Cartesian_mapping(feet_xy_pseudo_cart); ddc::select(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range( ddc::select(feet_coords_rp(irp)), IdxRangeTheta(idx_range_rp)); diff --git a/src/geometryRTheta/advection/spline_foot_finder.hpp b/src/geometryRTheta/advection/spline_foot_finder.hpp index 19817e21a..c6fbeea2f 100644 --- a/src/geometryRTheta/advection/spline_foot_finder.hpp +++ b/src/geometryRTheta/advection/spline_foot_finder.hpp @@ -34,7 +34,7 @@ class SplineFootFinder : public IFootFinder */ using Y_adv = typename AdvectionDomain::Y_adv; - using CircularToPseudoCartesian = CircularToCartesian; + using CircularToPseudoCartesian = CircularToCartesian; using AdvectionMapping = CartesianToPseudoCartesian; TimeStepper const& m_time_stepper; diff --git a/src/multipatch/connectivity/onion_patch_locator.hpp b/src/multipatch/connectivity/onion_patch_locator.hpp index 45bb61d97..28df61e3e 100644 --- a/src/multipatch/connectivity/onion_patch_locator.hpp +++ b/src/multipatch/connectivity/onion_patch_locator.hpp @@ -12,53 +12,65 @@ /** - * @brief Patch locator specialised for "onion" geometry. - * - * We define an "onion" geometry a set of patches mapping to the + * @brief Patch locator specialised for "onion" geometry. + * + * We define an "onion" geometry a set of patches mapping to the * physical domain in a shape of concentrical rings. The first patch - * is supposed to be a disk containing the O-point. The other patches - * are ordered as concentrical rings drawing away from the O-point. - * The order of patches is made by MultipatchIdxRanges and it is important + * is supposed to be a disk containing the O-point. The other patches + * are ordered as concentrical rings drawing away from the O-point. + * The order of patches is made by MultipatchIdxRanges and it is important * for the dichotomy method. - * - * We also suppose that we can define a global logical grid that we can split - * into the different logical grids of the patches. - * - * This operator locates on which patch a given physical coordinate is. - * - * @warning The operator can works on GPU or CPU according to the given - * execution space ExecSpace. The ExecSpace by default is device. - * The constructor will still need to be called from CPU, but the operator() - * needs to be called from the given ExecSpace. - * + * + * We also suppose that we can define a global logical grid that we can split + * into the different logical grids of the patches. + * + * This operator locates on which patch a given physical coordinate is. + * + * @warning The operator can works on GPU or CPU according to the given + * execution space ExecSpace. The ExecSpace by default is device. + * The constructor will still need to be called from CPU, but the operator() + * needs to be called from the given ExecSpace. + * * @anchor OnionPatchLocatorImplementation - * - * @tparam MultipatchIdxRanges A MultipatchType type containing the 2D index ranges - * on each patch. - * @tparam Mapping A mapping type for all the patches. + * + * @tparam MultipatchIdxRanges A MultipatchType type containing the 2D index ranges + * on each patch. + * @tparam LogicalToPhysicalMapping A mapping type for all the patches. * @tparam ExecSpace The space (CPU/GPU) where the calculations are carried out. - * By default it is on device. + * By default it is on device. */ -template +template < + class MultipatchIdxRanges, + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class ExecSpace = Kokkos::DefaultExecutionSpace> class OnionPatchLocator; /** - * @brief Patch locator specialised for "onion" geometry. - * + * @brief Patch locator specialised for "onion" geometry. + * * See @ref OnionPatchLocatorImplementation - * + * * @tparam ExecSpace The space (CPU/GPU) where the calculations are carried out. - * @tparam Patches Patch types. Their order is important. - * @tparam Mapping A mapping type for all the patches. + * @tparam Patches Patch types. Their order is important. + * @tparam LogicalToPhysicalMapping A mapping type for all the patches. */ -template -class OnionPatchLocator, Mapping, ExecSpace> +template < + class... Patches, + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class ExecSpace> +class OnionPatchLocator< + MultipatchType, + LogicalToPhysicalMapping, + PhysicalToLogicalMapping, + ExecSpace> { - using X = typename Mapping::cartesian_tag_x; - using Y = typename Mapping::cartesian_tag_y; - using R = typename Mapping::curvilinear_tag_r; - using Theta = typename Mapping::curvilinear_tag_theta; + using X = typename LogicalToPhysicalMapping::cartesian_tag_x; + using Y = typename LogicalToPhysicalMapping::cartesian_tag_y; + using R = typename LogicalToPhysicalMapping::curvilinear_tag_r; + using Theta = typename LogicalToPhysicalMapping::curvilinear_tag_theta; static_assert(Theta::PERIODIC, "Theta dimension must be periodic."); @@ -84,7 +96,7 @@ class OnionPatchLocator, Mapping, Ex static constexpr std::size_t n_patches = ddc::type_seq_size_v; static_assert( - std::is_invocable_r_v, Mapping, Coord>, + std::is_invocable_r_v, PhysicalToLogicalMapping, Coord>, "The mapping has to contain an operator from the physical domain to the logical " "domain."); static_assert( @@ -94,23 +106,29 @@ class OnionPatchLocator, Mapping, Ex "The mappings and the patches have to be defined on the same dimensions."); - Mapping const m_mapping; + LogicalToPhysicalMapping const m_to_cartesian_mapping; + PhysicalToLogicalMapping const m_to_curvilinear_mapping; MultipatchIdxRanges const m_all_idx_ranges; Kokkos::View*, typename ExecSpace::memory_space> m_radii; public: - /** - * @brief Instantiante the operator with MultipatchType of index ranges and - * a mapping on all the patches. - * + /** + * @brief Instantiante the operator with MultipatchType of index ranges and + * a mapping on all the patches. + * * The order of the elements in the tuple or the MultipatchType doesn't matter. - * + * * @param all_idx_ranges A MultipatchType of index ranges defined on the logical domain of each patch. - * @param mapping Mapping from the logical domains of every patch to the global physical domain. + * @param to_physical_mapping Mapping from the logical domains of every patch to the global physical domain. + * @param to_logical_mapping Mapping from the global physical domain to the logical domain of every patch. */ - OnionPatchLocator(MultipatchIdxRanges const& all_idx_ranges, Mapping const& mapping) - : m_mapping(mapping) + OnionPatchLocator( + MultipatchIdxRanges const& all_idx_ranges, + LogicalToPhysicalMapping const& to_physical_mapping, + PhysicalToLogicalMapping const& to_logical_mapping) + : m_to_cartesian_mapping(to_physical_mapping) + , m_to_curvilinear_mapping(to_logical_mapping) , m_all_idx_ranges(all_idx_ranges) , m_radii("m_radii", n_patches + 1) { @@ -121,15 +139,15 @@ class OnionPatchLocator, Mapping, Ex ~OnionPatchLocator() = default; /** - * @brief Get the patch where the given physical coordinate is. - * - * We use a dichotomy method to find the patch the physical coordinate is on. - * - * Each logical grid of the patches are defined on the same dimensions. + * @brief Get the patch where the given physical coordinate is. + * + * We use a dichotomy method to find the patch the physical coordinate is on. + * + * Each logical grid of the patches are defined on the same dimensions. * Knowing that, we can compare the logical coordinates between the patches - * in the dichotomy. - * - * @param coord [in] The given physical coordinate. + * in the dichotomy. + * + * @param coord [in] The given physical coordinate. * @return [int] The patch index where the physical coordinate. If the coordinate * is outside of the domain, it returns a negative value. */ @@ -140,7 +158,7 @@ class OnionPatchLocator, Mapping, Ex Coord r_min = m_radii(patch_index_min); Coord r_max = m_radii(patch_index_max + 1); - Coord r(m_mapping(coord)); + Coord r(m_to_curvilinear_mapping(coord)); KOKKOS_ASSERT(!Kokkos::isnan(double(r))); if (r > r_max) { @@ -172,43 +190,43 @@ class OnionPatchLocator, Mapping, Ex /** - * @brief Get the mapping on the given Patch. - * The function can run on device and host. - * @tparam Patch Patch type. - * @return The mapping on the given Patch. + * @brief Get the mapping on the given Patch. + * The function can run on device and host. + * @tparam Patch Patch type. + * @return The mapping on the given Patch. */ template - KOKKOS_FUNCTION Mapping get_mapping_on_patch() const + KOKKOS_FUNCTION LogicalToPhysicalMapping get_mapping_on_patch() const { - return m_mapping; + return m_to_cartesian_mapping; } /** - * @brief Get the mapping from given logical continuous dimensions. - * The function can run on device and host. - * @tparam Patch Patch type. - * @return The mapping on the given Patch. + * @brief Get the mapping from given logical continuous dimensions. + * The function can run on device and host. + * @tparam Patch Patch type. + * @return The mapping on the given Patch. */ template - KOKKOS_FUNCTION Mapping get_mapping_on_logical_dim() const + KOKKOS_FUNCTION LogicalToPhysicalMapping get_mapping_on_logical_dim() const { static_assert( (std::is_same_v)&&(std::is_same_v), "Wrong continuous dimensions."); - return m_mapping; + return m_to_cartesian_mapping; } /// @brief Get the type of the mapping on the given Patch. template - using get_mapping_on_patch_t = Mapping; + using get_mapping_on_patch_t = LogicalToPhysicalMapping; /// @brief Get the type of the mapping from given logical continuous dimensions. template - using get_mapping_on_logical_dim_t = Mapping; + using get_mapping_on_logical_dim_t = LogicalToPhysicalMapping; private: /** @brief Set the m_radii array containing all the boundary radial coordinates. - * Check if the patches are well ordered by comparing the radii. + * Check if the patches are well ordered by comparing the radii. */ void set_and_check_radii() { @@ -234,6 +252,17 @@ class OnionPatchLocator, Mapping, Ex // To help the template deduction. -template -OnionPatchLocator(MultipatchIdxRanges const& all_idx_ranges, Mapping const& mapping) - -> OnionPatchLocator; +template < + class MultipatchIdxRanges, + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class ExecSpace = Kokkos::DefaultExecutionSpace> +OnionPatchLocator( + MultipatchIdxRanges const& all_idx_ranges, + LogicalToPhysicalMapping const& to_physical_mapping, + PhysicalToLogicalMapping const& to_logical_mapping) + -> OnionPatchLocator< + MultipatchIdxRanges, + LogicalToPhysicalMapping, + PhysicalToLogicalMapping, + ExecSpace>; diff --git a/src/multipatch/connectivity/utils_patch_locators.hpp b/src/multipatch/connectivity/utils_patch_locators.hpp index ec8508e4f..606a7f282 100644 --- a/src/multipatch/connectivity/utils_patch_locators.hpp +++ b/src/multipatch/connectivity/utils_patch_locators.hpp @@ -16,8 +16,15 @@ template inline constexpr bool is_onion_patch_locator_v = is_onion_patch_locator::value; -template -struct is_onion_patch_locator> - : std::true_type +template < + class MultipatchIdxRanges, + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class ExecSpace> +struct is_onion_patch_locator> : std::true_type { }; diff --git a/src/multipatch/spline/multipatch_spline_evaluator_2d.hpp b/src/multipatch/spline/multipatch_spline_evaluator_2d.hpp index 0c0520b6e..74761153d 100644 --- a/src/multipatch/spline/multipatch_spline_evaluator_2d.hpp +++ b/src/multipatch/spline/multipatch_spline_evaluator_2d.hpp @@ -6,6 +6,8 @@ #include +#include + #include "ddc_aliases.hpp" #include "ddc_helper.hpp" #include "multipatch_type.hpp" diff --git a/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp b/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp index 557a24078..a536adde3 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp +++ b/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp @@ -11,6 +11,8 @@ #include +#include +#include #include #include #include @@ -47,30 +49,39 @@ namespace { namespace fs = std::filesystem; using CurvilinearMapping = Curvilinear2DToCartesian; -using CircularMapping = CircularToCartesian; -using CzarnyMapping = CzarnyToCartesian; +using CircularToCartMapping = CircularToCartesian; +using CzarnyToCartMapping = CzarnyToCartesian; +using CartToCircularMapping = CartesianToCircular; +using CartToCzarnyMapping = CartesianToCzarny; using DiscreteMappingBuilder = DiscreteToCartesianBuilder; } // end namespace -template +template < + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class AnalyticalMapping, + class AdvectionDomain> struct SimulationParameters { public: - Mapping const& mapping; + LogicalToPhysicalMapping const& to_physical_mapping; + PhysicalToLogicalMapping const to_logical_mapping; AnalyticalMapping const& analytical_mapping; AdvectionDomain const& advection_domain; std::string mapping_name; std::string domain_name; using IdxRangeSimulationAdvection = AdvectionDomain; SimulationParameters( - Mapping const& map, + LogicalToPhysicalMapping const& map, + PhysicalToLogicalMapping const& rev_map, AnalyticalMapping const& a_map, AdvectionDomain const& advection_dom, std::string m_name, std::string dom_name) - : mapping(map) + : to_physical_mapping(map) + , to_logical_mapping(rev_map) , analytical_mapping(a_map) , advection_domain(advection_dom) , mapping_name(m_name) @@ -198,7 +209,8 @@ void run_simulations_with_methods( std::string output_stem = output_stream.str(); simulate_the_3_simulations( - sim.mapping, + sim.to_physical_mapping, + sim.to_logical_mapping, sim.analytical_mapping, params.grid, num.time_stepper, @@ -281,7 +293,7 @@ int main(int argc, char** argv) // DEFINITION OF OPERATORS ------------------------------------------------------------------ - // --- Builders for the test function and the mapping: + // --- Builders for the test function and the to_physical_mapping: SplineRThetaBuilder const builder(grid); // --- Evaluator for the test function: @@ -310,41 +322,50 @@ int main(int argc, char** argv) // SET THE DIFFERENT PARAMETERS OF THE TESTS ------------------------------------------------ - CircularMapping const circ_map; - CzarnyMapping const czarny_map(0.3, 1.4); + CircularToCartMapping const from_circ_map; + CartesianToCircular to_circ_map; + CzarnyToCartMapping const from_czarny_map(0.3, 1.4); + CartesianToCzarny const to_czarny_map(0.3, 1.4); DiscreteMappingBuilder const discrete_czarny_map_builder( Kokkos::DefaultHostExecutionSpace(), - czarny_map, + from_czarny_map, builder, spline_evaluator_extrapol); DiscreteToCartesian const discrete_czarny_map = discrete_czarny_map_builder(); - AdvectionPhysicalDomain const physical_circular_mapping(circ_map); - AdvectionPhysicalDomain const physical_czarny_mapping(czarny_map); - AdvectionPseudoCartesianDomain const pseudo_cartesian_czarny_mapping(czarny_map); + AdvectionPhysicalDomain const + physical_circular_mapping(from_circ_map, to_circ_map); + AdvectionPhysicalDomain const + physical_czarny_mapping(from_czarny_map, to_czarny_map); + AdvectionPseudoCartesianDomain const pseudo_cartesian_czarny_mapping( + from_czarny_map); std::tuple simulations = std::make_tuple( SimulationParameters( - circ_map, - circ_map, + from_circ_map, + to_circ_map, + from_circ_map, physical_circular_mapping, "CIRCULAR", "PHYSICAL"), SimulationParameters( - czarny_map, - czarny_map, + from_czarny_map, + to_czarny_map, + from_czarny_map, physical_czarny_mapping, "CZARNY", "PHYSICAL"), SimulationParameters( - czarny_map, - czarny_map, + from_czarny_map, + to_czarny_map, + from_czarny_map, pseudo_cartesian_czarny_mapping, "CZARNY", "PSEUDO CARTESIAN"), SimulationParameters( discrete_czarny_map, - czarny_map, + to_czarny_map, + from_czarny_map, pseudo_cartesian_czarny_mapping, "DISCRETE", "PSEUDO CARTESIAN")); diff --git a/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp b/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp index 9ae606e05..d9c12a841 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp +++ b/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp @@ -8,6 +8,8 @@ #include +#include +#include #include #include #include @@ -40,15 +42,23 @@ namespace fs = std::filesystem; namespace { #if defined(CIRCULAR_MAPPING_PHYSICAL) -using X_adv = typename AdvectionPhysicalDomain>::X_adv; -using Y_adv = typename AdvectionPhysicalDomain>::Y_adv; +using X_adv = typename AdvectionPhysicalDomain< + CircularToCartesian, + CartesianToCircular>::X_adv; +using Y_adv = typename AdvectionPhysicalDomain< + CircularToCartesian, + CartesianToCircular>::Y_adv; #elif defined(CZARNY_MAPPING_PHYSICAL) -using X_adv = typename AdvectionPhysicalDomain>::X_adv; -using Y_adv = typename AdvectionPhysicalDomain>::Y_adv; +using X_adv = typename AdvectionPhysicalDomain< + CzarnyToCartesian, + CartesianToCzarny>::X_adv; +using Y_adv = typename AdvectionPhysicalDomain< + CzarnyToCartesian, + CartesianToCzarny>::Y_adv; #elif defined(CZARNY_MAPPING_PSEUDO_CARTESIAN) -using X_adv = typename AdvectionPseudoCartesianDomain>::X_adv; -using Y_adv = typename AdvectionPseudoCartesianDomain>::Y_adv; +using X_adv = typename AdvectionPseudoCartesianDomain>::X_adv; +using Y_adv = typename AdvectionPseudoCartesianDomain>::Y_adv; #elif defined(DISCRETE_MAPPING_PSEUDO_CARTESIAN) using X_adv = typename AdvectionPseudoCartesianDomain< @@ -132,7 +142,7 @@ int main(int argc, char** argv) // DEFINITION OF OPERATORS ------------------------------------------------------------------ - // --- Builders for the test function and the mapping: + // --- Builders for the test function and the to_physical_mapping: SplineRThetaBuilder const builder(grid); // --- Evaluator for the test function: @@ -162,9 +172,10 @@ int main(int argc, char** argv) // SELECTION OF THE MAPPING AND THE ADVECTION DOMAIN. #if defined(CIRCULAR_MAPPING_PHYSICAL) - CircularToCartesian analytical_mapping; - CircularToCartesian mapping; - AdvectionPhysicalDomain advection_domain(analytical_mapping); + CircularToCartesian analytical_mapping; + CircularToCartesian to_physical_mapping; + CartesianToCircular to_logical_mapping; + AdvectionPhysicalDomain advection_domain(analytical_mapping, to_logical_mapping); std::string const mapping_name = "CIRCULAR"; std::string const adv_domain_name = "PHYSICAL"; key += "circular_physical"; @@ -174,31 +185,34 @@ int main(int argc, char** argv) double const czarny_epsilon = 1.4; #if defined(CZARNY_MAPPING_PHYSICAL) - CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); - CzarnyToCartesian mapping(czarny_e, czarny_epsilon); - AdvectionPhysicalDomain advection_domain(analytical_mapping); + CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); + CzarnyToCartesian to_physical_mapping(czarny_e, czarny_epsilon); + CartesianToCzarny to_logical_mapping(czarny_e, czarny_epsilon); + AdvectionPhysicalDomain advection_domain(analytical_mapping, to_logical_mapping); std::string const mapping_name = "CZARNY"; std::string const adv_domain_name = "PHYSICAL"; key += "czarny_physical"; #elif defined(CZARNY_MAPPING_PSEUDO_CARTESIAN) - CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); - CzarnyToCartesian mapping(czarny_e, czarny_epsilon); - AdvectionPseudoCartesianDomain advection_domain(mapping); + CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); + CzarnyToCartesian to_physical_mapping(czarny_e, czarny_epsilon); + CartesianToCzarny to_logical_mapping(czarny_e, czarny_epsilon); + AdvectionPseudoCartesianDomain advection_domain(to_physical_mapping); std::string const mapping_name = "CZARNY"; std::string const adv_domain_name = "PSEUDO CARTESIAN"; key += "czarny_pseudo_cartesian"; #elif defined(DISCRETE_MAPPING_PSEUDO_CARTESIAN) - CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); + CzarnyToCartesian analytical_mapping(czarny_e, czarny_epsilon); + CartesianToCzarny to_logical_mapping(czarny_e, czarny_epsilon); DiscreteToCartesianBuilder mapping_builder( Kokkos::DefaultHostExecutionSpace(), analytical_mapping, builder, spline_evaluator_extrapol); - DiscreteToCartesian mapping = mapping_builder(); - AdvectionPseudoCartesianDomain advection_domain(mapping); + DiscreteToCartesian to_physical_mapping = mapping_builder(); + AdvectionPseudoCartesianDomain advection_domain(to_physical_mapping); std::string const mapping_name = "DISCRETE"; std::string const adv_domain_name = "PSEUDO CARTESIAN"; key += "discrete_pseudo_cartesian"; @@ -247,17 +261,17 @@ int main(int argc, char** argv) // SELECTION OF THE SIMULATION. #if defined(TRANSLATION_SIMULATION) - TranslationSimulation simulation(mapping, rmin, rmax); + TranslationSimulation simulation(to_physical_mapping, rmin, rmax); std::string const simu_type = "TRANSLATION"; key += "Translation"; #elif defined(ROTATION_SIMULATION) - RotationSimulation simulation(mapping, rmin, rmax); + RotationSimulation simulation(to_physical_mapping, rmin, rmax); std::string const simu_type = "ROTATION"; key += "Rotation"; #elif defined(DECENTRED_ROTATION_SIMULATION) - DecentredRotationSimulation simulation(mapping); + DecentredRotationSimulation simulation(to_physical_mapping); std::string const simu_type = "DECENTRED ROTATION"; key += "Decentred_rotation"; #endif @@ -270,7 +284,8 @@ int main(int argc, char** argv) std::cout << mapping_name << " MAPPING - " << adv_domain_name << " DOMAIN - " << method_name << " - " << simu_type << " : " << std::endl; simulate( - mapping, + to_physical_mapping, + to_logical_mapping, analytical_mapping, grid, time_stepper, diff --git a/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp b/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp index b110306f9..d1432b80b 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp +++ b/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp @@ -49,22 +49,22 @@ std::string to_lower(std::string s) * The stream to which the output is printed. * @param[in] coord_rp * The coordinate to be printed. - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical index range. * @param[in] idx_range_p * The index range to which the poloidal coordinate should be restricted. */ -template +template void print_coordinate( std::ofstream& out_file, CoordRTheta coord_rp, - Mapping const& mapping, + LogicalToPhysicalMapping const& to_physical_mapping, IdxRangeTheta idx_range_p) { double const r = ddc::get(coord_rp); double const th = ddcHelper::restrict_to_idx_range(ddc::select(coord_rp), idx_range_p); - CoordXY coord_xy(mapping(coord_rp)); + CoordXY coord_xy(to_physical_mapping(coord_rp)); double const x = ddc::get(coord_xy); double const y = ddc::get(coord_xy); @@ -76,7 +76,7 @@ void print_coordinate( * @brief Save the characteristic feet in the logical index range * and the physical index range. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] rp_index range @@ -86,9 +86,9 @@ void print_coordinate( * @param[in] name * The name of the file where the feet are saved. */ -template +template void save_feet( - Mapping const& mapping, + LogicalToPhysicalMapping const& to_physical_mapping, IdxRangeRTheta const& idx_range_rp, host_t> const& feet_coords_rp, std::string const& name) @@ -106,8 +106,8 @@ void save_feet( IdxTheta ip(irp); file_feet << std::setw(15) << (ir - ir_start).value() << std::setw(15) << (ip - ip_start).value(); - print_coordinate(file_feet, ddc::coordinate(irp), mapping, idx_range_p); - print_coordinate(file_feet, feet_coords_rp(irp), mapping, idx_range_p); + print_coordinate(file_feet, ddc::coordinate(irp), to_physical_mapping, idx_range_p); + print_coordinate(file_feet, feet_coords_rp(irp), to_physical_mapping, idx_range_p); file_feet << std::endl; }); file_feet.close(); @@ -117,7 +117,7 @@ void save_feet( /** * @brief Save the advected function. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] function @@ -125,8 +125,11 @@ void save_feet( * @param[in] name * The name of the file where the feet are saved. */ -template -void saving_computed(Mapping const& mapping, host_t function, std::string const& name) +template +void saving_computed( + LogicalToPhysicalMapping const& to_physical_mapping, + host_t function, + std::string const& name) { IdxRangeRTheta const grid = get_idx_range(function); std::ofstream out_file(name, std::ofstream::out); @@ -143,7 +146,7 @@ void saving_computed(Mapping const& mapping, host_t function, std: out_file << std::setw(15) << (ir - ir_start).value() << std::setw(15) << (ip - ip_start).value(); - print_coordinate(out_file, ddc::coordinate(irp), mapping, idx_range_p); + print_coordinate(out_file, ddc::coordinate(irp), to_physical_mapping, idx_range_p); out_file << std::setw(25) << function(irp); out_file << std::endl; }); @@ -156,7 +159,7 @@ void saving_computed(Mapping const& mapping, host_t function, std: * * @param[in] idx_range_rp * The logical index range where the characteristic feet are defined. - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] advection_field @@ -167,28 +170,30 @@ void saving_computed(Mapping const& mapping, host_t function, std: * * @return A FieldMem with the exact characteristic feet at the given time. */ -template +template host_t> compute_exact_feet_rp( IdxRangeRTheta const& idx_range_rp, - Mapping const& mapping, + LogicalToPhysicalMapping const& logical_to_physical_mapping, + PhysicalToLogicalMapping const& physical_to_logical_mapping, AdvectionField const& advection_field, double const time) { - static_assert( - !std::is_same_v>); + static_assert(!std::is_same_v< + LogicalToPhysicalMapping, + DiscreteToCartesian>); host_t> feet_coords_rp(idx_range_rp); - CoordXY const coord_xy_center = CoordXY(mapping(CoordRTheta(0, 0))); - + CoordXY const coord_xy_center = CoordXY(logical_to_physical_mapping(CoordRTheta(0, 0))); ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) { CoordRTheta const coord_rp = ddc::coordinate(irp); - CoordXY const coord_xy = advection_field.exact_feet(mapping(coord_rp), time); + CoordXY const coord_xy + = advection_field.exact_feet(logical_to_physical_mapping(coord_rp), time); CoordXY const coord_diff = coord_xy - coord_xy_center; if (norm_inf(coord_diff) < 1e-15) { feet_coords_rp(irp) = CoordRTheta(0, 0); } else { - feet_coords_rp(irp) = mapping(coord_xy); + feet_coords_rp(irp) = physical_to_logical_mapping(coord_xy); } }); @@ -201,7 +206,7 @@ host_t> compute_exact_feet_rp( * the computed advected function and the exact * solution. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] grid @@ -216,9 +221,9 @@ host_t> compute_exact_feet_rp( * @return The L2 norm of the difference between * the computed function and the exact solution. */ -template +template double compute_difference_L2_norm( - Mapping const& mapping, + LogicalToPhysicalMapping const& to_physical_mapping, IdxRangeRTheta const& grid, host_t allfdistribu_advected, Function& function_to_be_advected, @@ -233,7 +238,7 @@ double compute_difference_L2_norm( host_t const quadrature_coeffs = compute_coeffs_on_mapping( Kokkos::DefaultHostExecutionSpace(), - mapping, + to_physical_mapping, trapezoid_quadrature_coefficients(grid)); host_t> quadrature(get_const_field(quadrature_coeffs)); @@ -276,7 +281,7 @@ void display_time( /** * @brief Run an advection simulation. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] analytical_mapping @@ -313,7 +318,7 @@ void display_time( * @param[in] counter_function * A integer referring to a test case for the name of the saved files. * - * @tparam Mapping + * @tparam LogicalToPhysicalMapping * A child class of CurvilinearToCartesian. * @tparam AnalyticalMapping * A child class of AnalyticalInvertibleCurvilinearToCartesian. @@ -330,13 +335,15 @@ void display_time( * @see Simulation */ template < - class Mapping, + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, class AnalyticalMapping, class TimeStepper, class AdvectionDomain, class Simulation> void simulate( - Mapping const& mapping, + LogicalToPhysicalMapping const& to_physical_mapping, + PhysicalToLogicalMapping const& to_logical_mapping, AnalyticalMapping const& analytical_mapping, IdxRangeRTheta const& grid, TimeStepper const& time_stepper, @@ -352,14 +359,14 @@ void simulate( bool if_save_feet, std::string const& output_folder) { - SplineFootFinder const foot_finder( + SplineFootFinder const foot_finder( time_stepper, advection_domain, - mapping, + to_physical_mapping, advection_builder, advection_evaluator); - BslAdvectionRTheta advection_operator(function_interpolator, foot_finder, mapping); + BslAdvectionRTheta advection_operator(function_interpolator, foot_finder, to_physical_mapping); auto function_to_be_advected_test = simulation.get_test_function(); auto advection_field_test = simulation.get_advection_field(); @@ -399,7 +406,7 @@ void simulate( // Definition of advection field: ddc::for_each(grid, [&](IdxRTheta const irp) { // Moving the coordinates in the physical index range: - CoordXY const coord_xy = mapping(ddc::coordinate(irp)); + CoordXY const coord_xy = to_physical_mapping(ddc::coordinate(irp)); CoordXY const advection_field = advection_field_test(coord_xy, 0.); // Define the advection field on the physical index range: @@ -419,7 +426,7 @@ void simulate( // Save the advected function for each iteration: if (if_save_curves) { std::string const name = output_folder + "/after_" + std::to_string(i + 1) + ".txt"; - saving_computed(mapping, get_field(allfdistribu_advected_test), name); + saving_computed(to_physical_mapping, get_field(allfdistribu_advected_test), name); } } @@ -429,9 +436,18 @@ void simulate( // Compute the exact characteristic feet: host_t> feet_coords_rp_end_time(grid); host_t> feet_coords_rp_dt(grid); - feet_coords_rp_end_time - = compute_exact_feet_rp(grid, analytical_mapping, advection_field_test, end_time); - feet_coords_rp_dt = compute_exact_feet_rp(grid, analytical_mapping, advection_field_test, dt); + feet_coords_rp_end_time = compute_exact_feet_rp( + grid, + analytical_mapping, + to_logical_mapping, + advection_field_test, + end_time); + feet_coords_rp_dt = compute_exact_feet_rp( + grid, + analytical_mapping, + to_logical_mapping, + advection_field_test, + dt); // Compute the maximal absolute error on the space at the end of the simulation: @@ -452,7 +468,7 @@ void simulate( std::cout << " ... " << "Relative L2 norm error: " << compute_difference_L2_norm( - mapping, + to_physical_mapping, grid, allfdistribu_advected_test, function_to_be_advected_test, @@ -474,7 +490,7 @@ void simulate( ddc::for_each(grid, [&](const IdxRTheta irp) { feet(irp) = ddc::coordinate(irp); }); foot_finder(get_field(feet), get_const_field(advection_field_test_vec), dt); std::string const name = output_folder + "/feet_computed.txt"; - save_feet(mapping, grid, get_field(feet), name); + save_feet(to_physical_mapping, grid, get_field(feet), name); } // Save the values of the exact function at the initial and final states: @@ -491,8 +507,8 @@ void simulate( // Exact final state end_function(irp) = function_to_be_advected_test(feet_coords_rp_end_time(irp)); }); - saving_computed(mapping, get_field(initial_function), name_0); - saving_computed(mapping, get_field(end_function), name_1); + saving_computed(to_physical_mapping, get_field(initial_function), name_0); + saving_computed(to_physical_mapping, get_field(end_function), name_1); } @@ -500,7 +516,7 @@ void simulate( // Save the exact characteristic feet for a displacement on dt: if (if_save_feet) { std::string const name = output_folder + "/feet_exact.txt"; - save_feet(mapping, grid, get_field(feet_coords_rp_dt), name); + save_feet(to_physical_mapping, grid, get_field(feet_coords_rp_dt), name); } @@ -512,7 +528,7 @@ void simulate( /** * @brief Run three advection simulations for each test cases in the Simulation class. * - * @param[in] mapping + * @param[in] to_physical_mapping * The mapping function from the logical index range to the physical * index range. * @param[in] analytical_mapping @@ -554,9 +570,15 @@ void simulate( * @see AdvectionDomain * @see Simulation */ -template +template < + class LogicalToPhysicalMapping, + class PhysicalToLogicalMapping, + class AnalyticalMapping, + class TimeStepper, + class AdvectionDomain> void simulate_the_3_simulations( - Mapping const& mapping, + LogicalToPhysicalMapping const& to_physical_mapping, + PhysicalToLogicalMapping const& to_logical_mapping, AnalyticalMapping const& analytical_mapping, IdxRangeRTheta const& grid, TimeStepper& time_stepper, @@ -576,9 +598,9 @@ void simulate_the_3_simulations( double const rmin = ddc::coordinate(r_idx_range.front()); double const rmax = ddc::coordinate(r_idx_range.back()); - TranslationSimulation simulation_1(mapping, rmin, rmax); - RotationSimulation simulation_2(mapping, rmin, rmax); - DecentredRotationSimulation simulation_3(mapping); + TranslationSimulation simulation_1(to_physical_mapping, rmin, rmax); + RotationSimulation simulation_2(to_physical_mapping, rmin, rmax); + DecentredRotationSimulation simulation_3(to_physical_mapping); std::string const title_simu_1 = " TRANSLATION : "; std::string const title_simu_2 = " ROTATION : "; @@ -590,7 +612,8 @@ void simulate_the_3_simulations( } std::cout << title + title_simu_1 << std::endl; simulate( - mapping, + to_physical_mapping, + to_logical_mapping, analytical_mapping, grid, time_stepper, @@ -611,7 +634,8 @@ void simulate_the_3_simulations( } std::cout << title + title_simu_2 << std::endl; simulate( - mapping, + to_physical_mapping, + to_logical_mapping, analytical_mapping, grid, time_stepper, @@ -632,7 +656,8 @@ void simulate_the_3_simulations( } std::cout << title + title_simu_3 << std::endl; simulate( - mapping, + to_physical_mapping, + to_logical_mapping, analytical_mapping, grid, time_stepper, diff --git a/tests/geometryRTheta/advection_2d_rp/test_cases.hpp b/tests/geometryRTheta/advection_2d_rp/test_cases.hpp index 9dd3c621a..3951a8fdc 100644 --- a/tests/geometryRTheta/advection_2d_rp/test_cases.hpp +++ b/tests/geometryRTheta/advection_2d_rp/test_cases.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -347,7 +348,9 @@ class AdvectionField_rotation : public AdvectionField private: double const m_vr; double const m_vtheta; - CircularToCartesian const m_mapping; + CartesianToCircular const m_physical_to_logical_mapping; + CircularToCartesian const m_logical_to_physical_mapping; + public: /** @@ -361,15 +364,16 @@ class AdvectionField_rotation : public AdvectionField AdvectionField_rotation(CoordVr vr, CoordVtheta vtheta) : m_vr(vr) , m_vtheta(vtheta) - , m_mapping() + , m_physical_to_logical_mapping() + , m_logical_to_physical_mapping() { } CoordXY operator()(CoordXY const coord, double const t) const final { - CoordRTheta const coord_rp(m_mapping(coord)); + CoordRTheta const coord_rp(m_physical_to_logical_mapping(coord)); std::array, 2> jacobian; - m_mapping.jacobian_matrix(coord_rp, jacobian); + m_logical_to_physical_mapping.jacobian_matrix(coord_rp, jacobian); double const vx = m_vr * jacobian[0][0] + m_vtheta * jacobian[0][1]; double const vy = m_vr * jacobian[1][0] + m_vtheta * jacobian[1][1]; return CoordXY(vx, vy); @@ -377,9 +381,9 @@ class AdvectionField_rotation : public AdvectionField CoordXY exact_feet(CoordXY coord_xy, double const t) const final { - CoordRTheta const coord_rp(m_mapping(coord_xy)); + CoordRTheta const coord_rp(m_physical_to_logical_mapping(coord_xy)); CoordRTheta const velocity(m_vr, m_vtheta); - return m_mapping(coord_rp - t * velocity); + return m_logical_to_physical_mapping(coord_rp - t * velocity); } }; diff --git a/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp b/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp index 971850210..fd649c9cd 100644 --- a/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp +++ b/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp @@ -44,7 +44,8 @@ namespace { using DiscreteMappingBuilder = DiscreteToCartesianBuilder; -using Mapping = CircularToCartesian; +using LogicalToPhysicalMapping = CircularToCartesian; +using PhysicalToLogicalMapping = CartesianToCircular; namespace fs = std::filesystem; @@ -113,11 +114,12 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) PolarSplineEvaluator polar_spline_evaluator(r_extrapolation_rule); - // --- Define the mapping. ------------------------------------------------------------------------ - const Mapping mapping; + // --- Define the to_physical_mapping. ------------------------------------------------------------------------ + const LogicalToPhysicalMapping to_physical_mapping; + const PhysicalToLogicalMapping to_logical_mapping; DiscreteMappingBuilder const discrete_mapping_builder( Kokkos::DefaultHostExecutionSpace(), - mapping, + to_physical_mapping, builder, spline_evaluator_extrapol); DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); @@ -135,7 +137,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator); - AdvectionPhysicalDomain advection_idx_range(mapping); + AdvectionPhysicalDomain advection_idx_range(to_physical_mapping, to_logical_mapping); RK3>, host_t>, @@ -143,23 +145,23 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) SplineFootFinder find_feet( time_stepper, advection_idx_range, - mapping, + to_physical_mapping, builder, spline_evaluator_extrapol); - BslAdvectionRTheta advection_operator(interpolator, find_feet, mapping); + BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping); // --- Advection field finder --------------------------------------------------------------------- - AdvectionFieldFinder advection_field_computer(mapping); + AdvectionFieldFinder advection_field_computer(to_physical_mapping); // --- Choice of the simulation ------------------------------------------------------------------- #if defined(TRANSLATION) - TranslationAdvectionFieldSimulation simulation(mapping, rmin, rmax); + TranslationAdvectionFieldSimulation simulation(to_physical_mapping, rmin, rmax); #elif defined(ROTATION) - RotationAdvectionFieldSimulation simulation(mapping, rmin, rmax); + RotationAdvectionFieldSimulation simulation(to_physical_mapping, rmin, rmax); #elif defined(DECENTRED_ROTATION) - DecentredRotationAdvectionFieldSimulation simulation(mapping); + DecentredRotationAdvectionFieldSimulation simulation(to_physical_mapping); #endif // ================================================================================================ @@ -200,7 +202,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) auto advection_field = simulation.get_advection_field(); ddc::for_each(grid, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(ddc::coordinate(irp)); - CoordXY const coord_xy(mapping(coord_rp)); + CoordXY const coord_xy(to_physical_mapping(coord_rp)); allfdistribu_rp(irp) = function(coord_rp); allfdistribu_xy(irp) = allfdistribu_rp(irp); @@ -236,12 +238,12 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) // > Compare the advection field computed on RTheta to the advection field computed on XY host_t> difference_between_fields_xy_and_rp(grid); - MetricTensor metric_tensor(mapping); + MetricTensor metric_tensor(to_physical_mapping); ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(ddc::coordinate(irp)); std::array, 2> J; // Jacobian matrix - mapping.jacobian_matrix(coord_rp, J); + to_physical_mapping.jacobian_matrix(coord_rp, J); std::array, 2> G; // Metric tensor metric_tensor(G, coord_rp); diff --git a/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp b/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp index b2672eff2..f9ceca908 100644 --- a/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp +++ b/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp @@ -193,7 +193,7 @@ class ElectrostaticalPotentialSimulation_rotation : public ElectrostaticalPotent private: double const m_vr; double const m_vtheta; - CircularToCartesian const m_mapping; + CartesianToCircular const m_mapping; public: /** @@ -220,7 +220,8 @@ class ElectrostaticalPotentialSimulation_rotation : public ElectrostaticalPotent { CoordRTheta const coord_rp(m_mapping(coord_xy)); CoordRTheta const velocity(m_vr, m_vtheta); - return m_mapping(coord_rp - t * velocity); + CircularToCartesian logical_to_physical_mapping; + return logical_to_physical_mapping(coord_rp - t * velocity); } }; diff --git a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp index dcbed4955..9d60de5ee 100644 --- a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp +++ b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp @@ -29,9 +29,9 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver< SplineRThetaEvaluatorNullBound>; #if defined(CIRCULAR_MAPPING) -using Mapping = CircularToCartesian; +using Mapping = CircularToCartesian; #elif defined(CZARNY_MAPPING) -using Mapping = CzarnyToCartesian; +using Mapping = CzarnyToCartesian; #endif using DiscreteMappingBuilder = DiscreteToCartesianBuilder; diff --git a/tests/geometryRTheta/polar_poisson/test_cases.cpp b/tests/geometryRTheta/polar_poisson/test_cases.cpp index 17e787526..fc0d2c8de 100644 --- a/tests/geometryRTheta/polar_poisson/test_cases.cpp +++ b/tests/geometryRTheta/polar_poisson/test_cases.cpp @@ -2,14 +2,14 @@ #include "test_cases.hpp" template <> -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: solution_at_pole(Coord const& coord) const { return 0.0; } template <> -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: non_singular_solution(Coord const& coord) const { const double r = ddc::get(coord); @@ -75,14 +75,14 @@ double ManufacturedPoissonTest -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: solution_at_pole(Coord const& coord) const { return 0.0; } template <> -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: non_singular_solution(Coord const& coord) const { const double r = ddc::get(coord); @@ -1723,14 +1723,14 @@ double ManufacturedPoissonTest -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: solution_at_pole(Coord const& coord) const { return 0.0; } template <> -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: non_singular_solution(Coord const& coord) const { const double r = ddc::get(coord); @@ -1763,14 +1763,14 @@ double ManufacturedPoissonTest -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: solution_at_pole(Coord const& coord) const { return 0.0; } template <> -double ManufacturedPoissonTest>>:: +double ManufacturedPoissonTest>>:: non_singular_solution(Coord const& coord) const { const double r = ddc::get(coord); @@ -3315,7 +3315,7 @@ double ManufacturedPoissonTest>>; -template class ManufacturedPoissonTest>>; -template class ManufacturedPoissonTest>>; -template class ManufacturedPoissonTest>>; +template class ManufacturedPoissonTest>>; +template class ManufacturedPoissonTest>>; +template class ManufacturedPoissonTest>>; +template class ManufacturedPoissonTest>>; diff --git a/tests/geometryRTheta/polar_poisson/test_cases.hpp b/tests/geometryRTheta/polar_poisson/test_cases.hpp index d74d5d098..d8b85582f 100644 --- a/tests/geometryRTheta/polar_poisson/test_cases.hpp +++ b/tests/geometryRTheta/polar_poisson/test_cases.hpp @@ -6,6 +6,7 @@ #include #include +#include #include #include "ddc_aliases.hpp" diff --git a/tests/geometryRTheta/quadrature/tests_L1_and_L2_norms.cpp b/tests/geometryRTheta/quadrature/tests_L1_and_L2_norms.cpp index 6204cf747..e44debf80 100644 --- a/tests/geometryRTheta/quadrature/tests_L1_and_L2_norms.cpp +++ b/tests/geometryRTheta/quadrature/tests_L1_and_L2_norms.cpp @@ -227,7 +227,7 @@ TEST_P(SplineQuadrature, TestFunctions) // ------------------------------------------------------------------------------------------------ std::cout << "CIRCULAR MAPPING ---------------------------------------------------" << std::endl; - const CircularToCartesian mapping_1; + const CircularToCartesian mapping_1; std::array, 5> expected_norms; expected_norms[0][0] = M_PI; expected_norms[0][1] = std::sqrt(M_PI); diff --git a/tests/multipatch/connectivity/patch_locator_onion_shape_2patches.cpp b/tests/multipatch/connectivity/patch_locator_onion_shape_2patches.cpp index aefe9451e..7523ea493 100644 --- a/tests/multipatch/connectivity/patch_locator_onion_shape_2patches.cpp +++ b/tests/multipatch/connectivity/patch_locator_onion_shape_2patches.cpp @@ -2,6 +2,8 @@ #include +#include +#include #include #include @@ -32,7 +34,8 @@ class OnionPatchLocator2PatchesTest : public ::testing::Test public: using PatchLocator = OnionPatchLocator< MultipatchType, - CircularToCartesian>; + CircularToCartesian, + CartesianToCircular>; Patch1::IdxRange1 const idx_range_r1; Patch1::IdxRange2 const idx_range_theta1; @@ -102,9 +105,10 @@ TEST_F(OnionPatchLocator2PatchesTest, CheckPatchOrderingTest) Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_1, idx_range_2); - CzarnyToCartesian mapping(0.3, 1.4); + CzarnyToCartesian to_physical_mapping(0.3, 1.4); + CartesianToCzarny to_logical_mapping(0.3, 1.4); - EXPECT_NO_THROW(OnionPatchLocator(all_idx_ranges, mapping)); + EXPECT_NO_THROW((OnionPatchLocator(all_idx_ranges, to_physical_mapping, to_logical_mapping))); } TEST_F(OnionPatchLocator2PatchesTest, CheckPatchOrderingDeathTest) @@ -113,22 +117,36 @@ TEST_F(OnionPatchLocator2PatchesTest, CheckPatchOrderingDeathTest) Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_2, idx_range_1); - CircularToCartesian mapping; - - EXPECT_THROW(OnionPatchLocator(all_idx_ranges, mapping), std::invalid_argument); + CircularToCartesian to_physical_mapping; + CartesianToCircular to_logical_mapping; + + EXPECT_THROW( + (OnionPatchLocator< + MultipatchType, + CircularToCartesian, + CartesianToCircular, + Kokkos::DefaultExecutionSpace>( + all_idx_ranges, + to_physical_mapping, + to_logical_mapping)), + std::invalid_argument); } - - TEST_F(OnionPatchLocator2PatchesTest, DeviceCircularOnionPatchLocator2PatchesTest) { Patch1::IdxRange12 idx_range_1(idx_range_r1, idx_range_theta1); Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_1, idx_range_2); - CircularToCartesian mapping; + CircularToCartesian to_physical_mapping; + CartesianToCircular to_logical_mapping; - OnionPatchLocator patch_locator(all_idx_ranges, mapping); + OnionPatchLocator< + MultipatchType, + CircularToCartesian, + CartesianToCircular, + Kokkos::DefaultExecutionSpace> + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); int constexpr n_elements = 7; Kokkos::View*> coords("coords 1", n_elements); @@ -168,9 +186,15 @@ TEST_F(OnionPatchLocator2PatchesTest, DeviceCzarnyOnionPatchLocator2PatchesTest) Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_1, idx_range_2); - CzarnyToCartesian mapping(0.3, 1.4); + CzarnyToCartesian to_physical_mapping(0.3, 1.4); + CartesianToCzarny to_logical_mapping(0.3, 1.4); - OnionPatchLocator patch_locator(all_idx_ranges, mapping); + OnionPatchLocator< + MultipatchType, + CzarnyToCartesian, + CartesianToCzarny, + Kokkos::DefaultExecutionSpace> + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); int constexpr n_elements = 7; Kokkos::View*> coords("coords 1", n_elements); @@ -181,12 +205,12 @@ TEST_F(OnionPatchLocator2PatchesTest, DeviceCzarnyOnionPatchLocator2PatchesTest) Kokkos::View patches_host("patches_host 1", n_elements); - coords_host(0) = PhysicalCoordXY(mapping(Coord(0, 0))); - coords_host(1) = PhysicalCoordXY(mapping(Coord(0.2, 0))); + coords_host(0) = PhysicalCoordXY(to_physical_mapping(Coord(0, 0))); + coords_host(1) = PhysicalCoordXY(to_physical_mapping(Coord(0.2, 0))); coords_host(2) = PhysicalCoordXY(0.25, .03); coords_host(3) = PhysicalCoordXY(1, 1); - coords_host(4) = PhysicalCoordXY(mapping(Coord(1, 0))); - coords_host(5) = PhysicalCoordXY(mapping(Coord(1.5, 0))); + coords_host(4) = PhysicalCoordXY(to_physical_mapping(Coord(1, 0))); + coords_host(5) = PhysicalCoordXY(to_physical_mapping(Coord(1.5, 0))); coords_host(6) = PhysicalCoordXY(-2.1, 0); patches_host(0) = PatchLocator::outside_rmin_domain; @@ -210,13 +234,18 @@ TEST_F(OnionPatchLocator2PatchesTest, HostCircularOnionPatchLocator2PatchesTest) Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_1, idx_range_2); - CircularToCartesian mapping; - + CircularToCartesian to_physical_mapping; + CartesianToCircular to_logical_mapping; using MultipatchIdxRanges = MultipatchType; - using Mapping = CircularToCartesian; + using LogicalToPhysicalMapping = CircularToCartesian; + using PhysicalToLogicalMapping = CartesianToCircular; - OnionPatchLocator - patch_locator(all_idx_ranges, mapping); + OnionPatchLocator< + MultipatchIdxRanges, + LogicalToPhysicalMapping, + PhysicalToLogicalMapping, + Kokkos::DefaultHostExecutionSpace> + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); PhysicalCoordXY coord(0.15, .03); EXPECT_EQ(patch_locator(coord), PatchLocator::outside_rmin_domain); @@ -247,18 +276,23 @@ TEST_F(OnionPatchLocator2PatchesTest, HostCzarnyOnionPatchLocator2PatchesTest) Patch2::IdxRange12 idx_range_2(idx_range_r2, idx_range_theta2); MultipatchType all_idx_ranges(idx_range_1, idx_range_2); - CzarnyToCartesian mapping(0.3, 1.4); + CzarnyToCartesian to_physical_mapping(0.3, 1.4); + CartesianToCzarny to_logical_mapping(0.3, 1.4); using MultipatchIdxRanges = MultipatchType; - using Mapping = CzarnyToCartesian; - - OnionPatchLocator - patch_locator(all_idx_ranges, mapping); + using LogicalToPhysicalMapping = CzarnyToCartesian; + using PhysicalToLogicalMapping = CartesianToCzarny; + OnionPatchLocator< + MultipatchIdxRanges, + LogicalToPhysicalMapping, + PhysicalToLogicalMapping, + Kokkos::DefaultHostExecutionSpace> + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); PhysicalCoordXY coord(-0.05, .03); EXPECT_EQ(patch_locator(coord), PatchLocator::outside_rmin_domain); - coord = PhysicalCoordXY(mapping(Coord(0.2, 0))); + coord = PhysicalCoordXY(to_physical_mapping(Coord(0.2, 0))); EXPECT_EQ(patch_locator(coord), 0); coord = PhysicalCoordXY(0.25, .03); @@ -267,12 +301,12 @@ TEST_F(OnionPatchLocator2PatchesTest, HostCzarnyOnionPatchLocator2PatchesTest) coord = PhysicalCoordXY(1, 1); EXPECT_EQ(patch_locator(coord), 1); - coord = PhysicalCoordXY(mapping(Coord(1, 0))); + coord = PhysicalCoordXY(to_physical_mapping(Coord(1, 0))); EXPECT_EQ(patch_locator(coord), 1); - coord = PhysicalCoordXY(mapping(Coord(1.5, 0))); + coord = PhysicalCoordXY(to_physical_mapping(Coord(1.5, 0))); EXPECT_EQ(patch_locator(coord), 1); coord = PhysicalCoordXY(-2.1, 0); EXPECT_EQ(patch_locator(coord), PatchLocator::outside_rmax_domain); -} \ No newline at end of file +} diff --git a/tests/multipatch/spline/multipatch_spline_evaluator.cpp b/tests/multipatch/spline/multipatch_spline_evaluator.cpp index 14123243e..0d1368041 100644 --- a/tests/multipatch/spline/multipatch_spline_evaluator.cpp +++ b/tests/multipatch/spline/multipatch_spline_evaluator.cpp @@ -2,8 +2,6 @@ #include #include -#include - #include #include "2patches_2d_onion_shape_non_uniform.hpp" @@ -74,8 +72,8 @@ class MultipatchSplineEvaluatorTest : public MultipatchSplineOnionShapeTest DField, BSplinesTheta<1>>> const function_1_coef; DField, BSplinesTheta<2>>> const function_2_coef; - Mapping const mapping; - + LogicalToPhysicalMapping const to_physical_mapping; + PhysicalToLogicalMapping const to_logical_mapping; public: MultipatchSplineEvaluatorTest() @@ -96,7 +94,10 @@ class MultipatchSplineEvaluatorTest : public MultipatchSplineOnionShapeTest , evaluator_2_host(bc_r_min_2, bc_r_max_2, bc_theta, bc_theta) // Local splines on device and host , function_1_coef(get_field(function_1_coef_alloc)) - , function_2_coef(get_field(function_2_coef_alloc)) {}; + , function_2_coef(get_field(function_2_coef_alloc)) + , to_logical_mapping() + { + } template @@ -283,7 +284,8 @@ void test_deriv( TEST_F(MultipatchSplineEvaluatorTest, HostEvaluateOnSingleCoord) { // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); HostMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -344,7 +346,8 @@ TEST_F(MultipatchSplineEvaluatorTest, HostEvaluateOnSingleCoord) TEST_F(MultipatchSplineEvaluatorTest, DeviceEvaluateOnSingleCoord) { // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); DeviceMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -433,7 +436,8 @@ TEST_F(MultipatchSplineEvaluatorTest, EvaluateOnCoordField) // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); DeviceMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -484,7 +488,8 @@ TEST_F(MultipatchSplineEvaluatorTest, EvaluateOnCoordField) TEST_F(MultipatchSplineEvaluatorTest, HostDerivativesOnSingleCoord) { // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); HostMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -592,7 +597,8 @@ TEST_F(MultipatchSplineEvaluatorTest, HostDerivativesOnSingleCoord) TEST_F(MultipatchSplineEvaluatorTest, DeviceDerivativesOnSingleCoord) { // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); DeviceMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -677,7 +683,8 @@ TEST_F(MultipatchSplineEvaluatorTest, DeviceDerivativesOnSingleCoord) TEST_F(MultipatchSplineEvaluatorTest, DerivativesOnSingleCoordDeathTest) { // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); DeviceMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -798,7 +805,8 @@ TEST_F(MultipatchSplineEvaluatorTest, DerivativativesOnCoordField) // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); DeviceMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); @@ -904,7 +912,8 @@ TEST_F(MultipatchSplineEvaluatorTest, HostIntegrateOnCoordField) splines_host(get_field(function_1_coef_host), get_field(function_2_coef_host)); // Definition of MultipatchSplineEvaluator2D - PatchLocator const patch_locator(all_idx_ranges, mapping); + PatchLocator const + patch_locator(all_idx_ranges, to_physical_mapping, to_logical_mapping); ConstantExtrapolationRuleOnion> extrapolation_rule(r1_min, r2_max); HostMultipatchSplineRThetaEvaluator const evaluators(patch_locator, extrapolation_rule); diff --git a/tests/multipatch/spline/null_extrapolation_rules.cpp b/tests/multipatch/spline/null_extrapolation_rules.cpp index 584c8f355..f50abbe05 100644 --- a/tests/multipatch/spline/null_extrapolation_rules.cpp +++ b/tests/multipatch/spline/null_extrapolation_rules.cpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: MIT #include +#include #include #include diff --git a/tests/multipatch/spline/spline_testing_tools.hpp b/tests/multipatch/spline/spline_testing_tools.hpp index 66ac4afcd..d9b232ff8 100644 --- a/tests/multipatch/spline/spline_testing_tools.hpp +++ b/tests/multipatch/spline/spline_testing_tools.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -70,11 +71,16 @@ using SplineRThetaEvaluator = ddc::SplineEvaluator2D< GridTheta>; -using Mapping = CircularToCartesian; +using LogicalToPhysicalMapping = CircularToCartesian; +using PhysicalToLogicalMapping = CartesianToCircular; using MultipatchIdxRange = MultipatchType; template -using PatchLocator = OnionPatchLocator; +using PatchLocator = OnionPatchLocator< + MultipatchIdxRange, + LogicalToPhysicalMapping, + PhysicalToLogicalMapping, + ExecSpace>; /** diff --git a/vendor/sll/include/sll/mapping/cartesian_to_circular.hpp b/vendor/sll/include/sll/mapping/cartesian_to_circular.hpp new file mode 100644 index 000000000..c37fd46d6 --- /dev/null +++ b/vendor/sll/include/sll/mapping/cartesian_to_circular.hpp @@ -0,0 +1,233 @@ +// SPDX-License-Identifier: MIT +#pragma once + +#include + +#include + +#include "mapping_tools.hpp" + +/** + * @brief A class for describing the circular 2D mapping. + * + * The mapping @f$ (x,y)\mapsto (r,\theta) @f$ is defined as follow : + * + * @f$ r(x,y) = \sqrt x^2+y^2 ,@f$ + * + * @f$ \theta(x,y) = atan2(\frac{y}{x}) .@f$ + * + * It and its Jacobian matrix are invertible everywhere except for @f$ r = 0 @f$. + * + * The Jacobian matrix coefficients are defined as follow + * + * @f$ J_{11}(r,\theta) =\frac{2x}{\sqrt{x^2+y^2}} @f$ + * + * @f$ J_{12}(r,\theta) =\frac{2y}{\sqrt{x^2+y^2}} @f$ + * + * @f$ J_{21}(r,\theta) =\frac{-y}{(x^2+y^2)^2} @f$ + * + * @f$ J_{22}(r,\theta) =\frac{x}{(x^2+y^2)^2} @f$ + * + * and the matrix determinant: @f$ det(J) = r @f$. + * + */ +template +class CartesianToCircular +{ +public: + /// @brief Indicate the first physical coordinate. + using cartesian_tag_x = X; + /// @brief Indicate the second physical coordinate. + using cartesian_tag_y = Y; + /// @brief Indicate the first logical coordinate. + using curvilinear_tag_r = R; + /// @brief Indicate the second logical coordinate. + using curvilinear_tag_theta = Theta; + + /// The type of the argument of the function described by this mapping + using CoordArg = ddc::Coordinate; + /// The type of the result of the function described by this mapping + using CoordResult = ddc::Coordinate; + +public: + CartesianToCircular() = default; + + /** + * @brief Instantiate a CartesianToCircular from another CartesianToCircular (lvalue). + * + * @param[in] other + * CartesianToCircular mapping used to instantiate the new one. + */ + KOKKOS_FUNCTION CartesianToCircular(CartesianToCircular const& other) {} + + /** + * @brief Instantiate a Curvilinear2DToCartesian from another temporary CartesianToCircular (rvalue). + * + * @param[in] x + * Curvilinear2DToCartesian mapping used to instantiate the new one. + */ + CartesianToCircular(CartesianToCircular&& x) = default; + + ~CartesianToCircular() = default; + + /** + * @brief Assign a CartesianToCircular from another CartesianToCircular (lvalue). + * + * @param[in] x + * CartesianToCircular mapping used to assign. + * + * @return The CartesianToCircular assigned. + */ + CartesianToCircular& operator=(CartesianToCircular const& x) = default; + + /** + * @brief Assign a CartesianToCircular from another temporary CartesianToCircular (rvalue). + * + * @param[in] x + * CartesianToCircular mapping used to assign. + * + * @return The CartesianToCircular assigned. + */ + CartesianToCircular& operator=(CartesianToCircular&& x) = default; + + /** + * @brief Convert the coordinate (x,y) to the equivalent @f$ (r, \theta) @f$ coordinate. + * + * @param[in] coord The coordinate to be converted. + * + * @return The equivalent coordinate. + */ + KOKKOS_FUNCTION ddc::Coordinate operator()(ddc::Coordinate const& coord) const + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + const double r = Kokkos::sqrt(x * x + y * y); + const double theta = Kokkos::atan2(y, x); + return ddc::Coordinate(r, theta); + } + + /** + * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian. + * + * @return A double with the value of the determinant of the Jacobian matrix. + */ + KOKKOS_FUNCTION double jacobian(ddc::Coordinate const& coord) + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + return 2. * (x * x - y * y) / Kokkos::pow(x * x + y * y, 1.5); + } + + /** + * @brief Compute full Jacobian matrix. + * + * For some computations, we need the complete Jacobian matrix or just the + * coefficients. + * The coefficients can be given independently with the functions + * jacobian_11, jacobian_12, jacobian_21 and jacobian_22. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian matrix. + * @param[out] matrix + * The Jacobian matrix returned. + * + * + * @see Jacobian::jacobian_11 + * @see Jacobian::jacobian_12 + * @see Jacobian::jacobian_21 + * @see Jacobian::jacobian_22 + */ + KOKKOS_FUNCTION void jacobian_matrix(ddc::Coordinate const& coord, Matrix_2x2& matrix) + + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + matrix[0][0] = 2 * x / Kokkos::pow(x * x + y * y, 0.5); + matrix[0][1] = 2 * y / Kokkos::pow(x * x + y * y, 0.5); + matrix[1][0] = -y / Kokkos::pow(x * x + y * y, 2.); + matrix[1][1] = x / Kokkos::pow(x * x + y * y, 2.); + } + + /** + * @brief Compute the (1,1) coefficient of the Jacobian matrix. + * + * For a mapping given by @f$ \mathcal{F} : (x,y)\mapsto (r,\theta) @f$, the + * (1,1) coefficient of the Jacobian matrix is given by @f$ \frac{\partial r}{\partial x} @f$. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian matrix. + * + * @return A double with the value of the (1,1) coefficient of the Jacobian matrix. + */ + KOKKOS_FUNCTION double jacobian_11(ddc::Coordinate const& coord) + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + return 2 * x / Kokkos::pow(x * x + y * y, 0.5); + ; + } + + /** + * @brief Compute the (1,2) coefficient of the Jacobian matrix. + * + * For a mapping given by @f$ \mathcal{F} : (x,y)\mapsto (r,\theta) @f$, the + * (1,2) coefficient of the Jacobian matrix is given by @f$ \frac{\partial \theta}{\partial x} @f$. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian matrix. + * + * @return A double with the value of the (1,2) coefficient of the Jacobian matrix. + */ + KOKKOS_FUNCTION double jacobian_12(ddc::Coordinate const& coord) + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + return 2 * y / Kokkos::pow(x * x + y * y, 0.5); + } + + /** + * @brief Compute the (2,1) coefficient of the Jacobian matrix. + * + * For a mapping given by @f$ \mathcal{F} : (x,y)\mapsto (r,\theta) @f$, the + * (2,1) coefficient of the Jacobian matrix is given by @f$ \frac{\partial r}{\partial y} @f$. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian matrix. . + * + * @return A double with the value of the (2,1) coefficient of the Jacobian matrix. + */ + KOKKOS_FUNCTION double jacobian_21(ddc::Coordinate const& coord) + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + return -y / Kokkos::pow(x * x + y * y, 2.); + } + + /** + * @brief Compute the (2,2) coefficient of the Jacobian matrix. + * + * For a mapping given by @f$ \mathcal{F} : (x,y)\mapsto (r,\theta) @f$, the + * (2,2) coefficient of the Jacobian matrix is given by @f$ \frac{\partial \theta}{\partial y} @f$. + * + * @param[in] coord + * The coordinate where we evaluate the Jacobian matrix. + * + * @return A double with the value of the (2,2) coefficient of the Jacobian matrix. + */ + KOKKOS_FUNCTION double jacobian_22(ddc::Coordinate const& coord) + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + return x / Kokkos::pow(x * x + y * y, 2.); + } +}; + +namespace mapping_detail { +template +struct MappingAccessibility> : std::true_type +{ +}; +} // namespace mapping_detail diff --git a/vendor/sll/include/sll/mapping/cartesian_to_czarny.hpp b/vendor/sll/include/sll/mapping/cartesian_to_czarny.hpp new file mode 100644 index 000000000..d2f548c98 --- /dev/null +++ b/vendor/sll/include/sll/mapping/cartesian_to_czarny.hpp @@ -0,0 +1,153 @@ +// SPDX-License-Identifier: MIT +#pragma once + +#include + +#include + +#include "mapping_tools.hpp" + + +/** + * @brief A class for describing the Czarny 2D mapping. + * + * The mapping @f$ (x,y) \mapsto (r,\theta)@f$ is defined by + * + * @f$ r(x,y) = \sqrt{\frac{y^2 (1+\epsilon x)^2}{e^2\xi^2+0.25(\epsilon x^2-2x-\epsilon)^2}},@f$ + * + * @f$ \theta (x,y)) = atan2(2. y (1+\epsilon x), (e \xi (\epsilon x^2 - 2x-\epsilon))), @f$ + * + * with @f$ \xi = 1/\sqrt{1 - \epsilon^2 /4} @f$ and @f$ e @f$ and @f$ \epsilon @f$ given as parameters. + */ +template +class CartesianToCzarny +{ +public: + /// @brief Indicate the first physical coordinate. + using cartesian_tag_x = X; + /// @brief Indicate the second physical coordinate. + using cartesian_tag_y = Y; + /// @brief Indicate the first logical coordinate. + using curvilinear_tag_r = R; + /// @brief Indicate the second logical coordinate. + using curvilinear_tag_theta = Theta; + + /// The type of the argument of the function described by this mapping + using CoordArg = ddc::Coordinate; + /// The type of the result of the function described by this mapping + using CoordResult = ddc::Coordinate; + +private: + double m_epsilon; + double m_e; + +public: + /** + * @brief Instantiate a CartesianToCzarny from parameters. + * + * @param[in] epsilon + * The @f$ \epsilon @f$ parameter in the definition of the mapping CartesianToCzarny. + * + * @param[in] e + * The @f$ e @f$ parameter in the definition of the mapping CartesianToCzarny. + * + * @see CartesianToCzarny + */ + CartesianToCzarny(double epsilon, double e) : m_epsilon(epsilon), m_e(e) {} + + /** + * @brief Instantiate a CartesianToCzarny from another CartesianToCzarny (lvalue). + * + * @param[in] other + * CartesianToCzarny mapping used to instantiate the new one. + */ + KOKKOS_FUNCTION CartesianToCzarny(CartesianToCzarny const& other) + : m_epsilon(other.epsilon()) + , m_e(other.e()) + { + } + + /** + * @brief Instantiate a CartesianToCzarny from another temporary CartesianToCzarny (rvalue). + * + * @param[in] x + * CartesianToCzarny mapping used to instantiate the new one. + */ + CartesianToCzarny(CartesianToCzarny&& x) = default; + + ~CartesianToCzarny() = default; + + /** + * @brief Assign a CartesianToCzarny from another CartesianToCzarny (lvalue). + * + * @param[in] x + * CartesianToCzarny mapping used to assign. + * + * @return The CartesianToCzarny assigned. + */ + CartesianToCzarny& operator=(CartesianToCzarny const& x) = default; + + /** + * @brief Assign a CartesianToCzarny from another temporary CartesianToCzarny (rvalue). + * + * @param[in] x + * CartesianToCzarny mapping used to assign. + * + * @return The CartesianToCzarny assigned. + */ + CartesianToCzarny& operator=(CartesianToCzarny&& x) = default; + + /** + * @brief Return the @f$ \epsilon @f$ parameter. + * + * @return The value of @f$ \epsilon @f$. + * + * @see CartesianToCzarny + */ + KOKKOS_FUNCTION double epsilon() const + { + return m_epsilon; + } + + /** + * @brief Return the @f$ e @f$ parameter. + * + * @return The value of @f$ e @f$. + * + * @see CartesianToCzarny + */ + KOKKOS_FUNCTION double e() const + { + return m_e; + } + + /** + * @brief Convert the coordinate (x,y) to the equivalent @f$ (r, \theta) @f$ coordinate. + * + * @param[in] coord The coordinate to be converted. + * + * @return The equivalent coordinate. + */ + KOKKOS_FUNCTION ddc::Coordinate operator()(ddc::Coordinate const& coord) const + { + const double x = ddc::get(coord); + const double y = ddc::get(coord); + const double ex = 1. + m_epsilon * x; + const double ex2 = (m_epsilon * x * x - 2. * x - m_epsilon); + const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25); + const double xi = Kokkos::sqrt(xi2); + const double r = Kokkos::sqrt(y * y * ex * ex / (m_e * m_e * xi2) + ex2 * ex2 * 0.25); + double theta = Kokkos::atan2(2. * y * ex, (m_e * xi * ex2)); + if (theta < 0) { + theta = 2 * M_PI + theta; + } + return ddc::Coordinate(r, theta); + } +}; + +namespace mapping_detail { +template +struct MappingAccessibility> : std::true_type +{ +}; +} // namespace mapping_detail diff --git a/vendor/sll/include/sll/mapping/circular_to_cartesian.hpp b/vendor/sll/include/sll/mapping/circular_to_cartesian.hpp index 0581e0061..ad5c9ace3 100644 --- a/vendor/sll/include/sll/mapping/circular_to_cartesian.hpp +++ b/vendor/sll/include/sll/mapping/circular_to_cartesian.hpp @@ -8,7 +8,6 @@ #include -#include "coordinate_converter.hpp" #include "curvilinear2d_to_cartesian.hpp" #include "mapping_tools.hpp" #include "pseudo_cartesian_compatible_mapping.hpp" @@ -37,10 +36,9 @@ * and the matrix determinant: @f$ det(J) = r @f$. * */ -template +template class CircularToCartesian - : public CoordinateConverter, ddc::Coordinate> - , public PseudoCartesianCompatibleMapping + : public PseudoCartesianCompatibleMapping , public Curvilinear2DToCartesian { public: @@ -116,23 +114,6 @@ class CircularToCartesian return ddc::Coordinate(x, y); } - /** - * @brief Convert the coordinate (x,y) to the equivalent @f$ (r, \theta) @f$ coordinate. - * - * @param[in] coord The coordinate to be converted. - * - * @return The equivalent coordinate. - */ - KOKKOS_FUNCTION ddc::Coordinate operator()( - ddc::Coordinate const& coord) const final - { - const double x = ddc::get(coord); - const double y = ddc::get(coord); - const double r = Kokkos::sqrt(x * x + y * y); - const double theta = Kokkos::atan2(y, x); - return ddc::Coordinate(r, theta); - } - /** * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping. * @@ -147,13 +128,12 @@ class CircularToCartesian return r; } - /** * @brief Compute full Jacobian matrix. * * For some computations, we need the complete Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * jacobian_11, jacobian_12, jacobian_21 and jacobian_22. * * @param[in] coord @@ -248,7 +228,7 @@ class CircularToCartesian * * For some computations, we need the complete inverse Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22. * * @param[in] coord @@ -431,10 +411,9 @@ class CircularToCartesian } }; - namespace mapping_detail { template -struct MappingAccessibility> : std::true_type +struct MappingAccessibility> : std::true_type { }; } // namespace mapping_detail diff --git a/vendor/sll/include/sll/mapping/czarny_to_cartesian.hpp b/vendor/sll/include/sll/mapping/czarny_to_cartesian.hpp index a5ea2a6a4..fb580bb32 100644 --- a/vendor/sll/include/sll/mapping/czarny_to_cartesian.hpp +++ b/vendor/sll/include/sll/mapping/czarny_to_cartesian.hpp @@ -7,7 +7,6 @@ #include -#include "coordinate_converter.hpp" #include "curvilinear2d_to_cartesian.hpp" #include "mapping_tools.hpp" #include "pseudo_cartesian_compatible_mapping.hpp" @@ -45,11 +44,8 @@ * \frac{e\xi}{2 - \sqrt{1 + \epsilon(\epsilon + 2 r \cos(\theta))}}. @f$ * */ -template -class CzarnyToCartesian - : public CoordinateConverter, ddc::Coordinate> - , public PseudoCartesianCompatibleMapping - , public Curvilinear2DToCartesian +template +class CzarnyToCartesian : public PseudoCartesianCompatibleMapping { public: /// @brief Indicate the first physical coordinate. @@ -172,31 +168,6 @@ class CzarnyToCartesian return ddc::Coordinate(x, y); } - /** - * @brief Convert the coordinate (x,y) to the equivalent @f$ (r, \theta) @f$ coordinate. - * - * @param[in] coord The coordinate to be converted. - * - * @return The equivalent coordinate. - */ - KOKKOS_FUNCTION ddc::Coordinate operator()( - ddc::Coordinate const& coord) const final - { - const double x = ddc::get(coord); - const double y = ddc::get(coord); - const double ex = 1. + m_epsilon * x; - const double ex2 = (m_epsilon * x * x - 2. * x - m_epsilon); - const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25); - const double xi = Kokkos::sqrt(xi2); - const double r = Kokkos::sqrt(y * y * ex * ex / (m_e * m_e * xi2) + ex2 * ex2 * 0.25); - double theta - = Kokkos::atan2(2. * y * ex, (m_e * xi * (m_epsilon * x * x - 2. * x - m_epsilon))); - if (theta < 0) { - theta = 2 * M_PI + theta; - } - return ddc::Coordinate(r, theta); - } - /** * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping. * @@ -215,13 +186,12 @@ class CzarnyToCartesian / (2 - Kokkos::sqrt(1 + m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta)))); } - /** * @brief Compute full Jacobian matrix. * * For some computations, we need the complete Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * jacobian_11, jacobian_12, jacobian_21 and jacobian_22. * * @param[in] coord @@ -355,7 +325,7 @@ class CzarnyToCartesian * * For some computations, we need the complete inverse Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22. * * @param[in] coord @@ -393,8 +363,6 @@ class CzarnyToCartesian matrix[1][1] = 1 / r * cos_theta / fact_3; } - - /** * @brief Compute the (1,1) coefficient of the inverse Jacobian matrix. * @@ -500,7 +468,6 @@ class CzarnyToCartesian return 1 / r * cos_theta / fact_3; } - /** * @brief Compute the full Jacobian matrix from the mapping to the pseudo-Cartesian mapping at the central point. * @@ -590,7 +557,7 @@ class CzarnyToCartesian namespace mapping_detail { template -struct MappingAccessibility> : std::true_type +struct MappingAccessibility> : std::true_type { }; } // namespace mapping_detail diff --git a/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp b/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp index 6517f865d..6df177635 100644 --- a/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp +++ b/vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp @@ -150,7 +150,7 @@ class DiscreteToCartesian * * For some computations, we need the complete Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * jacobian_11, jacobian_12, jacobian_21 and jacobian_22. * * @param[in] coord diff --git a/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp b/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp index d6f877c59..395e634a2 100644 --- a/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp +++ b/vendor/sll/include/sll/mapping/inverse_jacobian_matrix.hpp @@ -31,7 +31,7 @@ class InverseJacobianMatrix * * For some computations, we need the complete inverse Jacobian matrix or just the * coefficients. - * The coefficients can be given indendently with the functions + * The coefficients can be given independently with the functions * inv_jacobian_11, inv_jacobian_12, inv_jacobian_21 and inv_jacobian_22. * * @param[in] coord The coordinate where we evaluate the Jacobian matrix. diff --git a/vendor/sll/tests/mapping_execution_space_access.cpp b/vendor/sll/tests/mapping_execution_space_access.cpp index 6f32a1c38..cab17fb7f 100644 --- a/vendor/sll/tests/mapping_execution_space_access.cpp +++ b/vendor/sll/tests/mapping_execution_space_access.cpp @@ -3,6 +3,8 @@ #include +#include "sll/mapping/cartesian_to_circular.hpp" +#include "sll/mapping/cartesian_to_czarny.hpp" #include "sll/mapping/circular_to_cartesian.hpp" #include "sll/mapping/curvilinear2d_to_cartesian.hpp" #include "sll/mapping/czarny_to_cartesian.hpp" @@ -162,7 +164,7 @@ class MappingMemoryAccess : public ::testing::Test template double check_logical_to_physical_coord_converter( - Mapping const& mapping, + Mapping const& to_physical_mapping, CoordRTheta const& coord_rtheta, CoordXY const& coord_xy) { @@ -174,7 +176,7 @@ double check_logical_to_physical_coord_converter( Kokkos::RangePolicy(exec, 0, 1), KOKKOS_LAMBDA(int const i, double& err) { // Coord converter: logical -> physical. - CoordXY diff_coord_xy = mapping(coord_rtheta) - coord_xy; + CoordXY diff_coord_xy = to_physical_mapping(coord_rtheta) - coord_xy; err = Kokkos::max(ddc::get(diff_coord_xy), err); err = Kokkos::max(ddc::get(diff_coord_xy), err); }, @@ -184,7 +186,7 @@ double check_logical_to_physical_coord_converter( template double check_physical_to_logical_coord_converter( - Mapping const& mapping, + Mapping const& to_physical_mapping, CoordRTheta const& coord_rtheta, CoordXY const& coord_xy) { @@ -196,7 +198,7 @@ double check_physical_to_logical_coord_converter( Kokkos::RangePolicy(exec, 0, 1), KOKKOS_LAMBDA(int const i, double& err) { // Coord converter: physical -> logical. - CoordRTheta diff_coord_rtheta = mapping(coord_xy) - coord_rtheta; + CoordRTheta diff_coord_rtheta = to_physical_mapping(coord_xy) - coord_rtheta; err = Kokkos::max(ddc::get(diff_coord_rtheta), err); err = Kokkos::max(ddc::get(diff_coord_rtheta), err); }, @@ -207,14 +209,19 @@ double check_physical_to_logical_coord_converter( } // namespace - TEST_F(MappingMemoryAccess, HostCircularCoordConverter) { // Mapping - CircularToCartesian const mapping; + CircularToCartesian const to_physical_mapping; + CartesianToCircular const to_logical_mapping; + static_assert(is_2d_mapping_v>); + static_assert(is_2d_mapping_v>); + static_assert(is_accessible_v< + Kokkos::DefaultHostExecutionSpace, + CartesianToCircular>); static_assert(is_accessible_v< Kokkos::DefaultHostExecutionSpace, - CircularToCartesian>); + CircularToCartesian>); // Test coordinates double const r = .75; @@ -226,12 +233,12 @@ TEST_F(MappingMemoryAccess, HostCircularCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - CoordXY diff_coord_xy = mapping(coord_rtheta) - coord_xy; + CoordXY diff_coord_xy = to_physical_mapping(coord_rtheta) - coord_xy; EXPECT_LE(ddc::get(diff_coord_xy), 1e-15); EXPECT_LE(ddc::get(diff_coord_xy), 1e-15); // Coord converter: physical--> logical. - CoordRTheta diff_coord_rtheta = mapping(coord_xy) - coord_rtheta; + CoordRTheta diff_coord_rtheta = to_logical_mapping(coord_xy) - coord_rtheta; EXPECT_LE(ddc::get(diff_coord_rtheta), 1e-15); EXPECT_LE(ddc::get(diff_coord_rtheta), 1e-15); } @@ -242,9 +249,14 @@ TEST_F(MappingMemoryAccess, HostCzarnyCoordConverter) // Mapping double const epsilon = 0.3; double const e = 1.4; - CzarnyToCartesian const mapping(epsilon, e); + CzarnyToCartesian const to_physical_mapping(epsilon, e); + CartesianToCzarny const to_logical_mapping(epsilon, e); + static_assert(is_2d_mapping_v>); + static_assert(is_2d_mapping_v>); + static_assert( + is_accessible_v>); static_assert( - is_accessible_v>); + is_accessible_v>); // Test coordinates double const r = .75; @@ -258,12 +270,12 @@ TEST_F(MappingMemoryAccess, HostCzarnyCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - CoordXY diff_coord_xy = mapping(coord_rtheta) - coord_xy; + CoordXY diff_coord_xy = to_physical_mapping(coord_rtheta) - coord_xy; EXPECT_LE(ddc::get(diff_coord_xy), 1e-15); EXPECT_LE(ddc::get(diff_coord_xy), 1e-15); // Coord converter: physical -> logical. - CoordRTheta diff_coord_rtheta = mapping(coord_xy) - coord_rtheta; + CoordRTheta diff_coord_rtheta = to_logical_mapping(coord_xy) - coord_rtheta; EXPECT_LE(ddc::get(diff_coord_rtheta), 1e-15); EXPECT_LE(ddc::get(diff_coord_rtheta), 1e-15); } @@ -274,9 +286,9 @@ TEST_F(MappingMemoryAccess, HostDiscreteCoordConverter) // Mapping double const epsilon = 0.3; double const e = 1.4; - CzarnyToCartesian const analytical_mapping(epsilon, e); + CzarnyToCartesian const analytical_mapping(epsilon, e); static_assert( - is_accessible_v>); + is_accessible_v>); SplineRThetaBuilder builder(interpolation_idx_range_rtheta); @@ -294,8 +306,9 @@ TEST_F(MappingMemoryAccess, HostDiscreteCoordConverter) SplineRThetaBuilder, SplineRThetaEvaluator> mapping_builder(HostExecSpace(), analytical_mapping, builder, evaluator); - DiscreteToCartesian mapping = mapping_builder(); - static_assert(is_accessible_v); + DiscreteToCartesian to_physical_mapping = mapping_builder(); + static_assert( + is_accessible_v); // Test coordinates double const r = .75; @@ -309,19 +322,19 @@ TEST_F(MappingMemoryAccess, HostDiscreteCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - CoordXY diff_coord_xy = mapping(coord_rtheta) - coord_xy; + CoordXY diff_coord_xy = to_physical_mapping(coord_rtheta) - coord_xy; EXPECT_LE(ddc::get(diff_coord_xy), 1e-7); EXPECT_LE(ddc::get(diff_coord_xy), 1e-7); } - TEST_F(MappingMemoryAccess, DeviceCircularCoordConverter) { // Mapping - CircularToCartesian const mapping; + CircularToCartesian const to_physical_mapping; + CartesianToCircular const to_logical_mapping; static_assert( - is_accessible_v>); + is_accessible_v>); // Test coordinates double const r = .75; @@ -333,11 +346,14 @@ TEST_F(MappingMemoryAccess, DeviceCircularCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - double err = check_logical_to_physical_coord_converter(mapping, coord_rtheta, coord_xy); + double err = check_logical_to_physical_coord_converter( + to_physical_mapping, + coord_rtheta, + coord_xy); EXPECT_LE(err, 1e-15); // Coord converter: physical--> logical. - err = check_physical_to_logical_coord_converter(mapping, coord_rtheta, coord_xy); + err = check_physical_to_logical_coord_converter(to_logical_mapping, coord_rtheta, coord_xy); EXPECT_LE(err, 1e-15); } @@ -347,9 +363,11 @@ TEST_F(MappingMemoryAccess, DeviceCzarnyCoordConverter) // Mapping double const epsilon = 0.3; double const e = 1.4; - CzarnyToCartesian const mapping(epsilon, e); + CzarnyToCartesian const to_physical_mapping(epsilon, e); + CartesianToCzarny const to_logical_mapping(epsilon, e); + static_assert( - is_accessible_v>); + is_accessible_v>); // Test coordinates double const r = .75; @@ -363,11 +381,14 @@ TEST_F(MappingMemoryAccess, DeviceCzarnyCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - double err = check_logical_to_physical_coord_converter(mapping, coord_rtheta, coord_xy); + double err = check_logical_to_physical_coord_converter( + to_physical_mapping, + coord_rtheta, + coord_xy); EXPECT_LE(err, 1e-15); // Coord converter: physical--> logical. - err = check_physical_to_logical_coord_converter(mapping, coord_rtheta, coord_xy); + err = check_physical_to_logical_coord_converter(to_logical_mapping, coord_rtheta, coord_xy); EXPECT_LE(err, 1e-15); } @@ -377,9 +398,9 @@ TEST_F(MappingMemoryAccess, DeviceDiscreteCoordConverter) // Mapping double const epsilon = 0.3; double const e = 1.4; - CzarnyToCartesian const analytical_mapping(epsilon, e); + CzarnyToCartesian const analytical_mapping(epsilon, e); static_assert( - is_accessible_v>); + is_accessible_v>); SplineRThetaBuilder builder(interpolation_idx_range_rtheta); @@ -397,8 +418,8 @@ TEST_F(MappingMemoryAccess, DeviceDiscreteCoordConverter) SplineRThetaBuilder, SplineRThetaEvaluator> mapping_builder(DeviceExecSpace(), analytical_mapping, builder, evaluator); - DiscreteToCartesian mapping = mapping_builder(); - static_assert(is_accessible_v); + DiscreteToCartesian to_physical_mapping = mapping_builder(); + static_assert(is_accessible_v); // Test coordinates double const r = .75; @@ -412,6 +433,9 @@ TEST_F(MappingMemoryAccess, DeviceDiscreteCoordConverter) CoordXY const coord_xy(x, y); // Coord converter: logical -> physical. - double err = check_logical_to_physical_coord_converter(mapping, coord_rtheta, coord_xy); + double err = check_logical_to_physical_coord_converter( + to_physical_mapping, + coord_rtheta, + coord_xy); EXPECT_LE(err, 1e-7); } diff --git a/vendor/sll/tests/mapping_jacobian.cpp b/vendor/sll/tests/mapping_jacobian.cpp index 1d45f703e..46d915425 100644 --- a/vendor/sll/tests/mapping_jacobian.cpp +++ b/vendor/sll/tests/mapping_jacobian.cpp @@ -152,7 +152,7 @@ class InvJacobianMatrix : public testing::TestWithParam mapping; + const CircularToCartesian mapping; CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -179,7 +179,7 @@ TEST_P(InvJacobianMatrix, InverseMatrixCircMap) theta_min + dtheta * ddc::select(irp).uid()); }); - static_assert(has_2d_jacobian_v, CoordRTheta>); + static_assert(has_2d_jacobian_v, CoordRTheta>); InverseJacobianMatrix inv_jacobian(mapping); // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix @@ -197,7 +197,7 @@ TEST_P(InvJacobianMatrix, InverseMatrixCircMap) TEST_P(InvJacobianMatrix, InverseMatrixCzarMap) { auto const [Nr, Nt] = GetParam(); - const CzarnyToCartesian mapping(0.3, 1.4); + const CzarnyToCartesian mapping(0.3, 1.4); CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -224,8 +224,8 @@ TEST_P(InvJacobianMatrix, InverseMatrixCzarMap) theta_min + dtheta * ddc::select(irp).uid()); }); - static_assert(has_2d_jacobian_v, CoordRTheta>); - static_assert(has_2d_inv_jacobian_v, CoordRTheta>); + static_assert(has_2d_jacobian_v, CoordRTheta>); + static_assert(has_2d_inv_jacobian_v, CoordRTheta>); // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix ddc::for_each(grid, [&](IdxRTheta const irp) { @@ -243,7 +243,7 @@ TEST_P(InvJacobianMatrix, InverseMatrixCzarMap) TEST_P(InvJacobianMatrix, InverseMatrixDiscCzarMap) { auto const [Nr, Nt] = GetParam(); - const CzarnyToCartesian analytical_mapping(0.3, 1.4); + const CzarnyToCartesian analytical_mapping(0.3, 1.4); CoordR const r_min(0.0); CoordR const r_max(1.0); diff --git a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp index e9ad28378..c4d9d634c 100644 --- a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp +++ b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp @@ -150,7 +150,7 @@ class JacobianMatrixAndJacobianCoefficients TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCircMap) { auto const [Nr, Nt] = GetParam(); - const CircularToCartesian mapping; + const CircularToCartesian mapping; CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -217,7 +217,7 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCircMap) TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCzarMap) { auto const [Nr, Nt] = GetParam(); - const CzarnyToCartesian mapping(0.3, 1.4); + const CzarnyToCartesian mapping(0.3, 1.4); CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -283,7 +283,7 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCzarMap) TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixDiscCzarMap) { auto const [Nr, Nt] = GetParam(); - const CzarnyToCartesian analytical_mapping(0.3, 1.4); + const CzarnyToCartesian analytical_mapping(0.3, 1.4); CoordR const r_min(0.0); CoordR const r_max(1.0); diff --git a/vendor/sll/tests/metric_tensor.cpp b/vendor/sll/tests/metric_tensor.cpp index 0930dbbdd..b42032439 100644 --- a/vendor/sll/tests/metric_tensor.cpp +++ b/vendor/sll/tests/metric_tensor.cpp @@ -108,7 +108,7 @@ class InverseMetricTensor : public testing::TestWithParam mapping; + const CircularToCartesian mapping; CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -135,7 +135,7 @@ TEST_P(InverseMetricTensor, InverseMatrixCircMap) theta_min + dp * ddc::select(irp).uid()); }); - MetricTensor, CoordRTheta> metric_tensor(mapping); + MetricTensor, CoordRTheta> metric_tensor(mapping); // Test for each coordinates if the inverse_metric_tensor is the inverse of the metric_tensor ddc::for_each(grid, [&](IdxRTheta const irp) { Matrix_2x2 matrix; @@ -153,7 +153,7 @@ TEST_P(InverseMetricTensor, InverseMatrixCircMap) TEST_P(InverseMetricTensor, InverseMatrixCzarMap) { auto const [Nr, Nt] = GetParam(); - const CzarnyToCartesian mapping(0.3, 1.4); + const CzarnyToCartesian mapping(0.3, 1.4); CoordR const r_min(0.0); CoordR const r_max(1.0); @@ -180,7 +180,7 @@ TEST_P(InverseMetricTensor, InverseMatrixCzarMap) theta_min + dp * ddc::select(irp).uid()); }); - MetricTensor, CoordRTheta> metric_tensor(mapping); + MetricTensor, CoordRTheta> metric_tensor(mapping); // Test for each coordinates if the inverse_metric_tensor is the inverse of the metric_tensor ddc::for_each(grid, [&](IdxRTheta const irp) { Matrix_2x2 matrix; diff --git a/vendor/sll/tests/polar_bsplines.cpp b/vendor/sll/tests/polar_bsplines.cpp index f6c3f7ab5..cc705c8e6 100644 --- a/vendor/sll/tests/polar_bsplines.cpp +++ b/vendor/sll/tests/polar_bsplines.cpp @@ -93,7 +93,7 @@ TYPED_TEST(PolarBsplineFixture, PartitionOfUnity) using PolarCoord = ddc::Coordinate; using BSplinesR = typename TestFixture::BSplinesR; using BSplinesTheta = typename TestFixture::BSplinesTheta; - using CircToCart = CircularToCartesian; + using CircToCart = CircularToCartesian; using SplineRThetaBuilder = ddc::SplineBuilder2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, diff --git a/vendor/sll/tests/polar_splines.cpp b/vendor/sll/tests/polar_splines.cpp index bcbe5d2b2..0bb57b0da 100644 --- a/vendor/sll/tests/polar_splines.cpp +++ b/vendor/sll/tests/polar_splines.cpp @@ -67,9 +67,9 @@ struct BSplines : PolarBSplines }; #if defined(CIRCULAR_MAPPING) -using CircToCart = CircularToCartesian; +using CircToCart = CircularToCartesian; #elif defined(CZARNY_MAPPING) -using CircToCart = CzarnyToCartesian; +using CircToCart = CzarnyToCartesian; #endif TEST(PolarSplineTest, ConstantEval) diff --git a/vendor/sll/tests/pseudo_cartesian_matrix.cpp b/vendor/sll/tests/pseudo_cartesian_matrix.cpp index c6184ff71..6a2270b3c 100644 --- a/vendor/sll/tests/pseudo_cartesian_matrix.cpp +++ b/vendor/sll/tests/pseudo_cartesian_matrix.cpp @@ -201,7 +201,7 @@ class PseudoCartesianJacobianMatrixTest // --- CIRCULAR MAPPING --------------------------------------------------------------------------- std::cout << " - Nr x Nt = " << Nr << " x " << Nt << std::endl << " - Circular mapping: "; - const CircularToCartesian analytical_mapping_circ; + const CircularToCartesian analytical_mapping_circ; DiscreteToCartesianBuilder mapping_builder_circ( Kokkos::DefaultHostExecutionSpace(), @@ -220,7 +220,7 @@ class PseudoCartesianJacobianMatrixTest // --- CZARNY MAPPING ----------------------------------------------------------------------------- std::cout << " - Czarny mapping: "; - const CzarnyToCartesian analytical_mapping_czar(0.3, 1.4); + const CzarnyToCartesian analytical_mapping_czar(0.3, 1.4); DiscreteToCartesianBuilder mapping_builder_czar( Kokkos::DefaultHostExecutionSpace(), diff --git a/vendor/sll/tests/refined_discrete_mapping.cpp b/vendor/sll/tests/refined_discrete_mapping.cpp index 536552c2a..425bf82c6 100644 --- a/vendor/sll/tests/refined_discrete_mapping.cpp +++ b/vendor/sll/tests/refined_discrete_mapping.cpp @@ -121,8 +121,8 @@ using FieldMemRTheta = ddc::Chunk; using IdxRangeRTheta = ddc::DiscreteDomain; -using CzarnyMapping = CzarnyToCartesian; -using CircularMapping = CircularToCartesian; +using CzarnyMapping = CzarnyToCartesian; +using CircularMapping = CircularToCartesian; template