diff --git a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp index 1fbf02d46..3755cbebd 100644 --- a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp +++ b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp @@ -53,30 +53,30 @@ int main(int argc, char** argv) // Reading config // --> Mesh info - IDomainX const mesh_x = init_spline_dependent_idx_range< - IDimX, + IdxRangeX const mesh_x = init_spline_dependent_idx_range< + GridX, BSplinesX, SplineInterpPointsX>(conf_voicexx, "x"); - IDomainY const mesh_y = init_spline_dependent_idx_range< - IDimY, + IdxRangeY const mesh_y = init_spline_dependent_idx_range< + GridY, BSplinesY, SplineInterpPointsY>(conf_voicexx, "y"); - IDomainVx const mesh_vx = init_spline_dependent_idx_range< - IDimVx, + IdxRangeVx const mesh_vx = init_spline_dependent_idx_range< + GridVx, BSplinesVx, SplineInterpPointsVx>(conf_voicexx, "vx"); - IDomainVy const mesh_vy = init_spline_dependent_idx_range< - IDimVy, + IdxRangeVy const mesh_vy = init_spline_dependent_idx_range< + GridVy, BSplinesVy, SplineInterpPointsVy>(conf_voicexx, "vy"); - IDomainXY const mesh_xy(mesh_x, mesh_y); - IDomainVxVy mesh_vxvy(mesh_vx, mesh_vy); - IDomainXYVxVy const meshXYVxVy(mesh_x, mesh_y, mesh_vx, mesh_vy); + IdxRangeXY const mesh_xy(mesh_x, mesh_y); + IdxRangeVxVy mesh_vxvy(mesh_vx, mesh_vy); + IdxRangeXYVxVy const meshXYVxVy(mesh_x, mesh_y, mesh_vx, mesh_vy); IdxRangeSp const dom_kinsp = init_species(conf_voicexx); - IDomainSpVxVy const meshSpVxVy(dom_kinsp, mesh_vx, mesh_vy); - IDomainSpXYVxVy const meshSpXYVxVy(dom_kinsp, meshXYVxVy); + IdxRangeSpVxVy const meshSpVxVy(dom_kinsp, mesh_vx, mesh_vy); + IdxRangeSpXYVxVy const meshSpXYVxVy(dom_kinsp, meshXYVxVy); SplineXBuilder const builder_x(meshXYVxVy); SplineYBuilder const builder_y(meshXYVxVy); @@ -86,15 +86,15 @@ int main(int argc, char** argv) SplineVyBuilder_1d const builder_vy_1d(mesh_vy); // Initialization of the distribution function - DFieldSpVxVy allfequilibrium(meshSpVxVy); + DFieldMemSpVxVy allfequilibrium(meshSpVxVy); MaxwellianEquilibrium const init_fequilibrium = MaxwellianEquilibrium::init_from_input(dom_kinsp, conf_voicexx); init_fequilibrium(allfequilibrium); - DFieldSpXYVxVy allfdistribu(meshSpXYVxVy); + DFieldMemSpXYVxVy allfdistribu(meshSpXYVxVy); SingleModePerturbInitialization const init = SingleModePerturbInitialization:: init_from_input(allfequilibrium, dom_kinsp, conf_voicexx); init(allfdistribu); - auto allfequilibrium_host = ddc::create_mirror_view_and_copy(allfequilibrium.span_view()); + auto allfequilibrium_host = ddc::create_mirror_view_and_copy(get_field(allfequilibrium)); // --> Algorithm info double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat"); @@ -105,44 +105,44 @@ int main(int argc, char** argv) int const nbstep_diag = int(time_diag / deltat); // Create spline evaluator - ddc::PeriodicExtrapolationRule bv_x_min; - ddc::PeriodicExtrapolationRule bv_x_max; + ddc::PeriodicExtrapolationRule bv_x_min; + ddc::PeriodicExtrapolationRule bv_x_max; SplineXEvaluator const spline_x_evaluator(bv_x_min, bv_x_max); PreallocatableSplineInterpolator const spline_x_interpolator(builder_x, spline_x_evaluator); - ddc::PeriodicExtrapolationRule bv_y_min; - ddc::PeriodicExtrapolationRule bv_y_max; + ddc::PeriodicExtrapolationRule bv_y_min; + ddc::PeriodicExtrapolationRule bv_y_max; SplineYEvaluator const spline_y_evaluator(bv_y_min, bv_y_max); PreallocatableSplineInterpolator const spline_y_interpolator(builder_y, spline_y_evaluator); - ddc::ConstantExtrapolationRule bv_vx_min(ddc::coordinate(mesh_vx.front())); - ddc::ConstantExtrapolationRule bv_vx_max(ddc::coordinate(mesh_vx.back())); + ddc::ConstantExtrapolationRule bv_vx_min(ddc::coordinate(mesh_vx.front())); + ddc::ConstantExtrapolationRule bv_vx_max(ddc::coordinate(mesh_vx.back())); SplineVxEvaluator const spline_vx_evaluator(bv_vx_min, bv_vx_max); PreallocatableSplineInterpolator const spline_vx_interpolator(builder_vx, spline_vx_evaluator); - ddc::ConstantExtrapolationRule bv_vy_min(ddc::coordinate(mesh_vy.front())); - ddc::ConstantExtrapolationRule bv_vy_max(ddc::coordinate(mesh_vy.back())); + ddc::ConstantExtrapolationRule bv_vy_min(ddc::coordinate(mesh_vy.front())); + ddc::ConstantExtrapolationRule bv_vy_max(ddc::coordinate(mesh_vy.back())); SplineVyEvaluator const spline_vy_evaluator(bv_vy_min, bv_vy_max); PreallocatableSplineInterpolator const spline_vy_interpolator(builder_vy, spline_vy_evaluator); // Create advection operator - BslAdvectionSpatial const advection_x(spline_x_interpolator); - BslAdvectionSpatial const advection_y(spline_y_interpolator); - BslAdvectionVelocity const advection_vx(spline_vx_interpolator); - BslAdvectionVelocity const advection_vy(spline_vy_interpolator); + BslAdvectionSpatial const advection_x(spline_x_interpolator); + BslAdvectionSpatial const advection_y(spline_y_interpolator); + BslAdvectionVelocity const advection_vx(spline_vx_interpolator); + BslAdvectionVelocity const advection_vy(spline_vy_interpolator); SplitVlasovSolver const vlasov(advection_x, advection_y, advection_vx, advection_vy); - host_t const quadrature_coeffs_host + host_t const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(mesh_vxvy, builder_vx_1d, builder_vy_1d); auto quadrature_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); - FFTPoissonSolver fft_poisson_solver( + get_field(quadrature_coeffs_host)); + FFTPoissonSolver fft_poisson_solver( mesh_xy); ChargeDensityCalculator const rhs(quadrature_coeffs); QNSolver const poisson(fft_poisson_solver, rhs); @@ -167,7 +167,7 @@ int main(int argc, char** argv) steady_clock::time_point const start = steady_clock::now(); - predcorr(allfdistribu.span_view(), deltat, nbiter); + predcorr(get_field(allfdistribu), deltat, nbiter); steady_clock::time_point const end = steady_clock::now(); diff --git a/src/geometryRTheta/poisson/polarpoissonlikesolver.hpp b/src/geometryRTheta/poisson/polarpoissonlikesolver.hpp index 5a8b6cfcd..f043ef172 100644 --- a/src/geometryRTheta/poisson/polarpoissonlikesolver.hpp +++ b/src/geometryRTheta/poisson/polarpoissonlikesolver.hpp @@ -159,7 +159,7 @@ class PolarSplineFEMPoissonLikeSolver const int nbasis_p; // Domains - BSDomainPolar fem_non_singular_domain; + BSDomainPolar fem_non_singular_idx_range; BSDomainR_Polar radial_bsplines; BSDomainP_Polar polar_bsplines; @@ -216,8 +216,8 @@ class PolarSplineFEMPoissonLikeSolver Mapping const& mapping) : nbasis_r(ddc::discrete_space().nbasis() - n_overlap_cells - 1) , nbasis_p(ddc::discrete_space().nbasis()) - , fem_non_singular_domain( - ddc::discrete_space().tensor_bspline_domain().remove_last( + , fem_non_singular_idx_range( + ddc::discrete_space().tensor_bspline_idx_range().remove_last( ddc::DiscreteVector {nbasis_p})) , radial_bsplines(ddc::discrete_space().full_domain().remove_first( ddc::DiscreteVector {n_overlap_cells})) @@ -240,7 +240,7 @@ class PolarSplineFEMPoissonLikeSolver , weights_r(quadrature_domain_r) , weights_p(quadrature_domain_p) , singular_basis_vals_and_derivs(ddc::DiscreteDomain( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), ddc::select(quadrature_domain_singular), ddc::select(quadrature_domain_singular))) , r_basis_vals_and_derivs(ddc::DiscreteDomain< @@ -316,7 +316,7 @@ class PolarSplineFEMPoissonLikeSolver } }); - auto singular_domain = PolarBSplinesRP::singular_domain(); + auto singular_idx_range = PolarBSplinesRP::singular_idx_range(); // Find value and derivative of 2D bsplines covering the singular point ddc::for_each(quadrature_domain_singular, [&](QuadratureIndexRP const irp) { @@ -335,20 +335,20 @@ class PolarSplineFEMPoissonLikeSolver // Calculate the value ddc::discrete_space().eval_basis(singular_vals, vals, coord); - for (auto ib : singular_domain) { + for (auto ib : singular_idx_range) { singular_basis_vals_and_derivs(ib, ir, ip).value = singular_vals[ib.uid()]; } // Calculate the radial derivative ddc::discrete_space().eval_deriv_r(singular_vals, vals, coord); - for (auto ib : singular_domain) { + for (auto ib : singular_idx_range) { singular_basis_vals_and_derivs(ib, ir, ip).radial_derivative = singular_vals[ib.uid()]; } // Calculate the poloidal derivative ddc::discrete_space().eval_deriv_p(singular_vals, vals, coord); - for (auto ib : singular_domain) { + for (auto ib : singular_idx_range) { singular_basis_vals_and_derivs(ib, ir, ip).poloidal_derivative = singular_vals[ib.uid()]; } @@ -399,8 +399,8 @@ class PolarSplineFEMPoissonLikeSolver int matrix_idx(0); // Calculate the matrix elements corresponding to the bsplines which cover the singular point - ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_test) { - ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_trial) { + ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_test) { + ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_trial) { // Calculate the weak integral row_coo_host(matrix_idx) = idx_test.uid(); col_coo_host(matrix_idx) = idx_trial.uid(); @@ -431,12 +431,13 @@ class PolarSplineFEMPoissonLikeSolver BSDomainR central_radial_bspline_domain(radial_bsplines.take_first( ddc::DiscreteVector {BSplinesR_Polar::degree()})); - BSDomainRP non_singular_domain_near_centre(central_radial_bspline_domain, polar_bsplines); + BSDomainRP + non_singular_idx_range_near_centre(central_radial_bspline_domain, polar_bsplines); // Calculate the matrix elements where bspline products overlap the bsplines which cover the singular point - ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_test) { + ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_test) { ddc::for_each( - non_singular_domain_near_centre, + non_singular_idx_range_near_centre, [&](IDimBSpline2D_Polar const idx_trial) { const IndexPolarBspl polar_idx_trial( PolarBSplinesRP::get_polar_index(idx_trial)); @@ -520,7 +521,7 @@ class PolarSplineFEMPoissonLikeSolver assert(matrix_idx == n_elements_singular + n_elements_overlap); // Calculate the matrix elements following a stencil - ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const polar_idx_test) { + ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const polar_idx_test) { const IDimBSpline2D_Polar idx_test(PolarBSplinesRP::get_2d_index(polar_idx_test)); const std::size_t r_idx_test( ddc::select(idx_test).uid()); @@ -640,7 +641,7 @@ class PolarSplineFEMPoissonLikeSolver // Fill b ddc::for_each( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), [&](IndexPolarBspl const idx) { b_host(0, idx.uid()) = ddc::transform_reduce( quadrature_domain_singular, @@ -656,7 +657,7 @@ class PolarSplineFEMPoissonLikeSolver }); }); const std::size_t ncells_r = ddc::discrete_space().ncells(); - ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const idx) { + ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const idx) { const IDimBSpline2D_Polar idx_2d(PolarBSplinesRP::get_2d_index(idx)); const std::size_t r_idx(ddc::select(idx_2d).uid()); const std::size_t p_idx(ddc::select(idx_2d).uid()); @@ -725,11 +726,11 @@ class PolarSplineFEMPoissonLikeSolver // Fill the spline ddc::for_each( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), [&](IndexPolarBspl const idx) { spline.singular_spline_coef(idx) = b_host(0, idx.uid()); }); - ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const idx) { + ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const idx) { const IDimBSpline2D_Polar idx_2d(PolarBSplinesRP::get_2d_index(idx)); spline.spline_coef(idx_2d) = b_host(0, idx.uid()); }); @@ -774,7 +775,7 @@ class PolarSplineFEMPoissonLikeSolver { BSDomainP_Polar polar_domain(ddc::discrete_space().full_domain()); SplinePolar - spline(PolarBSplinesRP::singular_domain(), + spline(PolarBSplinesRP::singular_idx_range(), BSDomainRP(radial_bsplines, polar_domain)); (*this)(rhs, spline); diff --git a/src/geometryRTheta/time_solver/bsl_predcorr.hpp b/src/geometryRTheta/time_solver/bsl_predcorr.hpp index f922facc3..816673d70 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr.hpp @@ -122,7 +122,7 @@ class BslPredCorrRP : public ITimeSolverRP BSDomainP polar_domain(ddc::discrete_space().full_domain()); SplinePolar electrostatic_potential_coef( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), BSDomainRP(radial_bsplines, polar_domain)); ddc::NullExtrapolationRule extrapolation_rule; PolarSplineEvaluator polar_spline_evaluator( diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp index da712866d..90c890b0f 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp @@ -155,7 +155,7 @@ class BslExplicitPredCorrRP : public ITimeSolverRP DFieldRP electrical_potential(grid); SplinePolar electrostatic_potential_coef( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), BSDomainRP(radial_bsplines, polar_domain)); ddc::NullExtrapolationRule extrapolation_rule; diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp index 568d545e5..0a9822cd8 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp @@ -155,7 +155,7 @@ class BslImplicitPredCorrRP : public ITimeSolverRP DFieldRP electrical_potential(grid); SplinePolar electrostatic_potential_coef( - PolarBSplinesRP::singular_domain(), + PolarBSplinesRP::singular_idx_range(), BSDomainRP(radial_bsplines, polar_domain)); ddc::NullExtrapolationRule extrapolation_rule; diff --git a/src/geometryXVx/boltzmann/splitvlasovsolver.hpp b/src/geometryXVx/boltzmann/splitvlasovsolver.hpp index 64b114df9..982a5d72d 100644 --- a/src/geometryXVx/boltzmann/splitvlasovsolver.hpp +++ b/src/geometryXVx/boltzmann/splitvlasovsolver.hpp @@ -19,7 +19,7 @@ class IAdvectionVelocity; * The Vlasov equation is split between two advection equations * along the X and Vx directions. The splitting involves solving * the x-direction advection first on a time interval of length dt/2, - * then the vx-direction advection on a tim dt, and then x-direction + * then the vx-direction advection on a time dt, and then x-direction * again on dt/2. */ class SplitVlasovSolver : public IBoltzmannSolver diff --git a/src/geometryXYVxVy/geometry/geometry.hpp b/src/geometryXYVxVy/geometry/geometry.hpp index 788814ac0..65e0ddabe 100644 --- a/src/geometryXYVxVy/geometry/geometry.hpp +++ b/src/geometryXYVxVy/geometry/geometry.hpp @@ -1,7 +1,6 @@ // SPDX-License-Identifier: MIT #pragma once - #include #include #include @@ -9,10 +8,12 @@ #include #include +#include "ddc_aliases.hpp" + /** * @brief A class which describes the real space in the first spatial direction X. */ -struct RDimX +struct X { /** * @brief A boolean indicating if the dimension is periodic. @@ -23,7 +24,7 @@ struct RDimX /** * @brief A class which describes the real space in the second spatial direction Y. */ -struct RDimY +struct Y { /** * @brief A boolean indicating if the dimension is periodic. @@ -34,7 +35,7 @@ struct RDimY /** * @brief A class which describes the real space in the second velocity direction X. */ -struct RDimVx +struct Vx { /** * @brief A boolean indicating if the dimension is periodic. @@ -45,7 +46,7 @@ struct RDimVx /** * @brief A class which describes the real space in the second velocity direction Y. */ -struct RDimVy +struct Vy { /** * @brief A boolean indicating if the dimension is periodic. @@ -53,12 +54,12 @@ struct RDimVy static bool constexpr PERIODIC = false; }; -using CoordX = ddc::Coordinate; -using CoordY = ddc::Coordinate; -using CoordXY = ddc::Coordinate; +using CoordX = Coord; +using CoordY = Coord; +using CoordXY = Coord; -using CoordVx = ddc::Coordinate; -using CoordVy = ddc::Coordinate; +using CoordVx = Coord; +using CoordVy = Coord; int constexpr BSDegreeX = 3; int constexpr BSDegreeY = 3; @@ -75,30 +76,30 @@ bool constexpr BsplineOnUniformCellsVy = true; struct BSplinesX : std::conditional_t< BsplineOnUniformCellsX, - ddc::UniformBSplines, - ddc::NonUniformBSplines> + ddc::UniformBSplines, + ddc::NonUniformBSplines> { }; struct BSplinesY : std::conditional_t< BsplineOnUniformCellsY, - ddc::UniformBSplines, - ddc::NonUniformBSplines> + ddc::UniformBSplines, + ddc::NonUniformBSplines> { }; struct BSplinesVx : std::conditional_t< BsplineOnUniformCellsVx, - ddc::UniformBSplines, - ddc::NonUniformBSplines> + ddc::UniformBSplines, + ddc::NonUniformBSplines> { }; struct BSplinesVy : std::conditional_t< BsplineOnUniformCellsVy, - ddc::UniformBSplines, - ddc::NonUniformBSplines> + ddc::UniformBSplines, + ddc::NonUniformBSplines> { }; @@ -118,16 +119,16 @@ using SplineInterpPointsVy = ddc::GrevilleInterpolationPoints; // IDim definition -struct IDimX : SplineInterpPointsX::interpolation_discrete_dimension_type +struct GridX : SplineInterpPointsX::interpolation_discrete_dimension_type { }; -struct IDimY : SplineInterpPointsY::interpolation_discrete_dimension_type +struct GridY : SplineInterpPointsY::interpolation_discrete_dimension_type { }; -struct IDimVx : SplineInterpPointsVx::interpolation_discrete_dimension_type +struct GridVx : SplineInterpPointsVx::interpolation_discrete_dimension_type { }; -struct IDimVy : SplineInterpPointsVy::interpolation_discrete_dimension_type +struct GridVy : SplineInterpPointsVy::interpolation_discrete_dimension_type { }; @@ -136,282 +137,282 @@ using SplineXBuilder = ddc::SplineBuilder< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesX, - IDimX, + GridX, SplineXBoundary, SplineXBoundary, ddc::SplineSolver::LAPACK, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridX, + GridY, + GridVx, + GridVy>; using SplineXEvaluator = ddc::SplineEvaluator< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesX, - IDimX, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridX, + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + GridX, + GridY, + GridVx, + GridVy>; using SplineYBuilder = ddc::SplineBuilder< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesY, - IDimY, + GridY, SplineYBoundary, SplineYBoundary, ddc::SplineSolver::LAPACK, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridX, + GridY, + GridVx, + GridVy>; using SplineYEvaluator = ddc::SplineEvaluator< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesY, - IDimY, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridY, + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + GridX, + GridY, + GridVx, + GridVy>; using SplineVxBuilder = ddc::SplineBuilder< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesVx, - IDimVx, + GridVx, SplineVxBoundary, SplineVxBoundary, ddc::SplineSolver::LAPACK, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridX, + GridY, + GridVx, + GridVy>; using SplineVxEvaluator = ddc::SplineEvaluator< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesVx, - IDimVx, - ddc::ConstantExtrapolationRule, - ddc::ConstantExtrapolationRule, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridVx, + ddc::ConstantExtrapolationRule, + ddc::ConstantExtrapolationRule, + GridX, + GridY, + GridVx, + GridVy>; using SplineVyBuilder = ddc::SplineBuilder< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesVy, - IDimVy, + GridVy, SplineVyBoundary, SplineVyBoundary, ddc::SplineSolver::LAPACK, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridX, + GridY, + GridVx, + GridVy>; using SplineVyEvaluator = ddc::SplineEvaluator< Kokkos::DefaultExecutionSpace, Kokkos::DefaultExecutionSpace::memory_space, BSplinesVy, - IDimVy, - ddc::ConstantExtrapolationRule, - ddc::ConstantExtrapolationRule, - IDimX, - IDimY, - IDimVx, - IDimVy>; + GridVy, + ddc::ConstantExtrapolationRule, + ddc::ConstantExtrapolationRule, + GridX, + GridY, + GridVx, + GridVy>; using SplineVxBuilder_1d = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesVx, - IDimVx, + GridVx, SplineVxBoundary, SplineVxBoundary, ddc::SplineSolver::LAPACK, - IDimVx>; + GridVx>; using SplineVyBuilder_1d = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesVy, - IDimVy, + GridVy, SplineVyBoundary, SplineVyBoundary, ddc::SplineSolver::LAPACK, - IDimVy>; + GridVy>; -using BSDomainX = ddc::DiscreteDomain; -using BSDomainY = ddc::DiscreteDomain; -using BSDomainXY = ddc::DiscreteDomain; -using BSDomainVx = ddc::DiscreteDomain; -using BSDomainVy = ddc::DiscreteDomain; -using BSDomainVxVy = ddc::DiscreteDomain; +using BSIdxRangeX = IdxRange; +using BSIdxRangeY = IdxRange; +using BSIdxRangeXY = IdxRange; +using BSIdxRangeVx = IdxRange; +using BSIdxRangeVy = IdxRange; +using BSIdxRangeVxVy = IdxRange; template -using BSViewXY = device_t>; -using DBSViewXY = BSViewXY; +using BSConstFieldXY = Field; +using DBSConstFieldXY = BSConstFieldXY; // Index -using IndexX = ddc::DiscreteElement; -using IndexY = ddc::DiscreteElement; -using IndexXY = ddc::DiscreteElement; -using IndexVx = ddc::DiscreteElement; -using IndexVy = ddc::DiscreteElement; -using IndexVxVy = ddc::DiscreteElement; -using IndexXYVxVy = ddc::DiscreteElement; -using IndexSpXYVxVy = ddc::DiscreteElement; +using IdxX = Idx; +using IdxY = Idx; +using IdxXY = Idx; +using IdxVx = Idx; +using IdxVy = Idx; +using IdxVxVy = Idx; +using IdxXYVxVy = Idx; +using IdxSpXYVxVy = Idx; // IVect definition -using IVectX = ddc::DiscreteVector; -using IVectY = ddc::DiscreteVector; -using IVectVx = ddc::DiscreteVector; -using IVectVy = ddc::DiscreteVector; - -// Idomain definition -using IDomainX = ddc::DiscreteDomain; -using IDomainY = ddc::DiscreteDomain; -using IDomainXY = ddc::DiscreteDomain; -using IDomainVx = ddc::DiscreteDomain; -using IDomainVy = ddc::DiscreteDomain; -using IDomainXYVxVy = ddc::DiscreteDomain; -using IDomainVxVy = ddc::DiscreteDomain; -using IDomainSpVxVy = ddc::DiscreteDomain; -using IDomainSpXYVxVy = ddc::DiscreteDomain; +using IdxStepX = IdxStep; +using IdxStepY = IdxStep; +using IdxStepVx = IdxStep; +using IdxStepVy = IdxStep; + +// Iindex range definition +using IdxRangeX = IdxRange; +using IdxRangeY = IdxRange; +using IdxRangeXY = IdxRange; +using IdxRangeVx = IdxRange; +using IdxRangeVy = IdxRange; +using IdxRangeXYVxVy = IdxRange; +using IdxRangeVxVy = IdxRange; +using IdxRangeSpVxVy = IdxRange; +using IdxRangeSpXYVxVy = IdxRange; template -using FieldX = device_t>; -using DFieldX = FieldX; +using FieldMemX = FieldMem; +using DFieldMemX = FieldMemX; template -using FieldY = device_t>; -using DFieldY = FieldY; +using FieldMemY = FieldMem; +using DFieldMemY = FieldMemY; template -using FieldXY = device_t>; -using DFieldXY = FieldXY; +using FieldMemXY = FieldMem; +using DFieldMemXY = FieldMemXY; template -using FieldVx = device_t>; +using FieldMemVx = FieldMem; template -using FieldVy = device_t>; +using FieldMemVy = FieldMem; template -using FieldVxVy = device_t>; -using DFieldVxVy = FieldVxVy; +using FieldMemVxVy = FieldMem; +using DFieldMemVxVy = FieldMemVxVy; template -using FieldXYVxVy = device_t>; -using DFieldXYVxVy = FieldXYVxVy; +using FieldMemXYVxVy = FieldMem; +using DFieldMemXYVxVy = FieldMemXYVxVy; template -using FieldSpVxVy = device_t>; -using DFieldSpVxVy = FieldSpVxVy; +using FieldMemSpVxVy = FieldMem; +using DFieldMemSpVxVy = FieldMemSpVxVy; template -using FieldSpXYVxVy = device_t>; -using DFieldSpXYVxVy = FieldSpXYVxVy; +using FieldMemSpXYVxVy = FieldMem; +using DFieldMemSpXYVxVy = FieldMemSpXYVxVy; // Span definition template -using SpanX = device_t>; -using DSpanX = SpanX; +using FieldX = Field; +using DFieldX = FieldX; template -using SpanY = device_t>; -using DSpanY = SpanY; +using FieldY = Field; +using DFieldY = FieldY; template -using SpanXY = device_t>; -using DSpanXY = SpanXY; +using FieldXY = Field; +using DFieldXY = FieldXY; template -using SpanVx = device_t>; -using DSpanVx = SpanVx; +using FieldVx = Field; +using DFieldVx = FieldVx; template -using SpanVy = device_t>; -using DSpanVy = SpanVy; +using FieldVy = Field; +using DFieldVy = FieldVy; template -using SpanVxVy = device_t>; -using DSpanVxVy = SpanVxVy; +using FieldVxVy = Field; +using DFieldVxVy = FieldVxVy; template -using SpanSpVxVy = device_t>; -using DSpanSpVxVy = SpanSpVxVy; +using FieldSpVxVy = Field; +using DFieldSpVxVy = FieldSpVxVy; template -using SpanSpXYVxVy = device_t>; -using DSpanSpXYVxVy = SpanSpXYVxVy; +using FieldSpXYVxVy = Field; +using DFieldSpXYVxVy = FieldSpXYVxVy; // View definition template -using ViewX = device_t>; +using ConstFieldX = Field; template -using ViewY = device_t>; +using ConstFieldY = Field; template -using ViewXY = device_t>; -using DViewXY = ViewXY; +using ConstFieldXY = Field; +using DConstFieldXY = ConstFieldXY; template -using ViewVx = device_t>; +using ConstFieldVx = Field; template -using ViewVy = device_t>; +using ConstFieldVy = Field; template -using ViewVxVy = device_t>; -using DViewVxVy = ViewVxVy; +using ConstFieldVxVy = Field; +using DConstFieldVxVy = ConstFieldVxVy; template -using ViewSpVxVy = device_t>; -using DViewSpVxVy = ViewSpVxVy; +using ConstFieldSpVxVy = Field; +using DConstFieldSpVxVy = ConstFieldSpVxVy; template -using ViewSpXYVxVy = device_t>; -using DViewSpXYVxVy = ViewSpXYVxVy; +using ConstFieldSpXYVxVy = Field; +using DConstFieldSpXYVxVy = ConstFieldSpXYVxVy; /** - * @brief A class providing aliases for useful subdomains of the geometry. It is used as template parameter for generic dimensionality-agnostic operat + * @brief A class providing aliases for useful subindex ranges of the geometry. It is used as template parameter for generic dimensionality-agnostic operat ors such as advections. */ class GeometryXYVxVy { public: /** - * @brief A templated type giving the velocity discrete dimension type associated to a spatial discrete dimension type. + * @brief A templated type giving the velocity discretised dimension type associated to a spatial discretised dimension type. */ template using velocity_dim_for = std::conditional_t< - std::is_same_v, - IDimVx, - std::conditional_t, IDimVy, void>>; + std::is_same_v, + GridVx, + std::conditional_t, GridVy, void>>; /** - * @brief A templated type giving the spatial discrete dimension type associated to a velocity discrete dimension type. + * @brief A templated type giving the spatial discretised dimension type associated to a velocity discretised dimension type. */ // template - // using spatial_dim_for = std::conditional_t, IDimX, std::conditional_t, IDimY, void>>; + // using spatial_dim_for = std::conditional_t, GridX, std::conditional_t, GridY, void>>; /** - * @brief An alias for the spatial discrete domain type. + * @brief An alias for the spatial discrete index range type. */ - using SpatialIdxRange = IDomainXY; + using SpatialIdxRange = IdxRangeXY; /** - * @brief An alias for the velocity discrete domain type. + * @brief An alias for the velocity discrete index range type. */ - using VelocityIdxRange = IDomainVxVy; + using VelocityIdxRange = IdxRangeVxVy; /** - * @brief An alias for the whole distribution function discrete domain type. + * @brief An alias for the whole distribution function discrete index range type. */ - using FdistribuIdxRange = IDomainSpXYVxVy; + using FdistribuIdxRange = IdxRangeSpXYVxVy; }; diff --git a/src/geometryXYVxVy/initialization/iequilibrium.hpp b/src/geometryXYVxVy/initialization/iequilibrium.hpp index e8fe48228..ab4f67978 100644 --- a/src/geometryXYVxVy/initialization/iequilibrium.hpp +++ b/src/geometryXYVxVy/initialization/iequilibrium.hpp @@ -4,10 +4,22 @@ #include +/** + * @brief An abstract class for initializing the equilibrium state of the distribution function. + * The equilibrium state does not depend on spatial dimensions. + */ class IEquilibrium { public: virtual ~IEquilibrium() = default; - virtual DSpanSpVxVy operator()(DSpanSpVxVy allfequilibrium) const = 0; + /** + * @brief Operator for initializing a distribution function that does not depend on space. + * + * @param[in, out] allfequilibrium On input: the uninitialized distribution function. + * On output: the initialized distribution function. + * @return The initialized distribution function. + */ + + virtual DFieldSpVxVy operator()(DFieldSpVxVy allfequilibrium) const = 0; }; diff --git a/src/geometryXYVxVy/initialization/iinitialization.hpp b/src/geometryXYVxVy/initialization/iinitialization.hpp index 078139f54..8a0ac7882 100644 --- a/src/geometryXYVxVy/initialization/iinitialization.hpp +++ b/src/geometryXYVxVy/initialization/iinitialization.hpp @@ -4,10 +4,21 @@ #include +/** + * @brief An abstract class that allows for initializing a distribution function. + */ class IInitialization { public: virtual ~IInitialization() = default; - virtual DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu) const = 0; + /** + * @brief Operator for initializing a distribution function. + * + * @param[in, out] allfdistribu On input: the uninitialized distribution function. + * On output: the initialized distribution function. + * @return The initialized distribution function. + */ + + virtual DFieldSpXYVxVy operator()(DFieldSpXYVxVy allfdistribu) const = 0; }; diff --git a/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp b/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp index 21cf031ef..d4f1f1a4c 100644 --- a/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp +++ b/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp @@ -14,14 +14,14 @@ MaxwellianEquilibrium::MaxwellianEquilibrium( { } -DSpanSpVxVy MaxwellianEquilibrium::operator()(DSpanSpVxVy const allfequilibrium) const +DFieldSpVxVy MaxwellianEquilibrium::operator()(DFieldSpVxVy const allfequilibrium) const { - IdxRangeSp const gridsp = allfequilibrium.domain(); - IDomainVxVy const gridvxvy = allfequilibrium.domain(); + IdxRangeSp const gridsp = get_idx_range(allfequilibrium); + IdxRangeVxVy const gridvxvy = get_idx_range(allfequilibrium); // Initialization of the maxwellian - DFieldVxVy maxwellian_alloc(gridvxvy); - DSpanVxVy maxwellian = maxwellian_alloc.span_view(); + DFieldMemVxVy maxwellian_alloc(gridvxvy); + DFieldVxVy maxwellian = get_field(maxwellian_alloc); ddc::for_each(gridsp, [&](IdxSp const isp) { compute_maxwellian( maxwellian, @@ -32,7 +32,7 @@ DSpanSpVxVy MaxwellianEquilibrium::operator()(DSpanSpVxVy const allfequilibrium) ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridvxvy, - KOKKOS_LAMBDA(IndexVxVy const ivxvy) { + KOKKOS_LAMBDA(IdxVxVy const ivxvy) { allfequilibrium(isp, ivxvy) = maxwellian(ivxvy); }); }); @@ -69,20 +69,20 @@ MaxwellianEquilibrium MaxwellianEquilibrium::init_from_input( with n the density and T the temperature and */ void MaxwellianEquilibrium::compute_maxwellian( - DSpanVxVy const fMaxwellian, + DFieldVxVy const fMaxwellian, double const density, double const temperature, double const mean_velocity) { double const inv_2pi = 1. / (2. * M_PI * temperature); - IDomainVxVy const gridvxvy = fMaxwellian.domain(); + IdxRangeVxVy const gridvxvy = get_idx_range(fMaxwellian); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridvxvy, - KOKKOS_LAMBDA(IndexVxVy const ivxvy) { - double const vx = ddc::coordinate(ddc::select(ivxvy)); - double const vy = ddc::coordinate(ddc::select(ivxvy)); + KOKKOS_LAMBDA(IdxVxVy const ivxvy) { + double const vx = ddc::coordinate(ddc::select(ivxvy)); + double const vy = ddc::coordinate(ddc::select(ivxvy)); fMaxwellian(ivxvy) = density * inv_2pi * Kokkos::exp( -((vx - mean_velocity) * (vx - mean_velocity) diff --git a/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp b/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp index a7a410300..83f871501 100644 --- a/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp +++ b/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp @@ -50,7 +50,7 @@ class MaxwellianEquilibrium : public IEquilibrium * @param[out] allfequilibrium A Span containing a Maxwellian distribution function. * @return A Span containing a Maxwellian distribution function. */ - DSpanSpVxVy operator()(DSpanSpVxVy allfequilibrium) const override; + DFieldSpVxVy operator()(DFieldSpVxVy allfequilibrium) const override; /** @@ -65,7 +65,7 @@ class MaxwellianEquilibrium : public IEquilibrium * @param[in] mean_velocity A parameter that represents the mean velocity of Maxwellian. */ static void compute_maxwellian( - DSpanVxVy const fMaxwellian, + DFieldVxVy const fMaxwellian, double const density, double const temperature, double const mean_velocity); diff --git a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp index a70016a79..fd6df873d 100644 --- a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp +++ b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp @@ -7,7 +7,7 @@ SingleModePerturbInitialization::SingleModePerturbInitialization( - DViewSpVxVy fequilibrium, + DConstFieldSpVxVy fequilibrium, host_t init_perturb_mode, host_t init_perturb_amplitude) : m_fequilibrium(fequilibrium) @@ -16,16 +16,16 @@ SingleModePerturbInitialization::SingleModePerturbInitialization( { } -DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const allfdistribu) const +DFieldSpXYVxVy SingleModePerturbInitialization::operator()(DFieldSpXYVxVy const allfdistribu) const { - IdxRangeSp const gridsp = allfdistribu.domain(); - IDomainXY const gridxy = allfdistribu.domain(); - IDomainXYVxVy const gridxyvxvy = allfdistribu.domain(); + IdxRangeSp const gridsp = get_idx_range(allfdistribu); + IdxRangeXY const gridxy = get_idx_range(allfdistribu); + IdxRangeXYVxVy const gridxyvxvy = get_idx_range(allfdistribu); // Initialization of the perturbation - DFieldXY perturbation_alloc(gridxy); - ddc::ChunkSpan fequilibrium_proxy = m_fequilibrium.span_view(); - ddc::ChunkSpan perturbation_proxy = perturbation_alloc.span_view(); + DFieldMemXY perturbation_alloc(gridxy); + ddc::ChunkSpan fequilibrium_proxy = get_field(m_fequilibrium); + ddc::ChunkSpan perturbation_proxy = get_field(perturbation_alloc); ddc::for_each(gridsp, [&](IdxSp const isp) { perturbation_initialization( perturbation_proxy, @@ -36,11 +36,11 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridxyvxvy, - KOKKOS_LAMBDA(IndexXYVxVy const ixyvxvy) { - IndexX const ix = ddc::select(ixyvxvy); - IndexY const iy = ddc::select(ixyvxvy); - IndexVx const ivx = ddc::select(ixyvxvy); - IndexVy const ivy = ddc::select(ixyvxvy); + KOKKOS_LAMBDA(IdxXYVxVy const ixyvxvy) { + IdxX const ix = ddc::select(ixyvxvy); + IdxY const iy = ddc::select(ixyvxvy); + IdxVx const ivx = ddc::select(ixyvxvy); + IdxVy const ivy = ddc::select(ixyvxvy); double fdistribu_val = fequilibrium_proxy(isp, ivx, ivy) * (1. + perturbation_proxy(ix, iy)); if (fdistribu_val < 1.e-60) { @@ -54,7 +54,7 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al SingleModePerturbInitialization SingleModePerturbInitialization::init_from_input( - DViewSpVxVy allfequilibrium, + DConstFieldSpVxVy allfequilibrium, IdxRangeSp dom_kinsp, PC_tree_t const& yaml_input_file) { @@ -76,23 +76,23 @@ SingleModePerturbInitialization SingleModePerturbInitialization::init_from_input void SingleModePerturbInitialization::perturbation_initialization( - DSpanXY const perturbation, + DFieldXY const perturbation, int const perturb_mode, double const perturb_amplitude) const { - IDomainXY const gridxy = perturbation.domain(); + IdxRangeXY const gridxy = get_idx_range(perturbation); double const kx = perturb_mode * 2. * M_PI - / ddcHelper::total_interval_length(ddc::select(gridxy)); + / ddcHelper::total_interval_length(ddc::select(gridxy)); double const ky = perturb_mode * 2. * M_PI - / ddcHelper::total_interval_length(ddc::select(gridxy)); + / ddcHelper::total_interval_length(ddc::select(gridxy)); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridxy, - KOKKOS_LAMBDA(IndexXY const ixy) { - IndexX const ix = ddc::select(ixy); + KOKKOS_LAMBDA(IdxXY const ixy) { + IdxX const ix = ddc::select(ixy); CoordX const x = ddc::coordinate(ix); - IndexY const iy = ddc::select(ixy); + IdxY const iy = ddc::select(ixy); CoordY const y = ddc::coordinate(iy); perturbation(ix, iy) = perturb_amplitude * (Kokkos::cos(kx * x) + Kokkos::cos(ky * y)); diff --git a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp index afa14408c..4e080df3b 100644 --- a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp +++ b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp @@ -12,7 +12,7 @@ /// Initialization operator with a sinusoidal perturbation of a Maxwellian. This initializes all species. class SingleModePerturbInitialization : public IInitialization { - DViewSpVxVy m_fequilibrium; + DConstFieldSpVxVy m_fequilibrium; host_t m_init_perturb_mode; @@ -28,7 +28,7 @@ class SingleModePerturbInitialization : public IInitialization * @param[in] perturb_amplitude The amplitude of the perturbation. */ void perturbation_initialization( - DSpanXY perturbation, + DFieldXY perturbation, int const perturb_mode, double const perturb_amplitude) const; @@ -39,7 +39,7 @@ class SingleModePerturbInitialization : public IInitialization * @param[in] init_perturb_amplitude The perturbation amplitude. */ SingleModePerturbInitialization( - DViewSpVxVy fequilibrium, + DConstFieldSpVxVy fequilibrium, host_t init_perturb_mode, host_t init_perturb_amplitude); @@ -50,7 +50,7 @@ class SingleModePerturbInitialization : public IInitialization * @param[in, out] allfdistribu The initialized distribution function. * @return The initialized distribution function. */ - DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu) const override; + DFieldSpXYVxVy operator()(DFieldSpXYVxVy allfdistribu) const override; /** * @brief Read init_perturb_mode and init_perturb amplitude in a YAML input file @@ -61,7 +61,7 @@ class SingleModePerturbInitialization : public IInitialization * @return an instance of SingleModePerturbInitialization class. */ static SingleModePerturbInitialization init_from_input( - DViewSpVxVy allfequilibrium, + DConstFieldSpVxVy allfequilibrium, IdxRangeSp dom_kinsp, PC_tree_t const& yaml_input_file); }; \ No newline at end of file diff --git a/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp b/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp index 6fe3620a4..b620db738 100644 --- a/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp +++ b/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp @@ -7,38 +7,39 @@ #include "chargedensitycalculator.hpp" -ChargeDensityCalculator::ChargeDensityCalculator(DViewVxVy coeffs) : m_quadrature(coeffs) {} +ChargeDensityCalculator::ChargeDensityCalculator(DConstFieldVxVy coeffs) : m_quadrature(coeffs) {} -void ChargeDensityCalculator::operator()(DSpanXY rho, DViewSpXYVxVy allfdistribu) const +void ChargeDensityCalculator::operator()(DFieldXY rho, DConstFieldSpXYVxVy allfdistribu) const { Kokkos::Profiling::pushRegion("ChargeDensityCalculator"); - IdxRangeSp const kin_species_domain = allfdistribu.domain(); + IdxRangeSp const kin_species_idx_range = get_idx_range(allfdistribu); host_t const charges_host = ddc::host_discrete_space().charges(); - host_t const kinetic_charges_host = charges_host[allfdistribu.domain()]; + host_t const kinetic_charges_host + = charges_host[get_idx_range(allfdistribu)]; auto const kinetic_charges_alloc = create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), kinetic_charges_host); - ddc::ChunkSpan kinetic_charges = kinetic_charges_alloc.span_view(); + ddc::ChunkSpan kinetic_charges = get_field(kinetic_charges_alloc); m_quadrature( Kokkos::DefaultExecutionSpace(), rho, - KOKKOS_LAMBDA(IndexXYVxVy idx) { + KOKKOS_LAMBDA(IdxXYVxVy idx) { double sum = 0.0; - for (auto isp : kinetic_charges.domain()) { + for (auto isp : get_idx_range(kinetic_charges)) { sum += kinetic_charges(isp) * allfdistribu(isp, idx); } return sum; }); - IdxSp const last_kin_species = kin_species_domain.back(); - IdxSp const last_species = charges_host.domain().back(); + IdxSp const last_kin_species = kin_species_idx_range.back(); + IdxSp const last_species = get_idx_range(charges_host).back(); if (last_kin_species != last_species) { double chargedens_adiabspecies = double(charge(last_species)); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), - rho.domain(), - KOKKOS_LAMBDA(IndexXY ixy) { rho(ixy) += chargedens_adiabspecies; }); + get_idx_range(rho), + KOKKOS_LAMBDA(IdxXY ixy) { rho(ixy) += chargedens_adiabspecies; }); } Kokkos::Profiling::popRegion(); diff --git a/src/geometryXYVxVy/poisson/chargedensitycalculator.hpp b/src/geometryXYVxVy/poisson/chargedensitycalculator.hpp index 73be945ab..6ed224a1f 100644 --- a/src/geometryXYVxVy/poisson/chargedensitycalculator.hpp +++ b/src/geometryXYVxVy/poisson/chargedensitycalculator.hpp @@ -21,7 +21,7 @@ class ChargeDensityCalculator : public IChargeDensityCalculator { private: - Quadrature m_quadrature; + Quadrature m_quadrature; public: /** @@ -29,12 +29,12 @@ class ChargeDensityCalculator : public IChargeDensityCalculator * @param[in] coeffs * The coefficients of the quadrature. */ - explicit ChargeDensityCalculator(DViewVxVy coeffs); + explicit ChargeDensityCalculator(DConstFieldVxVy coeffs); /** * @brief Computes the charge density rho from the distribution function. * @param[in, out] rho * @param[in] allfdistribu */ - void operator()(DSpanXY rho, DViewSpXYVxVy allfdistribu) const final; + void operator()(DFieldXY rho, DConstFieldSpXYVxVy allfdistribu) const final; }; diff --git a/src/geometryXYVxVy/poisson/ichargedensitycalculator.hpp b/src/geometryXYVxVy/poisson/ichargedensitycalculator.hpp index 0d3e8ad6b..8b9005a77 100644 --- a/src/geometryXYVxVy/poisson/ichargedensitycalculator.hpp +++ b/src/geometryXYVxVy/poisson/ichargedensitycalculator.hpp @@ -28,5 +28,5 @@ class IChargeDensityCalculator * @param[out] rho The charge density. * @param[in] allfdistribu The distribution function. */ - virtual void operator()(DSpanXY rho, DViewSpXYVxVy allfdistribu) const = 0; + virtual void operator()(DFieldXY rho, DConstFieldSpXYVxVy allfdistribu) const = 0; }; diff --git a/src/geometryXYVxVy/poisson/iqnsolver.hpp b/src/geometryXYVxVy/poisson/iqnsolver.hpp index 43f961ba5..b123cfe4a 100644 --- a/src/geometryXYVxVy/poisson/iqnsolver.hpp +++ b/src/geometryXYVxVy/poisson/iqnsolver.hpp @@ -1,18 +1,19 @@ // SPDX-License-Identifier: MIT #pragma once - #include #include +#include "ddc_aliases.hpp" + /** * @brief An operator which solves the Quasi-Neutrality equation using a fast * Fourier transform. * * An operator which solves the Quasi-Neutrality equation: * @f$ - \frac{d^2 \phi}{dx^2} = \rho @f$ - * using a fast Fourier transform on a periodic domain. + * using a fast Fourier transform on a periodic index range. * This operator only works for equidistant points. */ class IQNSolver @@ -29,8 +30,8 @@ class IQNSolver * @param[in] allfdistribu The distribution function. */ virtual void operator()( - DSpanXY electrostatic_potential, - DSpanXY electric_field_x, - DSpanXY electric_field_y, - DViewSpXYVxVy allfdistribu) const = 0; + DFieldXY electrostatic_potential, + DFieldXY electric_field_x, + DFieldXY electric_field_y, + DConstFieldSpXYVxVy allfdistribu) const = 0; }; diff --git a/src/geometryXYVxVy/poisson/nullqnsolver.cpp b/src/geometryXYVxVy/poisson/nullqnsolver.cpp index 65833a549..c8ed2e746 100644 --- a/src/geometryXYVxVy/poisson/nullqnsolver.cpp +++ b/src/geometryXYVxVy/poisson/nullqnsolver.cpp @@ -2,7 +2,10 @@ #include "nullqnsolver.hpp" -void NullQNSolver::operator()(DSpanXY const, DSpanXY const, DSpanXY const, DViewSpXYVxVy const) - const +void NullQNSolver::operator()( + DFieldXY const, + DFieldXY const, + DFieldXY const, + DConstFieldSpXYVxVy const) const { } diff --git a/src/geometryXYVxVy/poisson/nullqnsolver.hpp b/src/geometryXYVxVy/poisson/nullqnsolver.hpp index a5dee6369..1c2fbfd07 100644 --- a/src/geometryXYVxVy/poisson/nullqnsolver.hpp +++ b/src/geometryXYVxVy/poisson/nullqnsolver.hpp @@ -25,8 +25,8 @@ class NullQNSolver : public IQNSolver * @param[in] allfdistribu The distribution function. */ void operator()( - DSpanXY electrostatic_potential, - DSpanXY electric_field_x, - DSpanXY electric_field_y, - DViewSpXYVxVy allfdistribu) const override; + DFieldXY electrostatic_potential, + DFieldXY electric_field_x, + DFieldXY electric_field_y, + DConstFieldSpXYVxVy allfdistribu) const override; }; diff --git a/src/geometryXYVxVy/poisson/qnsolver.cpp b/src/geometryXYVxVy/poisson/qnsolver.cpp index 539c19313..76e2514a5 100644 --- a/src/geometryXYVxVy/poisson/qnsolver.cpp +++ b/src/geometryXYVxVy/poisson/qnsolver.cpp @@ -19,25 +19,25 @@ QNSolver::QNSolver(PoissonSolver const& solve_poisson, IChargeDensityCalculator } void QNSolver::operator()( - DSpanXY const electrostatic_potential, - DSpanXY const electric_field_x, - DSpanXY const electric_field_y, - DViewSpXYVxVy const allfdistribu) const + DFieldXY const electrostatic_potential, + DFieldXY const electric_field_x, + DFieldXY const electric_field_y, + DConstFieldSpXYVxVy const allfdistribu) const { Kokkos::Profiling::pushRegion("QNSolver"); - assert((electrostatic_potential.domain() == ddc::get_domain(allfdistribu))); - IDomainXY const xy_dom = electrostatic_potential.domain(); + assert((get_idx_range(electrostatic_potential) == get_idx_range(allfdistribu))); + IdxRangeXY const xy_dom = get_idx_range(electrostatic_potential); // Compute the RHS of the Quasi-Neutrality equation. - DFieldXY rho(xy_dom); - DFieldVxVy contiguous_slice_vxvy(allfdistribu.domain()); + DFieldMemXY rho(xy_dom); + DFieldMemVxVy contiguous_slice_vxvy(get_idx_range(allfdistribu)); m_compute_rho(rho, allfdistribu); VectorFieldSpan< double, - IDomainXY, - NDTag, - typename DFieldXY::layout_type, + IdxRangeXY, + NDTag, + typename DFieldMemXY::layout_type, Kokkos::DefaultExecutionSpace::memory_space> electric_field(electric_field_x, electric_field_y); m_solve_poisson(electrostatic_potential, electric_field, rho); diff --git a/src/geometryXYVxVy/poisson/qnsolver.hpp b/src/geometryXYVxVy/poisson/qnsolver.hpp index c0acec17c..3f09dca98 100644 --- a/src/geometryXYVxVy/poisson/qnsolver.hpp +++ b/src/geometryXYVxVy/poisson/qnsolver.hpp @@ -1,8 +1,8 @@ // SPDX-License-Identifier: MIT #pragma once - #include "chargedensitycalculator.hpp" +#include "ddc_aliases.hpp" #include "ipoisson_solver.hpp" #include "iqnsolver.hpp" @@ -12,7 +12,7 @@ * * An operator which solves the Quasi-Neutrality equation: * @f$ - \frac{d^2 \phi}{dx^2} = \rho @f$ - * using a fast Fourier transform on a periodic domain. + * using a fast Fourier transform on a periodic index range. * This operator only works for equidistant points. * * The electric field, @f$ \frac{d \phi}{dx} @f$ is calculated using @@ -21,8 +21,8 @@ class QNSolver : public IQNSolver { using PoissonSolver = IPoissonSolver< - IDomainXY, - IDomainXY, + IdxRangeXY, + IdxRangeXY, std::experimental::layout_right, typename Kokkos::DefaultExecutionSpace::memory_space>; PoissonSolver const& m_solve_poisson; @@ -48,8 +48,8 @@ class QNSolver : public IQNSolver * @param[in] allfdistribu The distribution function. */ void operator()( - DSpanXY electrostatic_potential, - DSpanXY electric_field_x, - DSpanXY electric_field_y, - DViewSpXYVxVy allfdistribu) const override; + DFieldXY electrostatic_potential, + DFieldXY electric_field_x, + DFieldXY electric_field_y, + DConstFieldSpXYVxVy allfdistribu) const override; }; diff --git a/src/geometryXYVxVy/time_integration/CMakeLists.txt b/src/geometryXYVxVy/time_integration/CMakeLists.txt index abc78f85c..c0ebff827 100644 --- a/src/geometryXYVxVy/time_integration/CMakeLists.txt +++ b/src/geometryXYVxVy/time_integration/CMakeLists.txt @@ -15,6 +15,8 @@ target_link_libraries("time_integration_xyvxvy" gslx::poisson_xy gslx::speciesinfo gslx::vlasov_xyvxvy + gslx::utils + ) add_library("gslx::time_integration_xyvxvy" ALIAS "time_integration_xyvxvy") diff --git a/src/geometryXYVxVy/time_integration/itimesolver.hpp b/src/geometryXYVxVy/time_integration/itimesolver.hpp index b8eb6d13b..a08462523 100644 --- a/src/geometryXYVxVy/time_integration/itimesolver.hpp +++ b/src/geometryXYVxVy/time_integration/itimesolver.hpp @@ -4,11 +4,23 @@ #include +/** + * @brief An abstract class for solving a Vlasov-Poisson system of equations. + */ class ITimeSolver { public: virtual ~ITimeSolver() = default; - virtual DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu, double dt, int steps = 1) + /** + * @brief Solves the Vlasov-Poisson system. + * @param[in, out] allfdistribu On input : the initial value of the distribution function. + * On output : the value of the distribution function after solving + * the Vlasov-Poisson system a given number of iterations. + * @param[in] dt The timestep. + * @param[in] steps The number of iterations to be performed by the predictor-corrector. + * @return The distribution function after solving the system. + */ + virtual DFieldSpXYVxVy operator()(DFieldSpXYVxVy allfdistribu, double dt, int steps = 1) const = 0; }; diff --git a/src/geometryXYVxVy/time_integration/predcorr.cpp b/src/geometryXYVxVy/time_integration/predcorr.cpp index a764e9736..3f9b62859 100644 --- a/src/geometryXYVxVy/time_integration/predcorr.cpp +++ b/src/geometryXYVxVy/time_integration/predcorr.cpp @@ -16,29 +16,29 @@ PredCorr::PredCorr(IVlasovSolver const& vlasov_solver, IQNSolver const& poisson_ { } -DSpanSpXYVxVy PredCorr::operator()( - DSpanSpXYVxVy const allfdistribu, +DFieldSpXYVxVy PredCorr::operator()( + DFieldSpXYVxVy const allfdistribu, double const dt, int const steps) const { auto allfdistribu_host_alloc = ddc::create_mirror_view_and_copy(allfdistribu); - ddc::ChunkSpan allfdistribu_host = allfdistribu_host_alloc.span_view(); + ddc::ChunkSpan allfdistribu_host = get_field(allfdistribu_host_alloc); // electrostatic potential and electric field (depending only on x) - DFieldXY electrostatic_potential(allfdistribu.domain()); - DFieldXY electric_field_x(allfdistribu.domain()); - DFieldXY electric_field_y(allfdistribu.domain()); + DFieldMemXY electrostatic_potential(get_idx_range(allfdistribu)); + DFieldMemXY electric_field_x(get_idx_range(allfdistribu)); + DFieldMemXY electric_field_y(get_idx_range(allfdistribu)); - host_t electrostatic_potential_host(allfdistribu.domain()); + host_t electrostatic_potential_host(get_idx_range(allfdistribu)); // a 2D chunck of the same size as fdistribu - DFieldSpXYVxVy allfdistribu_half_t(allfdistribu.domain()); + DFieldMemSpXYVxVy allfdistribu_half_t(get_idx_range(allfdistribu)); m_poisson_solver( electrostatic_potential, electric_field_x, electric_field_y, - allfdistribu.span_cview()); + get_const_field(allfdistribu)); int iter = 0; for (; iter < steps; ++iter) { @@ -50,7 +50,7 @@ DSpanSpXYVxVy PredCorr::operator()( electrostatic_potential, electric_field_x, electric_field_y, - allfdistribu.span_cview()); + get_const_field(allfdistribu)); // copies necessary to PDI ddc::parallel_deepcopy(allfdistribu_host, allfdistribu); ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); diff --git a/src/geometryXYVxVy/time_integration/predcorr.hpp b/src/geometryXYVxVy/time_integration/predcorr.hpp index ea55cb2db..0cc1569bc 100644 --- a/src/geometryXYVxVy/time_integration/predcorr.hpp +++ b/src/geometryXYVxVy/time_integration/predcorr.hpp @@ -45,5 +45,5 @@ class PredCorr : public ITimeSolver * @param[in] steps The number of iterations to be performed by the predictor-corrector. * @return The distribution function after solving the system. */ - DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu, double dt, int steps = 1) const override; + DFieldSpXYVxVy operator()(DFieldSpXYVxVy allfdistribu, double dt, int steps = 1) const override; }; diff --git a/src/geometryXYVxVy/vlasov/ivlasovsolver.hpp b/src/geometryXYVxVy/vlasov/ivlasovsolver.hpp index d3d239761..b1a9d6e2a 100644 --- a/src/geometryXYVxVy/vlasov/ivlasovsolver.hpp +++ b/src/geometryXYVxVy/vlasov/ivlasovsolver.hpp @@ -4,14 +4,29 @@ #include +/** + * @brief An abstract class for solving a Vlasov equation. + */ class IVlasovSolver { public: virtual ~IVlasovSolver() = default; - virtual DSpanSpXYVxVy operator()( - DSpanSpXYVxVy allfdistribu, - DViewXY efield_x, - DViewXY efield_y, + /** + * @brief Solves a Vlasov equation on a timestep dt. + * + * @param[in, out] allfdistribu On input : the initial value of the distribution function. + * On output : the value of the distribution function after solving + * the Vlasov equation. + * @param[in] efield_x The electric field in the x direction computed at all spatial positions. + * @param[in] efield_y The electric field in the y direction computed at all spatial positions. + * @param[in] dt The timestep. + * + * @return The distribution function after solving the Vlasov equation. + */ + virtual DFieldSpXYVxVy operator()( + DFieldSpXYVxVy allfdistribu, + DConstFieldXY efield_x, + DConstFieldXY efield_y, double dt) const = 0; }; diff --git a/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp b/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp index 119bcf53f..d31396754 100644 --- a/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp +++ b/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp @@ -5,10 +5,10 @@ #include "splitvlasovsolver.hpp" SplitVlasovSolver::SplitVlasovSolver( - IAdvectionSpatial const& advec_x, - IAdvectionSpatial const& advec_y, - IAdvectionVelocity const& advec_vx, - IAdvectionVelocity const& advec_vy) + IAdvectionSpatial const& advec_x, + IAdvectionSpatial const& advec_y, + IAdvectionVelocity const& advec_vx, + IAdvectionVelocity const& advec_vy) : m_advec_x(advec_x) , m_advec_y(advec_y) , m_advec_vx(advec_vx) @@ -16,10 +16,10 @@ SplitVlasovSolver::SplitVlasovSolver( { } -DSpanSpXYVxVy SplitVlasovSolver::operator()( - DSpanSpXYVxVy const allfdistribu, - DViewXY const electric_field_x, - DViewXY const electric_field_y, +DFieldSpXYVxVy SplitVlasovSolver::operator()( + DFieldSpXYVxVy const allfdistribu, + DConstFieldXY const electric_field_x, + DConstFieldXY const electric_field_y, double const dt) const { m_advec_x(allfdistribu, dt / 2); diff --git a/src/geometryXYVxVy/vlasov/splitvlasovsolver.hpp b/src/geometryXYVxVy/vlasov/splitvlasovsolver.hpp index cb42f397c..49e708270 100644 --- a/src/geometryXYVxVy/vlasov/splitvlasovsolver.hpp +++ b/src/geometryXYVxVy/vlasov/splitvlasovsolver.hpp @@ -6,31 +6,63 @@ #include "ivlasovsolver.hpp" -template +template class IAdvectionSpatial; -template +template class IAdvectionVelocity; +/** + * @brief A class that solves a Vlasov equation using Strang's splitting. + * + * The Vlasov equation is split between four advection equations + * along the X, Y, Vx and Vy directions. The splitting involves solving + * the advections in the X, Y, and Vx directions first on a time interval + * of length dt/2, then the Vy-direction advection on a time dt, and + * finally the X, Y, and Vx directions again in reverse order on dt/2. + */ class SplitVlasovSolver : public IVlasovSolver { - IAdvectionSpatial const& m_advec_x; - IAdvectionSpatial const& m_advec_y; + /// Advection operator in the x direction + IAdvectionSpatial const& m_advec_x; + /// Advection operator in the y direction + IAdvectionSpatial const& m_advec_y; - IAdvectionVelocity const& m_advec_vx; - IAdvectionVelocity const& m_advec_vy; + /// Advection operator in the vx direction + IAdvectionVelocity const& m_advec_vx; + /// Advection operator in the vy direction + IAdvectionVelocity const& m_advec_vy; public: + /** + * @brief Creates an instance of the split vlasov solver class. + * @param[in] advec_x An advection operator along the x direction. + * @param[in] advec_y An advection operator along the y direction. + * @param[in] advec_vx An advection operator along the vx direction. + * @param[in] advec_vy An advection operator along the vy direction. + */ SplitVlasovSolver( - IAdvectionSpatial const& advec_x, - IAdvectionSpatial const& advec_y, - IAdvectionVelocity const& advec_vx, - IAdvectionVelocity const& advec_vy); + IAdvectionSpatial const& advec_x, + IAdvectionSpatial const& advec_y, + IAdvectionVelocity const& advec_vx, + IAdvectionVelocity const& advec_vy); ~SplitVlasovSolver() override = default; - DSpanSpXYVxVy operator()( - DSpanSpXYVxVy allfdistribu, - DViewXY electric_field_x, - DViewXY electric_field_y, + /** + * @brief Solves a Vlasov equation on a timestep dt. + * + * @param[in, out] allfdistribu On input : the initial value of the distribution function. + * On output : the value of the distribution function after solving + * the Vlasov equation. + * @param[in] electric_field_x The electric field in the x direction computed at all spatial positions. + * @param[in] electric_field_y The electric field in the y direction computed at all spatial positions. + * @param[in] dt The timestep. + * + * @return The distribution function after solving the Vlasov equation. + */ + DFieldSpXYVxVy operator()( + DFieldSpXYVxVy allfdistribu, + DConstFieldXY electric_field_x, + DConstFieldXY electric_field_y, double dt) const override; }; diff --git a/src/multipatch/CMakeLists.txt b/src/multipatch/CMakeLists.txt index fca360501..8c03fbe8b 100644 --- a/src/multipatch/CMakeLists.txt +++ b/src/multipatch/CMakeLists.txt @@ -1,6 +1,5 @@ add_subdirectory(connectivity) -add_subdirectory(interfaces) add_subdirectory(spline) diff --git a/src/multipatch/README.md b/src/multipatch/README.md index 9da37255e..3d68e5da5 100644 --- a/src/multipatch/README.md +++ b/src/multipatch/README.md @@ -4,6 +4,4 @@ The `multipatch` folder contains all the code describing methods and classes whi - [connectivity](./connectivity/README.md) - Defines structures and methods to describe patches and how they are connected to one another. -- [interface](./interfaces/README.md) - Defines structures and methods for sticking patches together. - - [spline](./spline/README.md) - Defines structures and methods to deal with different splines on every patches. diff --git a/src/multipatch/connectivity/CMakeLists.txt b/src/multipatch/connectivity/CMakeLists.txt index 97008a0c2..94bddf490 100644 --- a/src/multipatch/connectivity/CMakeLists.txt +++ b/src/multipatch/connectivity/CMakeLists.txt @@ -7,5 +7,6 @@ target_include_directories("multipatch_connectivity" ) target_link_libraries("multipatch_connectivity" INTERFACE DDC::DDC + gslx::utils ) add_library("gslx::multipatch_connectivity" ALIAS "multipatch_connectivity") diff --git a/src/multipatch/connectivity/README.md b/src/multipatch/connectivity/README.md index 7da162a5c..78f66265b 100644 --- a/src/multipatch/connectivity/README.md +++ b/src/multipatch/connectivity/README.md @@ -1,5 +1,7 @@ # Multipatch connectivity +This directory defines structures and methods to describe patches and how they are connected to one another. + ## Patch tag The tag Patch refers to a single 2D patch geometry. The tag contains aliases (or shortcuts) to the DDC geometry elements, such as: @@ -12,4 +14,106 @@ The tag Patch refers to a single 2D patch geometry. The tag contains aliases (or * The type which describes the index range of a grid (e.g. `IdxRange1`). * The type which describes the index range of a spline coefficients grid (e.g. `BSIdxRange1`). -The domain defined on this patch is called *logical domain*. \ No newline at end of file +The domain defined on this patch is called *logical domain*. + + +## Interfaces +### Sticking of Two Edges + +We follow the idea to represent the topology of a multipatch domain using +the idea of *domain manifolds* (see [1]). This means that we have several +independent tensor-product logical domains and define their sticking via coordinate +transformations. + + +**Remark:** In the referenced paper, the authors use affine isometries of the entire +patch instead of coordinate transformations of the faces. This is equivalent to our approach in +the sense that, when the scaling of the domains is removed and we consider only +unit squares as logical domains, the coordinate transformations of the faces and the information +about which edges are identified determine the +isometries used in the paper and vice versa. This is easy to show. + +#### Multi-patch domain + +For simplicity, we will constrain ourselves to the 2D case, but this approach +could be generalized to arbitrary dimensions. + +Let $\Omega$ be the domain of interest and assume that we have patches $\Omega^{(i)}$, $i=1,...,K$, which are disjoint, s.t. +```math +\Omega = \dot{\bigcup}_{i=1}^K \Omega^{(i)}. +``` + +For every patch $\Omega^{(i)}$ we have a mapping +```math +F^{(i)}: [a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}] \rightarrow \Omega^{(i)}, +``` + +where the rectangular domain $`[a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}]`$ defines +the logical coordinates for the patch. We call this logical coordinate domain +the *logical patch*. + +#### Edges + +So we obtain four logical edges +```math +[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}, \quad +[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}, \quad +\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}], \quad +\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}] +``` + +In the code, we define edges as follows. +Every edge of a logical patch is identified via the patch it belongs to, the dimension +and whether it is at the front or the back of the domain. +So e.g.the edge $`[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}`$ would be identified with patch $i$, dimensions `RDimXi` and `FRONT`. +$`[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}`$ would be identified with patch $i$, dimensions `RDimXi` and `BACK` and +$`\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$ would be identified with patch $i$, dimensions `RDimYi` and `BACK`. + +#### Sticking and Coordinate Transformation + +Any edge can then be associated with an edge on another patch. +This corresponds to the 'sticking'. +So for a different patch $`\Omega^{(j)}`$, we have the logical coordinate domain +$`[a_x^{(j)}, b_x^{(j)}] \times [a_y^{(j)} b_y^{(j)}]`$. +If we want to stick patch $i$ to patch $j$, we have to determine the edges that +are identified and how they are identified. +The way that they are identified is mathematically determined via the coordinate transformation from one edge +to the other. +Since the transformations are supposed to be affine and bijective, there are only two options: +the transformation can be order preserving or +order reversing (this corresponds to the orientation of the phyiscal edge where two parametrizations +coming from the two patches can have either the same or the opposite orientation respectively). + +So for example, we want to stick the edge $`\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$ +on patch $i$ to the edge $`[a_x^{(j)}, b_x^{(j)}] \times \{ b_y^{(j)} \}`$ on patch $j$. +If the transformation is order-preserving (i.e. the orientations of the parametrizations +of the physical edge agree), then the transformation from the first edge to the second is +```math +t \mapsto a_x^{(j)} + \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}). +``` + +If the transformation is order-reversing (i.e. the orientations of the parametrizations +of the physical edge are opposite), then it is +```math +t \mapsto b_x^{(j)} - \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}). +``` + + +In the code, the sticking of two edges is represented by an `Interface` structure which contains tags +to the first patch and the second patch as well as the boolean member `orientations_agree`. +The transformation from one edge to the other is done using the `EdgeTransformation` operator. + +## References +[1] Buchegger, F., Jüttler, B., Mantzaflaris, A., +*Adaptively refined multi-patch B-splines with enhanced smoothness*. +Applied Mathematics and Computation, +Volume 272, Part 1 +(2016). +https://www.sciencedirect.com/science/article/pii/S009630031500836X. + + +## Contents +- `edge_transformation.hpp`: Defines the `EdgeTransformation` operator to transform the coordinate from one edge to the other (see above). +- `edge.hpp`: Defines the `Edge` structure. +- `interface.hpp`: Defines the `Interface` structure. +- `patch.hpp`: Defines the `Patch` tag. diff --git a/src/multipatch/connectivity/edge.hpp b/src/multipatch/connectivity/edge.hpp new file mode 100644 index 000000000..849dec70c --- /dev/null +++ b/src/multipatch/connectivity/edge.hpp @@ -0,0 +1,45 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include + + +enum Extremity { FRONT, BACK }; + + +/** + * @brief Define an edge of a given patch. + * + * An edge is defined by a patch, a dimension and an extremity. + * For example, in the patch defined on logical domain + * @f$ [a_x, b_x]\times[a_y, b_y] @f$, + * + * * the edge IDimX, BACK refers to the set @f$ [a_x, b_x]\times\{b_y\} @f$, + * * and the edge IDimX, FRONT refers to the set @f$ [a_x, b_x]\times\{a_y\} @f$. + * + * @tparam Patch Patch where the edge is defined. + * @tparam Grid Discrete dimension where the edge is defined. + * @tparam extremity_val The BACK or FRONT value. +*/ +template +struct Edge +{ + /// @brief Patch where the edge is defined. + using associated_patch = Patch; + /// @brief Discrete dimension where the edge is defined. + using grid = Grid; + /// @brief Design if the edge is on the BACK or the FRONT of the other dimension. + static constexpr Extremity extremity = extremity_val; +}; + + +/** + * @brief Define an edge for the outside domain. + * OutsideEdge is a pseudo-edge outside the domain + * used to define interfaces between patches and the + * outside domaine. +*/ +struct OutsideEdge +{ +}; diff --git a/src/multipatch/connectivity/edge_transformation.hpp b/src/multipatch/connectivity/edge_transformation.hpp new file mode 100644 index 000000000..4ac53b1ad --- /dev/null +++ b/src/multipatch/connectivity/edge_transformation.hpp @@ -0,0 +1,110 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include + +#include + +#include "interface.hpp" + +/** + * @brief Transform a coordinate or an index from one edge to the one on the other edge. + * + * According to the orientation of the interface, we compute + * * if True, + * @f$ \\ t \mapsto min_2 + \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$ + * + * * if False, + * @f$ \\ t \mapsto max_2 - \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$ + * + * where @f$ min_i @f$ and @f$ max_i @f$ are the minimum and maximum + * coordinates of the edge @f$ i @f$. + * + * @tparam Interface The Interface type where we want to compute the transformation. +*/ +template +class EdgeTransformation +{ + using EdgeGrid1 = typename Interface::Edge1::grid; + using EdgeGrid2 = typename Interface::Edge2::grid; + + using EdgeDim1 = typename EdgeGrid1::continuous_dimension_type; + using EdgeDim2 = typename EdgeGrid2::continuous_dimension_type; + + using Patch1 = typename Interface::Patch1; + using Patch2 = typename Interface::Patch2; + + using IdxRangeEdge1 = ddc::DiscreteDomain; + using IdxRangeEdge2 = ddc::DiscreteDomain; + + IdxRangeEdge1 const& m_idx_range_patch_1; + IdxRangeEdge2 const& m_idx_range_patch_2; + + +public: + /** + * @brief Instantiate an EdgeTransformation. + * @param idx_range_patch_1 1D index range on the patch 1 of the interface. + * @param idx_range_patch_2 1D index range on the patch 2 of the interface. + */ + EdgeTransformation( + IdxRangeEdge1 const& idx_range_patch_1, + IdxRangeEdge2 const& idx_range_patch_2) + : m_idx_range_patch_1(idx_range_patch_1) + , m_idx_range_patch_2(idx_range_patch_2) {}; + + ~EdgeTransformation() = default; + + + /** + * @brief Transform a coordinate on the edge in the dimension + * of the current patch to the analogous coordinate on the target patch. + * + * @param current_coord + * A coordinate on the edge of the current patch. + * + * @tparam CurrentDim + * The current continuous dimension of the given coordinate coord. + * + * @return The analogous coordinate on the target patch. + */ + template + ddc::Coordinate< + std::conditional_t, EdgeDim2, EdgeDim1>> const + operator()(ddc::Coordinate const& current_coord) const + { + // The other continuous dimension + using ODim = std::conditional_t, EdgeDim2, EdgeDim1>; + + // The index range of CurrentDim + using CurrentIdxRange = std:: + conditional_t, IdxRangeEdge1, IdxRangeEdge2>; + // The index range of other dimension + using OIdxRange = std:: + conditional_t, IdxRangeEdge2, IdxRangeEdge1>; + + // Get min and length on each 1D domains: + CurrentIdxRange const current_idx_range( + ddc::DiscreteDomain< + EdgeGrid1, + EdgeGrid2>(m_idx_range_patch_1, m_idx_range_patch_2)); + OIdxRange const target_idx_range(ddc::DiscreteDomain< + EdgeGrid1, + EdgeGrid2>(m_idx_range_patch_1, m_idx_range_patch_2)); + + ddc::Coordinate const current_min = ddc::coordinate(current_idx_range.front()); + double const current_length = ddcHelper::total_interval_length(current_idx_range); + + ddc::Coordinate const target_min = ddc::coordinate(target_idx_range.front()); + double const target_length = ddcHelper::total_interval_length(target_idx_range); + + double rescale_x = (current_coord - current_min) / current_length * target_length; + + bool constexpr orientations_agree = Interface::orientations_agree; + if constexpr (!orientations_agree) { + rescale_x = target_length - rescale_x; + } + return target_min + rescale_x; // This has type ddc::Coordinate. + }; +}; diff --git a/src/multipatch/connectivity/interface.hpp b/src/multipatch/connectivity/interface.hpp new file mode 100644 index 000000000..4e2db9e1a --- /dev/null +++ b/src/multipatch/connectivity/interface.hpp @@ -0,0 +1,38 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include + +#include "edge.hpp" + +/** + * @brief Represent a simple sticking of two edges. + * + * @tparam Edge1Type Edge type defined on first patch. + * @tparam Edge2Type Edge type defined on second patch. + * @tparam orientations_agree_bool Boolean, + * if true, the parametrisations of the physical edge have the same orientation; + * else, the parametrisations of the physical edge have the opposite orientation. + * @see EdgeTransformation +*/ +template +struct Interface +{ + /// @brief Edge type of the first patch. + using Edge1 = Edge1Type; + /// @brief Edge type of the second patch. + using Edge2 = Edge2Type; + + /** + * If True, the parametrisations of the physical edge have the same orientation. + * If False, the parametrisations of the physical edge have the opposite orientation. + * (See EdgeTransformation). + */ + static constexpr bool orientations_agree = orientations_agree_bool; + + /// @brief Type of first patch. + using Patch1 = typename Edge1::associated_patch; + /// @brief Type of second patch. + using Patch2 = typename Edge2::associated_patch; +}; diff --git a/src/multipatch/interfaces/CMakeLists.txt b/src/multipatch/interfaces/CMakeLists.txt deleted file mode 100644 index bd9a65069..000000000 --- a/src/multipatch/interfaces/CMakeLists.txt +++ /dev/null @@ -1,12 +0,0 @@ - - -add_library("multipatch_interface" INTERFACE) -target_include_directories("multipatch_interface" - INTERFACE - "$" -) -target_link_libraries("multipatch_interface" INTERFACE - DDC::DDC - gslx::utils -) -add_library("gslx::multipatch_interface" ALIAS "multipatch_interface") diff --git a/src/multipatch/interfaces/README.md b/src/multipatch/interfaces/README.md deleted file mode 100644 index 133e0cd25..000000000 --- a/src/multipatch/interfaces/README.md +++ /dev/null @@ -1,106 +0,0 @@ -# Interfaces - -This directory defines structures and methods for the sticking of -patches of a multi-patch domain. - -## Sticking of Two Edges - -We follow the idea to represent the topology of a multipatch domain using -the idea of *domain manifolds* (see [1]). This means that we have several -independent tensor-product logical domains and define their sticking via coordinate -transformations. - - -**Remark:** In the referenced paper, the authors use affine isometries of the entire -patch instead of coordinate transformations of the faces. This is equivalent to our approach in -the sense that, when the scaling of the domains is removed and we consider only -unit squares as logical domains, the coordinate transformations of the faces and the information -about which edges are identified determine the -isometries used in the paper and vice versa. This is easy to show. - -### Multi-patch domain - -For simplicity, we will constrain ourselves to the 2D case, but this approach -could be generalized to arbitrary dimensions. - -Let $\Omega$ be the domain of interest and assume that we have patches -$\Omega^{(i)}$, $i=1,...,K$, which are disjoint, s.t. -$$ -\Omega = \dot{\bigcup}_{i=1}^K \Omega^{(i)}. -$$ -For every patch $\Omega^{(i)}$ we have a mapping -$$ -F^{(i)}: [a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}] \rightarrow \Omega^{(i)}, -$$ -where the rectangular domain $`[a_x^{(i)}, b_x^{(i)}]\times[a_y^{(i)} b_y^{(i)}]`$ defines -the logical coordinates for the patch. We call this logical coordinate domain -the *logical patch*. - -### Edges - -So we obtain four logical edges -$$ -[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}, \quad -[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}, \quad -\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}], \quad -\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}] -$$ - -In the code, we define edges as follows. Every edge of a -logical patch is identified via the patch it belongs to, the dimension -and whether it is at the front or the back of the domain. So e.g. -the edge $`[a_x^{(i)}, b_x^{(i)}] \times \{ a_y^{(i)} \}`$ would be -identified with patch $i$, dimensions `RDimXi` and `FRONT`. -$`[a_x^{(i)}, b_x^{(i)}] \times \{ b_y^{(i)} \}`$ would be -identified with patch $i$, dimensions `RDimXi` and `BACK` and -$`\{ b_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$ would be -identified with patch $i$, dimensions `RDimYi` and `BACK`. - -### Sticking and Coordinate Transformation - -Any edge can then be identified with an edge on another patch. This corresponds -to the 'sticking'. So for a different patch $`\Omega^{(j)}`$, we have the logical coordinate domain -$`[a_x^{(j)}, b_x^{(j)}] \times [a_y^{(j)} b_y^{(j)}]`$. -If we want to stick patch $i$ to patch $j$, we have to determine the edges that -are identified and how they are identified. The way that they are identified -is mathematically determined via the coordinate transformation from one edge -to the other. Since the transformations are supposed to be affine and bijective, -there are only two options: The transformation can be order preserving or -order reversing (this corresponds to the orientation of the phyiscal edge where -two parametrizations coming from the two patches can have either the same or the -opposite orientation respectively). - -So for example, we want to stick the edge $`\{ a_x^{(i)} \} \times [a_y^{(i)}, b_y^{(i)}]`$ -on patch $i$ to the edge $`[a_x^{(j)}, b_x^{(j)}] \times \{ b_y^{(j)} \}`$ on patch -$j$. If the transformation is order-preserving (i.e. the orientations of the parametrizations -of the physical edge agree), then the transformation from the first edge to the second is -$$ -t \mapsto a_x^{(j)} + \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}). -$$ -If the transformation is order-reversing (i.e. the orientations of the parametrizations -of the physical edge are opposite), then it is -$$ -t \mapsto b_x^{(j)} - \frac{t - a_y^{(i)}}{b_y^{(i)} - a_y^{(i)}} \, (b_x^{(j)} - a_x^{(j)}). -$$ - - -In the code, the sticking of two edges is represented by an `Interface` structure -which contains references to -the first patch and the second patch as well as the boolean member -`orientations_agree`. The transformation from one edge to the other -is done using the -`EdgeCoordinatesTransformation` operator. - -## References -[1] Buchegger, F., Jüttler, B., Mantzaflaris, A., -*Adaptively refined multi-patch B-splines with enhanced smoothness*. -Applied Mathematics and Computation, -Volume 272, Part 1 -(2016). -https://www.sciencedirect.com/science/article/pii/S009630031500836X. - - -## Contents -- `coord_transformation.hpp`: Defines the `EdgeCoordinatesTransformation` operator to transform the coordinate from one edge to the other (see above). -- `edge.hpp`: Defines the `Edge` structure. -- `interface.hpp`: Defines the `Interface` structure. diff --git a/src/multipatch/interfaces/coord_transformation.hpp b/src/multipatch/interfaces/coord_transformation.hpp deleted file mode 100644 index da39f8743..000000000 --- a/src/multipatch/interfaces/coord_transformation.hpp +++ /dev/null @@ -1,114 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once - -#include - -#include - -#include "interface.hpp" - -/** - * @brief Transform a coordinate from one edge to a coordinate on the other edge. - * - * According to the orientation stored in the interface, we compute - * * if True, - * - * @f$ t \mapsto min_2 + \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$ - * - * * if False, - * - * @f$ t \mapsto max_2 - \frac{t - min_1}{max_1 - min_1}(max_2 - min_2) @f$ - * - * where @f$ min_i @f$ and @f$ max_i @f$ are the minimum and maximum - * coordinates of the edge @f$ i @f$. - * - * @tparam IDim1 - * The discrete dimension of the first edge. - * @tparam IDim2 - * The discrete dimension of the second edge. - * -*/ -template -class EdgeCoordinatesTransformation -{ - using RDim1 = typename IDim1::continuous_dimension_type; - using RDim2 = typename IDim2::continuous_dimension_type; - - Interface const m_interface; - -public: - /** - * @brief Instantiate an EdgeCoordinatesTransformation. - * - * @param interface - * An Interface object between the two patches. - */ - EdgeCoordinatesTransformation(Interface const interface) - : m_interface(interface) {}; - ~EdgeCoordinatesTransformation() = default; - - - /** - * @brief Transform a coordinate on the edge in the dimension - * of the current patch to the analogous coordinate on the target patch. - * - * @param current_coord - * A coordinate on the edge of the current patch. - * - * @tparam CurrentRDim - * The current continuous dimension of the given coordinate coord. - * - * @return The analogous coordinate on the target patch. - */ - template - ddc::Coordinate, RDim2, RDim1>> - operator()(ddc::Coordinate current_coord) const - { - // The discrete dimension of CurrentRDim - using CurrentIDim = std::conditional_t, IDim1, IDim2>; - // The other continuous dimension - using ORDim = std::conditional_t, RDim2, RDim1>; - // The other discrete dimension - using OIDim = std::conditional_t, IDim2, IDim1>; - - - bool const orientations_agree = m_interface.orientations_agree; - - ddc::Coordinate current_min; - double current_length; - get_min_and_length(current_min, current_length); - - ddc::Coordinate target_min; - double target_length; - get_min_and_length(target_min, target_length); - - double rescale_x = (current_coord - current_min) / current_length * target_length; - - if (!orientations_agree) { - rescale_x = target_length - rescale_x; - } - return target_min + rescale_x; // This has type ddc::Coordinate. - }; - - -private: - /** - * @brief Get the minimum coordinate of a patch on the edge and its length. - * - * @tparam IDim - * The discrete dimension of the edge of the current patch. - */ - template - void get_min_and_length( - ddc::Coordinate& min, - double& length) const - { - Edge const edge = m_interface.template get_edge(); - - ddc::DiscreteDomain const dom = edge.domain; - - min = ddc::coordinate(dom.front()); - length = ddcHelper::total_interval_length(dom); - }; -}; diff --git a/src/multipatch/interfaces/edge.hpp b/src/multipatch/interfaces/edge.hpp deleted file mode 100644 index 61f38e130..000000000 --- a/src/multipatch/interfaces/edge.hpp +++ /dev/null @@ -1,52 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once - -#include - - -enum Extremity { FRONT, BACK }; - - - -/** - * @brief Define an edge of a given patch. - * - * An edge is defined by a patch index, a dimension - * and an extremity. - * For example, in the patch defined on logical domain - * @f$ [a_x, b_x]\times[a_y, b_y] @f$, - * - * * the edge IDimX, BACK refers to the set @f$ [a_x, b_x]\times\{b_y\} @f$, - * * and the edge IDimX, FRONT refers to the set @f$ [a_x, b_x]\times\{a_y\} @f$. - * - * Here, the domain given as input corresponds to @f$ [a_x, b_x] @f$ - * in both case. - * - * - * @tparam IDim - * The discrete dimension of the edge. -*/ -template -struct Edge -{ -public: - /** - * The index of the patch. - */ - std::size_t patch_index; - /** - * The discrete domain of the edge. - * See example with @f$ [a_x, b_x] @f$. - */ - ddc::DiscreteDomain domain; - /** - * A Extremity type containing "BACK" or "FRONT" tag. - */ - Extremity extremity; - - /** - * The discrete dimension of the edge. - */ - using discrete_dimension = IDim; -}; diff --git a/src/multipatch/interfaces/interface.hpp b/src/multipatch/interfaces/interface.hpp deleted file mode 100644 index 4d75e5b41..000000000 --- a/src/multipatch/interfaces/interface.hpp +++ /dev/null @@ -1,56 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once - -#include - -#include "edge.hpp" - -/** - * @brief Represent a simple sticking of two edges. - * - * @tparam IDim1 - * The discrete dimension of the first edge. - * @tparam IDim2 - * The discrete dimension of the second edge. - * - * @see EdgeCoordinatesTransformation -*/ -template -struct Interface -{ - /** - * The Edge of the first patch. - */ - Edge edge_1; - /** - * The Edge of the second patch. - */ - Edge edge_2; - /** - * If True, the parametrisations of the physical edge have the same orientation. - * If False, the parametrisations of the physical edge have the opposite orientation. - * (See EdgeCoordinatesTransformation). - */ - bool orientations_agree; - - /** - * @brief Get the edge on the given dimension. - * - * @tparam IDim - * The discrete dimension of the edge. - * - * @return The edge of the interface defined on the given dimension. - */ - template - constexpr Edge get_edge() const - { - static_assert(((std::is_same_v) || (std::is_same_v))); - - if constexpr (std::is_same_v) { - return edge_1; - } else { - return edge_2; - } - }; -}; diff --git a/tests/multipatch/CMakeLists.txt b/tests/multipatch/CMakeLists.txt index 3b574b8dc..856f13953 100644 --- a/tests/multipatch/CMakeLists.txt +++ b/tests/multipatch/CMakeLists.txt @@ -1,29 +1,7 @@ - - -set(test_name_gtest "interface") - -add_executable("${test_name_gtest}" - coord_transformation.cpp - ../main.cpp -) -target_link_libraries("${test_name_gtest}" -PUBLIC - DDC::DDC - GTest::gtest - GTest::gmock - - gslx::multipatch_connectivity - gslx::multipatch_geometries - gslx::multipatch_interface -) -target_compile_definitions("${test_name_gtest}" PUBLIC) -gtest_discover_tests("${test_name_gtest}" DISCOVERY_MODE PRE_TEST) - - - +add_subdirectory(connectivity) add_subdirectory(geometries) add_subdirectory(spline) diff --git a/tests/multipatch/connectivity/CMakeLists.txt b/tests/multipatch/connectivity/CMakeLists.txt new file mode 100644 index 000000000..9fc0336b5 --- /dev/null +++ b/tests/multipatch/connectivity/CMakeLists.txt @@ -0,0 +1,17 @@ + + +add_executable("coord_transformation_test" + coord_transformation.cpp + ../../main.cpp +) +target_link_libraries("coord_transformation_test" +PUBLIC + DDC::DDC + GTest::gtest + GTest::gmock + + gslx::multipatch_connectivity + gslx::multipatch_geometries +) +target_compile_definitions("coord_transformation_test" PUBLIC) +gtest_discover_tests("coord_transformation_test" DISCOVERY_MODE PRE_TEST) diff --git a/tests/multipatch/coord_transformation.cpp b/tests/multipatch/connectivity/coord_transformation.cpp similarity index 72% rename from tests/multipatch/coord_transformation.cpp rename to tests/multipatch/connectivity/coord_transformation.cpp index 655ac4608..4e829f64d 100644 --- a/tests/multipatch/coord_transformation.cpp +++ b/tests/multipatch/connectivity/coord_transformation.cpp @@ -5,8 +5,8 @@ #include #include "2patches_2d_non_periodic_uniform.hpp" -#include "coord_transformation.hpp" #include "edge.hpp" +#include "edge_transformation.hpp" #include "interface.hpp" #include "patch.hpp" @@ -66,15 +66,12 @@ class CoordinateTransformationTest : public ::testing::Test TEST_F(CoordinateTransformationTest, InvertedOrientation) { - // --- TEST 1 --- - // Creating interfaces ....................................................................... - Edge edge_1 = {1, idx_range_x1, BACK}; - Edge edge_2 = {2, idx_range_x2, FRONT}; - - Interface interface = {edge_1, edge_2, false}; + using EgdeX1B = Edge; + using EgdeX2F = Edge; + using Interface12 = Interface; // Coordinate transformation ................................................................. - EdgeCoordinatesTransformation coord_transformation(interface); + EdgeTransformation coord_transformation(idx_range_x1, idx_range_x2); Patch1::Coord1 test_coord_x1(1.5); Patch2::Coord1 test_coord_x2(coord_transformation(test_coord_x1)); @@ -86,15 +83,12 @@ TEST_F(CoordinateTransformationTest, InvertedOrientation) TEST_F(CoordinateTransformationTest, StickingDifferentDimensions) { - // --- TEST 2 --- - // Creating interfaces ....................................................................... - Edge edge_1 = {1, idx_range_x1, BACK}; - Edge edge_2 = {2, idx_range_y2, FRONT}; - - Interface interface = {edge_1, edge_2, true}; + using EgdeX1B = Edge; + using EgdeY2F = Edge; + using Interface12 = Interface; // Coordinate transformation ................................................................. - EdgeCoordinatesTransformation coord_transformation(interface); + EdgeTransformation coord_transformation(idx_range_x1, idx_range_y2); Patch1::Coord1 test_coord_x1(0.7); Patch2::Coord2 test_coord_y2(coord_transformation(test_coord_x1)); @@ -106,15 +100,12 @@ TEST_F(CoordinateTransformationTest, StickingDifferentDimensions) TEST_F(CoordinateTransformationTest, ReverseTransformation) { - // --- TEST 3 --- - // Creating interfaces ....................................................................... - Edge edge_1 = {1, idx_range_x1, BACK}; - Edge edge_2 = {2, idx_range_x2, FRONT}; - - Interface interface = {edge_1, edge_2, false}; + using EgdeX1B = Edge; + using EgdeX2F = Edge; + using Interface12 = Interface; // Coordinate transformation ................................................................. - EdgeCoordinatesTransformation coord_transformation(interface); + EdgeTransformation coord_transformation(idx_range_x1, idx_range_x2); Patch2::Coord1 test_coord_x2(1.5); Patch1::Coord1 test_coord_x1 = coord_transformation(test_coord_x2); diff --git a/vendor/sll/include/sll/polar_bsplines.hpp b/vendor/sll/include/sll/polar_bsplines.hpp index da4b2577d..9c5353ba3 100644 --- a/vendor/sll/include/sll/polar_bsplines.hpp +++ b/vendor/sll/include/sll/polar_bsplines.hpp @@ -14,7 +14,7 @@ * A class containing all information describing polar bsplines. * * Polar bsplines are 2D bsplines with a special treatment for the central singular point - * of a polar domain. At this singular point new bsplines are created which traverse the + * of a polar index range. At this singular point new bsplines are created which traverse the * singular point and ensure the desired continuity condition. * * @tparam BSplinesR The basis of radial bsplines from which the polar bsplines are constructed. @@ -74,25 +74,25 @@ class PolarBSplines * The type of a 2D index for the subset of the polar bsplines which can be expressed as a tensor * product of 1D bsplines. */ - using tensor_product_discrete_element_type = ddc::DiscreteElement; + using tensor_product_index_type = ddc::DiscreteElement; /** - * The type of the 2D domain for the subset of the polar bsplines which can be expressed as a tensor + * The type of the 2D idx_range for the subset of the polar bsplines which can be expressed as a tensor * product of 1D bsplines. */ - using tensor_product_discrete_domain_type = ddc::DiscreteDomain; + using tensor_product_idx_range_type = ddc::DiscreteDomain; /** * The type of a 2D vector for the subset of the polar bsplines which can be expressed as a tensor * product of 1D bsplines. */ - using tensor_product_discrete_vector_type = ddc::DiscreteVector; + using tensor_product_idx_step_type = ddc::DiscreteVector; private: - using IndexR = ddc::DiscreteElement; - using IndexP = ddc::DiscreteElement; - using LengthR = ddc::DiscreteVector; - using LengthP = ddc::DiscreteVector; + using IdxR = ddc::DiscreteElement; + using IdxP = ddc::DiscreteElement; + using IdxStepR = ddc::DiscreteVector; + using IdxStepP = ddc::DiscreteVector; public: /** @@ -114,7 +114,7 @@ class PolarBSplines * the singular point. */ template - static constexpr ddc::DiscreteDomain singular_domain() + static constexpr ddc::DiscreteDomain singular_idx_range() { return ddc::DiscreteDomain( ddc::DiscreteElement {0}, @@ -127,11 +127,10 @@ class PolarBSplines * * @param idx The index of a 2D BSpline which is expressed as a tensor product of 1D BSplines. * - * @returns The index of the basis spline in the PolarBSpline domain. + * @returns The index of the basis spline in the PolarBSpline index range. */ template - static ddc::DiscreteElement get_polar_index( - tensor_product_discrete_element_type const& idx) + static ddc::DiscreteElement get_polar_index(tensor_product_index_type const& idx) { int const r_idx = ddc::select(idx).uid(); int const p_idx = ddc::select(idx).uid(); @@ -144,12 +143,12 @@ class PolarBSplines * Get the 2D index of the tensor product bspline which, when evaluated at the same point, * returns the same values as the polar bspline indicated by the index passed as an argument. * - * @param idx The index of the basis spline in the PolarBSpline domain. + * @param idx The index of the basis spline in the PolarBSpline index range. * * @returns The index of the equivalent 2D BSpline expressed as a 2D tensor product of 1D BSplines. */ template - static tensor_product_discrete_element_type get_2d_index(ddc::DiscreteElement const& idx) + static tensor_product_index_type get_2d_index(ddc::DiscreteElement const& idx) { assert(idx.uid() >= n_singular_basis()); int const idx_2d = idx.uid() - n_singular_basis(); @@ -161,12 +160,12 @@ class PolarBSplines } private: - using Spline2D = ddc::Chunk; + using Spline2D = ddc::Chunk; public: /** * The Impl class holds the implementation of the PolarBSplines. The implementation is specific to the - * memory space so that the Chunks can be defined with domains related to instances of this class. + * memory space so that the Chunks can be defined with index ranges related to instances of this class. * * @tparam MemorySpace Indicates where the object is saved. This is either on the host or the device. */ @@ -178,19 +177,19 @@ class PolarBSplines private: /** - * The type of the domain for the linear combinations defining the bsplines which traverse the singular point. + * The type of the index range for the linear combinations defining the bsplines which traverse the singular point. * * The bsplines which traverse the singular O-point are constructed from a linear combination of 2D * bsplines. These 2D bsplines can be expressed as a tensor product of 1D bsplines. This type - * describes the domain on which the coefficients of these linear combinations are defined. There is + * describes the index range on which the coefficients of these linear combinations are defined. There is * an index for the polar bspline being constructed, and 2 indices for the 2D bspline. */ - using singular_basis_linear_combination_domain_type + using singular_basis_linear_combination_idx_range_type = ddc::DiscreteDomain; ddc::Chunk< double, - singular_basis_linear_combination_domain_type, + singular_basis_linear_combination_idx_range_type, ddc::HostAllocator> m_singular_basis_elements; @@ -226,7 +225,7 @@ class PolarBSplines /// The type of an index associated with a PolarBSpline. using discrete_element_type = ddc::DiscreteElement; - /// The type of a domain of PolarBSplines. + /// The type of a index range of PolarBSplines. using discrete_domain_type = ddc::DiscreteDomain; /// The type of a vector associated with a PolarBSpline. @@ -244,7 +243,7 @@ class PolarBSplines { using DimX = typename DiscreteMapping::cartesian_tag_x; using DimY = typename DiscreteMapping::cartesian_tag_y; - using mapping_tensor_product_discrete_element_type = ddc::DiscreteElement< + using mapping_tensor_product_index_type = ddc::DiscreteElement< typename DiscreteMapping::BSplineR, typename DiscreteMapping::BSplineP>; if constexpr (C > -1) { @@ -256,7 +255,7 @@ class PolarBSplines for (std::size_t i(0); i < ddc::discrete_space().size(); ++i) { const ddc::Coordinate point = curvilinear_to_cartesian.control_point( - mapping_tensor_product_discrete_element_type(1, i)); + mapping_tensor_product_index_type(1, i)); const double c_x = ddc::get(point); const double c_y = ddc::get(point); @@ -288,53 +287,53 @@ class PolarBSplines ddc::init_discrete_space(barycentric_coordinate_converter); // The number of radial bases used to construct the bsplines traversing the singular point. - constexpr LengthR nr_in_singular(C + 1); + constexpr IdxStepR nr_in_singular(C + 1); assert(nr_in_singular.value() < int(ddc::discrete_space().size())); // The number of poloidal bases used to construct the bsplines traversing the singular point. - const LengthP np_in_singular(ddc::discrete_space().nbasis()); + const IdxStepP np_in_singular(ddc::discrete_space().nbasis()); // The number of elements of the poloidal basis which will have an associated coefficient // (This will be larger than np_in_singular as it includes the periodicity) - const LengthP np_tot(ddc::discrete_space().size()); + const IdxStepP np_tot(ddc::discrete_space().size()); - // The domain of the 2D bsplines in the innermost circles from which the polar bsplines + // The index range of the 2D bsplines in the innermost circles from which the polar bsplines // traversing the singular point will be constructed. - tensor_product_discrete_domain_type const dom_bsplines_inner( - tensor_product_discrete_element_type(0, 0), - tensor_product_discrete_vector_type(nr_in_singular, np_tot)); + tensor_product_idx_range_type const dom_bsplines_inner( + tensor_product_index_type(0, 0), + tensor_product_idx_step_type(nr_in_singular, np_tot)); // Initialise memory m_singular_basis_elements - = ddc::Chunk( - singular_basis_linear_combination_domain_type( - singular_domain(), + = ddc::Chunk( + singular_basis_linear_combination_idx_range_type( + singular_idx_range(), dom_bsplines_inner)); - ddc::DiscreteDomain bernstein_domain( + ddc::DiscreteDomain bernstein_idx_range( ddc::DiscreteElement {0}, ddc::DiscreteVector {n_singular_basis()}); - ddc::DiscreteDomain poloidal_spline_domain + ddc::DiscreteDomain poloidal_spline_idx_range = ddc::discrete_space().full_domain(); - for (IndexR const ir : ddc::DiscreteDomain(IndexR(0), LengthR(C + 1))) { - for (IndexP const ip : poloidal_spline_domain.take_first(np_in_singular)) { + for (IdxR const ir : ddc::DiscreteDomain(IdxR(0), IdxStepR(C + 1))) { + for (IdxP const ip : poloidal_spline_idx_range.take_first(np_in_singular)) { const ddc::Coordinate point = curvilinear_to_cartesian.control_point( - mapping_tensor_product_discrete_element_type(ir, ip)); + mapping_tensor_product_index_type(ir, ip)); ddc::Chunk> bernstein_vals( - bernstein_domain); + bernstein_idx_range); ddc::discrete_space().eval_basis(bernstein_vals, point); // Fill spline coefficients - for (auto k : bernstein_domain) { + for (auto k : bernstein_idx_range) { m_singular_basis_elements(discrete_element_type {k.uid()}, ir, ip) = bernstein_vals(k); } } - for (discrete_element_type k : singular_domain()) { - for (IndexP const ip : - poloidal_spline_domain.take_first(LengthP {BSplinesP::degree()})) { + for (discrete_element_type k : singular_idx_range()) { + for (IdxP const ip : + poloidal_spline_idx_range.take_first(IdxStepP {BSplinesP::degree()})) { m_singular_basis_elements(k, ir, ip + np_in_singular) = m_singular_basis_elements(k, ir, ip); } @@ -342,13 +341,13 @@ class PolarBSplines } } else { // Initialise m_singular_basis_elements to avoid any problems in the copy constructor - tensor_product_discrete_domain_type const empty_dom_bsplines( - tensor_product_discrete_element_type(0, 0), - tensor_product_discrete_vector_type(0, 0)); + tensor_product_idx_range_type const empty_dom_bsplines( + tensor_product_index_type(0, 0), + tensor_product_idx_step_type(0, 0)); m_singular_basis_elements - = ddc::Chunk( - singular_basis_linear_combination_domain_type( - singular_domain(), + = ddc::Chunk( + singular_basis_linear_combination_idx_range_type( + singular_idx_range(), empty_dom_bsplines)); } } @@ -420,7 +419,7 @@ class PolarBSplines * * @returns The 2D tensor product index of the first b-spline element in the values array. */ - tensor_product_discrete_element_type eval_basis( + tensor_product_index_type eval_basis( DSpan1D singular_values, DSpan2D values, ddc::Coordinate p) const; @@ -441,7 +440,7 @@ class PolarBSplines * * @returns The 2D tensor product index of the first b-spline element in the values array. */ - tensor_product_discrete_element_type eval_deriv_r( + tensor_product_index_type eval_deriv_r( DSpan1D singular_derivs, DSpan2D derivs, ddc::Coordinate p) const; @@ -462,7 +461,7 @@ class PolarBSplines * * @returns The 2D tensor product index of the first b-spline element in the values array. */ - tensor_product_discrete_element_type eval_deriv_p( + tensor_product_index_type eval_deriv_p( DSpan1D singular_derivs, DSpan2D derivs, ddc::Coordinate p) const; @@ -484,7 +483,7 @@ class PolarBSplines * * @returns The 2D tensor product index of the first b-spline element in the values array. */ - tensor_product_discrete_element_type eval_deriv_r_and_p( + tensor_product_index_type eval_deriv_r_and_p( DSpan1D singular_derivs, DSpan2D derivs, ddc::Coordinate p) const; @@ -509,9 +508,9 @@ class PolarBSplines } /** - * Returns the domain containing the indices of all the polar b-splines. + * Returns the index range containing the indices of all the polar b-splines. * - * @returns The domain containing the indices of all the polar b-splines. + * @returns The index range containing the indices of all the polar b-splines. */ discrete_domain_type full_domain() const noexcept { @@ -525,7 +524,7 @@ class PolarBSplines * @returns The ddc::DiscreteDomain containing the indices of the b-splines which don't traverse * the singular point. */ - discrete_domain_type tensor_bspline_domain() const noexcept + discrete_domain_type tensor_bspline_idx_range() const noexcept { return full_domain().remove_first(discrete_vector_type {n_singular_basis()}); } @@ -630,7 +629,7 @@ ddc::DiscreteElement PolarBSplines()) { + for (discrete_element_type k : singular_idx_range()) { singular_values(k.uid()) = 0.0; for (std::size_t i(0); i < nr_done; ++i) { for (std::size_t j(0); j < np; ++j) { @@ -685,14 +684,14 @@ void PolarBSplines::Impl::integrals( r_bspl_space.integrals(r_integrals); p_bspl_space.integrals(p_integrals); - ddc::for_each(singular_domain(), [&](auto k) { + ddc::for_each(singular_idx_range(), [&](auto k) { int_vals.singular_spline_coef(k) = ddc::transform_reduce( ddc::select(m_singular_basis_elements.domain()), 0.0, ddc::reducer::sum(), - [&](tensor_product_discrete_element_type const idx) { - IndexR i = ddc::select(idx); - IndexP j = ddc::select(idx); + [&](tensor_product_index_type const idx) { + IdxR i = ddc::select(idx); + IdxP j = ddc::select(idx); return m_singular_basis_elements(k, i, j) * r_integrals(i) * p_integrals(j); }); }); @@ -700,10 +699,10 @@ void PolarBSplines::Impl::integrals( ddc::DiscreteDomain r_tensor_product_dom( ddc::select(int_vals.spline_coef.domain())); - tensor_product_discrete_domain_type - tensor_bspline_domain(r_tensor_product_dom, p_integrals.domain()); + tensor_product_idx_range_type + tensor_bspline_idx_range(r_tensor_product_dom, p_integrals.domain()); - ddc::for_each(tensor_bspline_domain, [&](auto idx) { + ddc::for_each(tensor_bspline_idx_range, [&](auto idx) { int_vals.spline_coef(idx) = r_integrals(ddc::select(idx)) * p_integrals(ddc::select(idx)); }); @@ -711,7 +710,7 @@ void PolarBSplines::Impl::integrals( if (int_vals.spline_coef.domain().template extent() == p_bspl_space.size()) { ddc::DiscreteDomain periodic_points(p_integrals.domain().take_last( ddc::DiscreteVector {BSplinesP::degree()})); - tensor_product_discrete_domain_type repeat_domain(r_tensor_product_dom, periodic_points); - ddc::for_each(repeat_domain, [&](auto idx) { int_vals.spline_coef(idx) = 0.0; }); + tensor_product_idx_range_type repeat_idx_range(r_tensor_product_dom, periodic_points); + ddc::for_each(repeat_idx_range, [&](auto idx) { int_vals.spline_coef(idx) = 0.0; }); } } diff --git a/vendor/sll/tests/bernstein.cpp b/vendor/sll/tests/bernstein.cpp index 965169933..a1a0e4398 100644 --- a/vendor/sll/tests/bernstein.cpp +++ b/vendor/sll/tests/bernstein.cpp @@ -41,11 +41,11 @@ struct BernsteinFixture; template struct BernsteinFixture>> : public testing::Test { - struct DimX + struct X { static constexpr bool PERIODIC = false; }; - struct DimY + struct Y { static constexpr bool PERIODIC = false; }; @@ -59,7 +59,7 @@ struct BernsteinFixture>> : pu { }; static constexpr std::size_t poly_degree = D; - struct Bernstein : BernsteinPolynomialBasis + struct Bernstein : BernsteinPolynomialBasis { }; }; @@ -72,34 +72,34 @@ TYPED_TEST_SUITE(BernsteinFixture, Cases); TYPED_TEST(BernsteinFixture, PartitionOfUnity) { - using DimX = typename TestFixture::DimX; - using DimY = typename TestFixture::DimY; + using X = typename TestFixture::X; + using Y = typename TestFixture::Y; using Corner1 = typename TestFixture::Corner1; using Corner2 = typename TestFixture::Corner2; using Corner3 = typename TestFixture::Corner3; using Bernstein = typename TestFixture::Bernstein; - using CoordXY = ddc::Coordinate; + using CoordXY = ddc::Coordinate; const CoordXY c1(-1.0, -1.0); const CoordXY c2(0.0, 1.0); const CoordXY c3(1.0, -1.0); - CartesianToBarycentricCoordinates + CartesianToBarycentricCoordinates coordinate_converter(c1, c2, c3); ddc::init_discrete_space(coordinate_converter); - ddc::DiscreteDomain - domain(ddc::DiscreteElement(0), - ddc::DiscreteVector(Bernstein::nbasis())); + ddc::DiscreteDomain idx_range( + ddc::DiscreteElement(0), + ddc::DiscreteVector(Bernstein::nbasis())); - ddc::Chunk> values(domain); + ddc::Chunk> values(idx_range); std::size_t const n_test_points = 100; for (std::size_t i(0); i < n_test_points; ++i) { CoordXY const test_point = generate_random_point_in_triangle(c1, c2, c3); ddc::discrete_space().eval_basis(values, test_point); double total = ddc::transform_reduce( - domain, + idx_range, 0.0, ddc::reducer::sum(), [&](ddc::DiscreteElement const ix) { return values(ix); }); diff --git a/vendor/sll/tests/cosine_evaluator.hpp b/vendor/sll/tests/cosine_evaluator.hpp index 681761180..2a6d2e1c5 100644 --- a/vendor/sll/tests/cosine_evaluator.hpp +++ b/vendor/sll/tests/cosine_evaluator.hpp @@ -1,18 +1,26 @@ #pragma once - #include #include #include +#include "ddc_aliases.hpp" + +/** + * A wrapper around a cosine Evaluator that can be used in tests. + */ struct CosineEvaluator { - template + /** + * @brief An analytical evaluator to be used for exact comparisons in tests. + */ + template class Evaluator { public: - using Dim = DDim; + /// The grid dimension. + using Dim = Grid1D; private: static inline constexpr double s_2_pi = 2. * M_PI; @@ -23,43 +31,83 @@ struct CosineEvaluator double m_c1; public: - template - Evaluator(Domain domain) : m_c0(1.0) - , m_c1(0.0) + /** + * @brief A constructor taking the range over which the evaluator will be applied. + * This sets the equation to: + * @f$ cos(2 \pi x) @f$ + * + * @param[in] idx_range The range of the grid over which the evaluator will be applied. + */ + template + Evaluator(IdxRange idx_range) : m_c0(1.0) + , m_c1(0.0) { } + /** + * @brief A constructor parametrising the equation. + * This sets the equation to: + * @f$ cos(2 \pi (c_0 x + c_1)) @f$ + * + * @param[in] c0 The first parameter of the equation. + * @param[in] c1 The second parameter of the equation. + */ Evaluator(double c0, double c1) : m_c0(c0), m_c1(c1) {} + /** + * Evaluate the equation at x + * @param[in] x The position where the equation is evaluated. + * @returns The result of the equation. + */ double operator()(double const x) const noexcept { return eval(x, 0); } - void operator()(ddc::ChunkSpan> chunk) const + /** + * Evaluate the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + */ + void operator()(ddc::ChunkSpan> chunk) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - for (ddc::DiscreteElement const i : domain) { + for (ddc::DiscreteElement const i : idx_range) { chunk(i) = eval(ddc::coordinate(i), 0); } } + /** + * Evaluate the derivative of the equation at x + * @param[in] x The position where the equation is evaluated. + * @param[in] derivative The order of the derivative. + * @returns The result of the equation. + */ double deriv(double const x, int const derivative) const noexcept { return eval(x, derivative); } - void deriv(ddc::ChunkSpan> chunk, int const derivative) + /** + * Evaluate the derivative of the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + * @param[in] derivative The order of the derivative. + */ + void deriv(ddc::ChunkSpan> chunk, int const derivative) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - for (ddc::DiscreteElement const i : domain) { + for (ddc::DiscreteElement const i : idx_range) { chunk(i) = eval(ddc::coordinate(i), derivative); } } + /** + * The maximum norm of this equation. Used for normalisation. + * @param[in] diff The derivative whose max norm should be returned. + * @returns The maximum value of the infinite norm. + */ double max_norm(int diff = 0) const { return ipow(s_2_pi * m_c0, diff); diff --git a/vendor/sll/tests/evaluator_2d.hpp b/vendor/sll/tests/evaluator_2d.hpp index 1d3802643..ec18d9714 100644 --- a/vendor/sll/tests/evaluator_2d.hpp +++ b/vendor/sll/tests/evaluator_2d.hpp @@ -1,9 +1,16 @@ #pragma once - #include +#include "ddc_aliases.hpp" + +/** + * A wrapper around a 2D Evaluator that can be used in tests. + */ struct Evaluator2D { + /** + * @brief An analytical 2D evaluator combining 2 1D evaluators, to be used for exact comparisons in tests. + */ template class Evaluator { @@ -12,65 +19,113 @@ struct Evaluator2D Eval2 eval_func2; public: - template - Evaluator(Domain domain) - : eval_func1(ddc::select(domain)) - , eval_func2(ddc::select(domain)) + /** + * @brief A constructor taking the range over which the evaluator will be applied. + * The polynomial is initialised with random values between 0.0 and 1.0 + * + * @param[in] idx_range The range of the grid over which the evaluator will be applied. + */ + template + Evaluator(IdxRange idx_range) + : eval_func1(ddc::select(idx_range)) + , eval_func2(ddc::select(idx_range)) { } + /** + * Evaluate the equation at (x,y) + * @param[in] x The x-position where the equation is evaluated. + * @param[in] y The y-position where the equation is evaluated. + * @returns The result of the equation. + */ double operator()(double const x, double const y) const noexcept { return eval_func1(x) * eval_func2(y); } - template - double operator()(ddc::Coordinate const x) const noexcept + /** + * Evaluate the equation at @f$(x_1,x_2)@f$. + * @param[in] x The position where the equation is evaluated. + * @returns The result of the equation. + */ + template + double operator()(ddc::Coordinate const x) const noexcept { - return eval_func1(ddc::get(x)) * eval_func2(ddc::get(x)); + return eval_func1(ddc::get(x)) * eval_func2(ddc::get(x)); } - template - void operator()(ddc::ChunkSpan> chunk) const + /** + * Evaluate the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + */ + template + void operator()(ddc::ChunkSpan> chunk) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - ddc::for_each(domain, [=](ddc::DiscreteElement const i) { - chunk(i) = eval_func1(ddc::coordinate(ddc::select(i))) - * eval_func2(ddc::coordinate(ddc::select(i))); + ddc::for_each(idx_range, [=](ddc::DiscreteElement const i) { + chunk(i) = eval_func1(ddc::coordinate(ddc::select(i))) + * eval_func2(ddc::coordinate(ddc::select(i))); }); } + /** + * Evaluate the derivative of the equation at (x,y) + * @param[in] x The x-position where the equation is evaluated. + * @param[in] y The y-position where the equation is evaluated. + * @param[in] derivative_x The order of the x-derivative. + * @param[in] derivative_y The order of the x-derivative. + * @returns The result of the equation. + */ double deriv(double const x, double const y, int const derivative_x, int const derivative_y) const noexcept { return eval_func1.deriv(x, derivative_x) * eval_func2.deriv(y, derivative_y); } - template + /** + * Evaluate the derivative of the equation at @f$(x_1,x_2)@f$. + * @param[in] x The position where the equation is evaluated. + * @param[in] derivative_x The order of the x-derivative. + * @param[in] derivative_y The order of the x-derivative. + * @returns The result of the equation. + */ + template double deriv( - ddc::Coordinate const x, + ddc::Coordinate const x, int const derivative_x, int const derivative_y) const noexcept { - return eval_func1.deriv(ddc::get(x), derivative_x) - * eval_func2.deriv(ddc::get(x), derivative_y); + return eval_func1.deriv(ddc::get(x), derivative_x) + * eval_func2.deriv(ddc::get(x), derivative_y); } - template + /** + * Evaluate the derivative of the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + * @param[in] derivative_x The order of the x-derivative. + * @param[in] derivative_y The order of the x-derivative. + */ + template void deriv( - ddc::ChunkSpan> chunk, + ddc::ChunkSpan> chunk, int const derivative_x, int const derivative_y) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - ddc::for_each(domain, [=](ddc::DiscreteElement const i) { - chunk(i) = eval_func1.deriv(ddc::coordinate(ddc::select(i)), derivative_x) - * eval_func2.deriv(ddc::coordinate(ddc::select(i)), derivative_y); + ddc::for_each(idx_range, [=](ddc::DiscreteElement const i) { + chunk(i) = eval_func1.deriv(ddc::coordinate(ddc::select(i)), derivative_x) + * eval_func2.deriv(ddc::coordinate(ddc::select(i)), derivative_y); }); } + /** + * The maximum norm of this equation. Used for normalisation. + * @param[in] diff1 The x-derivative whose max norm should be returned. + * @param[in] diff2 The y-derivative whose max norm should be returned. + * @returns The maximum value of the infinite norm. + */ double max_norm(int diff1 = 0, int diff2 = 0) const { return eval_func1.max_norm(diff1) * eval_func2.max_norm(diff2); diff --git a/vendor/sll/tests/gauss_legendre_integration.cpp b/vendor/sll/tests/gauss_legendre_integration.cpp index 6a0695314..eb41185ec 100644 --- a/vendor/sll/tests/gauss_legendre_integration.cpp +++ b/vendor/sll/tests/gauss_legendre_integration.cpp @@ -26,7 +26,7 @@ class fn std::size_t m_n; }; -struct DimX +struct X { }; @@ -42,24 +42,24 @@ int test_integrate() std::stringstream oss; oss << std::scientific << std::hexfloat; - std::vector> domains + std::vector> idx_ranges = {{0.0, 1.0}, {1.0, 2.0}, {-0.2, 1.5}, {-1.5, -1.0}}; bool test_passed = true; - for (std::size_t order = 1; order <= GaussLegendre::max_order(); ++order) { - GaussLegendre const gl(order); + for (std::size_t order = 1; order <= GaussLegendre::max_order(); ++order) { + GaussLegendre const gl(order); std::cout << "integration at order " << order; std::cout << std::endl; for (std::size_t p = 0; p < 2 * order; ++p) { fn const f(p); - for (std::size_t i = 0; i < domains.size(); ++i) { + for (std::size_t i = 0; i < idx_ranges.size(); ++i) { double const sol_exact = 1.0 / (p + 1) - * (std::pow(domains[i].second, p + 1) - - std::pow(domains[i].first, p + 1)); - double const sol_num = gl.integrate(f, domains[i].first, domains[i].second); + * (std::pow(idx_ranges[i].second, p + 1) + - std::pow(idx_ranges[i].first, p + 1)); + double const sol_num = gl.integrate(f, idx_ranges[i].first, idx_ranges[i].second); double const err = std::fabs((sol_num - sol_exact) / sol_exact); bool ok = true; @@ -73,8 +73,8 @@ int test_integrate() oss << " of x^" << std::setw(2) << std::left << p; oss << ' '; oss << std::fixed << std::setprecision(1) << std::right; - oss << " on the domain [" << std::setw(4) << domains[i].first << ", " - << std::setw(4) << domains[i].second << "]"; + oss << " on the idx_range [" << std::setw(4) << idx_ranges[i].first << ", " + << std::setw(4) << idx_ranges[i].second << "]"; oss << std::scientific << std::hexfloat; oss << ' '; oss << std::setw(25) << std::left << sol_num; @@ -103,38 +103,38 @@ int test_compute_points_and_weights() std::stringstream oss; oss << std::scientific << std::hexfloat; - std::vector> domains + std::vector> idx_ranges = {{0.0, 1.0}, {1.0, 2.0}, {-0.2, 1.5}, {-1.5, -1.0}}; bool test_passed = true; - struct IDimX + struct GridX { }; - for (std::size_t order = 1; order <= GaussLegendre::max_order(); ++order) { - ddc::DiscreteDomain - domain(ddc::DiscreteElement {0}, ddc::DiscreteVector {order}); - ddc::Chunk, ddc::DiscreteDomain> gl_points(domain); - ddc::Chunk> gl_weights(domain); - GaussLegendre const gl(order); + for (std::size_t order = 1; order <= GaussLegendre::max_order(); ++order) { + ddc::DiscreteDomain + idx_range(ddc::DiscreteElement {0}, ddc::DiscreteVector {order}); + ddc::Chunk, ddc::DiscreteDomain> gl_points(idx_range); + ddc::Chunk> gl_weights(idx_range); + GaussLegendre const gl(order); std::cout << "integration at order " << order; std::cout << std::endl; for (std::size_t p = 0; p < 2 * order; ++p) { fn const f(p); - for (std::size_t i = 0; i < domains.size(); ++i) { + for (std::size_t i = 0; i < idx_ranges.size(); ++i) { gl.compute_points_and_weights( gl_points.span_view(), gl_weights.span_view(), - ddc::Coordinate(domains[i].first), - ddc::Coordinate(domains[i].second)); + ddc::Coordinate(idx_ranges[i].first), + ddc::Coordinate(idx_ranges[i].second)); double const sol_exact = 1.0 / (p + 1) - * (std::pow(domains[i].second, p + 1) - - std::pow(domains[i].first, p + 1)); + * (std::pow(idx_ranges[i].second, p + 1) + - std::pow(idx_ranges[i].first, p + 1)); double sol_num = 0.0; - for (auto xi : domain) { + for (auto xi : idx_range) { sol_num += gl_weights(xi) * f(gl_points(xi)); } double const err = std::fabs((sol_num - sol_exact) / sol_exact); @@ -150,8 +150,8 @@ int test_compute_points_and_weights() oss << " of x^" << std::setw(2) << std::left << p; oss << ' '; oss << std::fixed << std::setprecision(1) << std::right; - oss << " on the domain [" << std::setw(4) << domains[i].first << ", " - << std::setw(4) << domains[i].second << "]"; + oss << " on the idx_range [" << std::setw(4) << idx_ranges[i].first << ", " + << std::setw(4) << idx_ranges[i].second << "]"; oss << std::scientific << std::hexfloat; oss << ' '; oss << std::setw(25) << std::left << sol_num; diff --git a/vendor/sll/tests/mapping_jacobian.cpp b/vendor/sll/tests/mapping_jacobian.cpp index c1abf356e..ac5f9e631 100644 --- a/vendor/sll/tests/mapping_jacobian.cpp +++ b/vendor/sll/tests/mapping_jacobian.cpp @@ -14,32 +14,32 @@ namespace { -struct DimX +struct X { }; -struct DimY +struct Y { }; -struct DimR +struct R { static bool constexpr PERIODIC = false; }; -struct DimP +struct P { static bool constexpr PERIODIC = true; }; -using CoordR = ddc::Coordinate; -using CoordP = ddc::Coordinate; -using CoordRP = ddc::Coordinate; +using CoordR = ddc::Coordinate; +using CoordP = ddc::Coordinate

; +using CoordRP = ddc::Coordinate; int constexpr BSDegree = 3; -struct BSplinesR : ddc::NonUniformBSplines +struct BSplinesR : ddc::NonUniformBSplines { }; -struct BSplinesP : ddc::NonUniformBSplines +struct BSplinesP : ddc::NonUniformBSplines { }; @@ -48,10 +48,10 @@ using InterpPointsR = ddc:: using InterpPointsP = ddc:: GrevilleInterpolationPoints; -struct IDimR : InterpPointsR::interpolation_discrete_dimension_type +struct GridR : InterpPointsR::interpolation_discrete_dimension_type { }; -struct IDimP : InterpPointsP::interpolation_discrete_dimension_type +struct GridP : InterpPointsP::interpolation_discrete_dimension_type { }; @@ -60,51 +60,51 @@ using SplineRPBuilder = ddc::SplineBuilder2D< Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using SplineRPEvaluator = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::NullExtrapolationRule, ddc::NullExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; -using BSDomainR = ddc::DiscreteDomain; -using BSDomainP = ddc::DiscreteDomain; -using BSDomainRP = ddc::DiscreteDomain; +using BSIdxRangeR = ddc::DiscreteDomain; +using BSIdxRangeP = ddc::DiscreteDomain; +using BSIdxRangeRP = ddc::DiscreteDomain; -using IDomainR = ddc::DiscreteDomain; -using IDomainP = ddc::DiscreteDomain; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeR = ddc::DiscreteDomain; +using IdxRangeP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; -using IndexR = ddc::DiscreteElement; -using IndexP = ddc::DiscreteElement; -using IndexRP = ddc::DiscreteElement; +using IdxR = ddc::DiscreteElement; +using IdxP = ddc::DiscreteElement; +using IdxRP = ddc::DiscreteElement; -using IVectR = ddc::DiscreteVector; -using IVectP = ddc::DiscreteVector; -using IVectRP = ddc::DiscreteVector; +using IdxStepR = ddc::DiscreteVector; +using IdxStepP = ddc::DiscreteVector; +using IdxStepRP = ddc::DiscreteVector; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; template -using FieldRP = ddc::Chunk; +using FieldMemRP = ddc::Chunk; using Matrix_2x2 = std::array, 2>; @@ -149,35 +149,35 @@ class InverseJacobianMatrix : public testing::TestWithParam mapping; + const CircularToCartesian mapping; CoordR const r_min(0.0); CoordR const r_max(1.0); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point at r = 0. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point at r = 0. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 Jacobian_matrix; Matrix_2x2 inv_Jacobian_matrix; @@ -192,35 +192,35 @@ TEST_P(InverseJacobianMatrix, InverseMatrixCircMap) TEST_P(InverseJacobianMatrix, 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point at r = 0. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point at r = 0. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 Jacobian_matrix; Matrix_2x2 inv_Jacobian_matrix; @@ -235,15 +235,15 @@ TEST_P(InverseJacobianMatrix, InverseMatrixCzarMap) TEST_P(InverseJacobianMatrix, 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); @@ -262,29 +262,28 @@ TEST_P(InverseJacobianMatrix, InverseMatrixDiscCzarMap) ddc::init_discrete_space(r_knots); ddc::init_discrete_space(p_knots); - ddc::init_discrete_space(InterpPointsR::get_sampling()); - ddc::init_discrete_space(InterpPointsP::get_sampling()); + ddc::init_discrete_space(InterpPointsR::get_sampling()); + ddc::init_discrete_space(InterpPointsP::get_sampling()); - IDomainR interpolation_domain_R(InterpPointsR::get_domain()); - IDomainP interpolation_domain_P(InterpPointsP::get_domain()); - IDomainRP grid(interpolation_domain_R, interpolation_domain_P); + IdxRangeR interpolation_idx_range_R(InterpPointsR::get_domain()); + IdxRangeP interpolation_idx_range_P(InterpPointsP::get_domain()); + IdxRangeRP grid(interpolation_idx_range_R, interpolation_idx_range_P); SplineRPBuilder builder(grid); ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; + ddc::PeriodicExtrapolationRule

p_extrapolation_rule; SplineRPEvaluator evaluator( r_extrapolation_rule, r_extrapolation_rule, p_extrapolation_rule, p_extrapolation_rule); - DiscreteToCartesian mapping - = DiscreteToCartesian:: - analytical_to_discrete(analytical_mapping, builder, evaluator); + DiscreteToCartesian mapping = DiscreteToCartesian:: + analytical_to_discrete(analytical_mapping, builder, evaluator); // Test for each coordinates if the inv_Jacobian_matrix is the inverse of the Jacobian_matrix - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { const CoordRP coord_rp(ddc::coordinate(irp)); - const double r = ddc::get(coord_rp); + const double r = ddc::get(coord_rp); if (fabs(r) > 1e-15) { Matrix_2x2 Jacobian_matrix; Matrix_2x2 inv_Jacobian_matrix; diff --git a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp index 7c935bc60..d9b33dfc2 100644 --- a/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp +++ b/vendor/sll/tests/mapping_jacobian_matrix_coef.cpp @@ -14,32 +14,32 @@ namespace { -struct DimX +struct X { }; -struct DimY +struct Y { }; -struct DimR +struct R { static bool constexpr PERIODIC = false; }; -struct DimP +struct P { static bool constexpr PERIODIC = true; }; -using CoordR = ddc::Coordinate; -using CoordP = ddc::Coordinate; -using CoordRP = ddc::Coordinate; +using CoordR = ddc::Coordinate; +using CoordP = ddc::Coordinate

; +using CoordRP = ddc::Coordinate; int constexpr BSDegree = 3; -struct BSplinesR : ddc::NonUniformBSplines +struct BSplinesR : ddc::NonUniformBSplines { }; -struct BSplinesP : ddc::NonUniformBSplines +struct BSplinesP : ddc::NonUniformBSplines { }; @@ -48,10 +48,10 @@ using InterpPointsR = ddc:: using InterpPointsP = ddc:: GrevilleInterpolationPoints; -struct IDimR : InterpPointsR::interpolation_discrete_dimension_type +struct GridR : InterpPointsR::interpolation_discrete_dimension_type { }; -struct IDimP : InterpPointsP::interpolation_discrete_dimension_type +struct GridP : InterpPointsP::interpolation_discrete_dimension_type { }; @@ -60,51 +60,51 @@ using SplineRPBuilder = ddc::SplineBuilder2D< Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using SplineRPEvaluator = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::NullExtrapolationRule, ddc::NullExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; -using BSDomainR = ddc::DiscreteDomain; -using BSDomainP = ddc::DiscreteDomain; -using BSDomainRP = ddc::DiscreteDomain; +using BSIdxRangeR = ddc::DiscreteDomain; +using BSIdxRangeP = ddc::DiscreteDomain; +using BSIdxRangeRP = ddc::DiscreteDomain; -using IDomainR = ddc::DiscreteDomain; -using IDomainP = ddc::DiscreteDomain; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeR = ddc::DiscreteDomain; +using IdxRangeP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; -using IndexR = ddc::DiscreteElement; -using IndexP = ddc::DiscreteElement; -using IndexRP = ddc::DiscreteElement; +using IdxR = ddc::DiscreteElement; +using IdxP = ddc::DiscreteElement; +using IdxRP = ddc::DiscreteElement; -using IVectR = ddc::DiscreteVector; -using IVectP = ddc::DiscreteVector; -using IVectRP = ddc::DiscreteVector; +using IdxStepR = ddc::DiscreteVector; +using IdxStepP = ddc::DiscreteVector; +using IdxStepRP = ddc::DiscreteVector; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; template -using FieldRP = ddc::Chunk; +using FieldMemRP = ddc::Chunk; using Matrix_2x2 = std::array, 2>; @@ -147,37 +147,37 @@ 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point at r = 0. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point at r = 0. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the coefficients defined by the coefficients functions //are the same as the coefficients in the matrix function. // --- for the Jacobian matrix: - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 Jacobian_matrix; Matrix_2x2 Jacobian_matrix_coeff; @@ -191,7 +191,7 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCircMap) }); // --- for the inverse Jacobian matrix: - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 inv_Jacobian_matrix; Matrix_2x2 inv_Jacobian_matrix_coeff; @@ -211,37 +211,37 @@ 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point at r = 0. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point at r = 0. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the coefficients defined by the coefficients functions //are the same as the coefficients in the matrix function. // --- for the Jacobian matrix: - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 Jacobian_matrix; Matrix_2x2 Jacobian_matrix_coeff; @@ -255,7 +255,7 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixCzarMap) }); // --- for the inverseJacobian matrix: - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 inv_Jacobian_matrix; Matrix_2x2 inv_Jacobian_matrix_coeff; @@ -275,15 +275,15 @@ 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); @@ -302,30 +302,29 @@ TEST_P(JacobianMatrixAndJacobianCoefficients, MatrixDiscCzarMap) ddc::init_discrete_space(r_knots); ddc::init_discrete_space(p_knots); - ddc::init_discrete_space(InterpPointsR::get_sampling()); - ddc::init_discrete_space(InterpPointsP::get_sampling()); + ddc::init_discrete_space(InterpPointsR::get_sampling()); + ddc::init_discrete_space(InterpPointsP::get_sampling()); - IDomainR interpolation_domain_R(InterpPointsR::get_domain()); - IDomainP interpolation_domain_P(InterpPointsP::get_domain()); - IDomainRP grid(interpolation_domain_R, interpolation_domain_P); + IdxRangeR interpolation_idx_range_R(InterpPointsR::get_domain()); + IdxRangeP interpolation_idx_range_P(InterpPointsP::get_domain()); + IdxRangeRP grid(interpolation_idx_range_R, interpolation_idx_range_P); SplineRPBuilder builder(grid); ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; + ddc::PeriodicExtrapolationRule

p_extrapolation_rule; SplineRPEvaluator evaluator( r_extrapolation_rule, r_extrapolation_rule, p_extrapolation_rule, p_extrapolation_rule); - DiscreteToCartesian mapping - = DiscreteToCartesian:: - analytical_to_discrete(analytical_mapping, builder, evaluator); + DiscreteToCartesian mapping = DiscreteToCartesian:: + analytical_to_discrete(analytical_mapping, builder, evaluator); // Test for each coordinates if the coefficients defined by the coefficients functions //are the same as the coefficients in the matrix function. - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { const CoordRP coord_rp(ddc::coordinate(irp)); - const double r = ddc::get(coord_rp); + const double r = ddc::get(coord_rp); if (fabs(r) > 1e-15) { // --- for the Jacobian matrix: Matrix_2x2 Jacobian_matrix; diff --git a/vendor/sll/tests/matrix_batch_tridiag.cpp b/vendor/sll/tests/matrix_batch_tridiag.cpp index d24d7e85b..bbf02ecdc 100644 --- a/vendor/sll/tests/matrix_batch_tridiag.cpp +++ b/vendor/sll/tests/matrix_batch_tridiag.cpp @@ -2,7 +2,7 @@ #include "sll/matrix_batch_tridiag.hpp" -using View2d = Kokkos::View; +using ConstField2d = Kokkos::View; void solve_batched_tridiag_system( int const batch_size, @@ -13,10 +13,10 @@ void solve_batched_tridiag_system( Kokkos::View Rhs_view_host) { - View2d A_view("A", batch_size, mat_size); - View2d B_view("B", batch_size, mat_size); - View2d C_view("C", batch_size, mat_size); - View2d Rhs_view("R", batch_size, mat_size); + ConstField2d A_view("A", batch_size, mat_size); + ConstField2d B_view("B", batch_size, mat_size); + ConstField2d C_view("C", batch_size, mat_size); + ConstField2d Rhs_view("R", batch_size, mat_size); Kokkos::deep_copy(A_view, sub_diag); Kokkos::deep_copy(B_view, diag); Kokkos::deep_copy(C_view, up_diag); @@ -68,9 +68,9 @@ TEST(MatrixBatchTridiag, NotValidMatrix) int const batch_size = 2; int const mat_size = 4; - View2d A_view("A", batch_size, mat_size); - View2d B_view("B", batch_size, mat_size); - View2d C_view("C", batch_size, mat_size); + ConstField2d A_view("A", batch_size, mat_size); + ConstField2d B_view("B", batch_size, mat_size); + ConstField2d C_view("C", batch_size, mat_size); //neither positive definite symmetric nor diagonal dominant Kokkos::deep_copy(A_view, 0.93); Kokkos::deep_copy(B_view, -0.367); @@ -102,10 +102,10 @@ TEST(MatrixBatchTridiag, GeneralValidMatrix) Kokkos::View Rhs_view_host(r, batch_size, mat_size); - View2d A_view("A", batch_size, mat_size); - View2d B_view("B", batch_size, mat_size); - View2d C_view("C", batch_size, mat_size); - View2d Rhs_view("R", batch_size, mat_size); + ConstField2d A_view("A", batch_size, mat_size); + ConstField2d B_view("B", batch_size, mat_size); + ConstField2d C_view("C", batch_size, mat_size); + ConstField2d Rhs_view("R", batch_size, mat_size); Kokkos::deep_copy(A_view, A_view_host); Kokkos::deep_copy(B_view, B_view_host); Kokkos::deep_copy(C_view, C_view_host); diff --git a/vendor/sll/tests/metric_tensor.cpp b/vendor/sll/tests/metric_tensor.cpp index 4d7fd48dc..348cf0612 100644 --- a/vendor/sll/tests/metric_tensor.cpp +++ b/vendor/sll/tests/metric_tensor.cpp @@ -6,34 +6,34 @@ #include "test_utils.hpp" -struct DimX +struct X { static bool constexpr PERIODIC = false; }; -struct DimY +struct Y { static bool constexpr PERIODIC = false; }; -struct DimR +struct R { static bool constexpr PERIODIC = false; }; -struct DimP +struct P { static bool constexpr PERIODIC = true; }; -using CoordR = ddc::Coordinate; -using CoordP = ddc::Coordinate; -using CoordRP = ddc::Coordinate; +using CoordR = ddc::Coordinate; +using CoordP = ddc::Coordinate

; +using CoordRP = ddc::Coordinate; int constexpr BSDegree = 3; -struct BSplinesR : ddc::NonUniformBSplines +struct BSplinesR : ddc::NonUniformBSplines { }; -struct BSplinesP : ddc::NonUniformBSplines +struct BSplinesP : ddc::NonUniformBSplines { }; struct PolarBSplinesRP : PolarBSplines @@ -45,31 +45,31 @@ using InterpPointsR = ddc:: using InterpPointsP = ddc:: GrevilleInterpolationPoints; -struct IDimR : InterpPointsR::interpolation_discrete_dimension_type +struct GridR : InterpPointsR::interpolation_discrete_dimension_type { }; -struct IDimP : InterpPointsP::interpolation_discrete_dimension_type +struct GridP : InterpPointsP::interpolation_discrete_dimension_type { }; -using BSDomainR = ddc::DiscreteDomain; -using BSDomainP = ddc::DiscreteDomain; -using BSDomainRP = ddc::DiscreteDomain; -using BSDomainPolar = ddc::DiscreteDomain; +using BSIdxRangeR = ddc::DiscreteDomain; +using BSIdxRangeP = ddc::DiscreteDomain; +using BSIdxRangeRP = ddc::DiscreteDomain; +using BSIdxRangePolar = ddc::DiscreteDomain; -using IndexR = ddc::DiscreteElement; -using IndexP = ddc::DiscreteElement; -using IndexRP = ddc::DiscreteElement; +using IdxR = ddc::DiscreteElement; +using IdxP = ddc::DiscreteElement; +using IdxRP = ddc::DiscreteElement; -using IVectR = ddc::DiscreteVector; -using IVectP = ddc::DiscreteVector; -using IVectRP = ddc::DiscreteVector; +using IdxStepR = ddc::DiscreteVector; +using IdxStepP = ddc::DiscreteVector; +using IdxStepRP = ddc::DiscreteVector; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; template -using FieldRP = ddc::Chunk; +using FieldMemRP = ddc::Chunk; using Matrix_2x2 = std::array, 2>; @@ -102,35 +102,35 @@ class InverseMetricTensor : public testing::TestWithParam mapping; + const CircularToCartesian mapping; CoordR const r_min(0.0); CoordR const r_max(1.0); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the inverse_metric_tensor is the inverse of the metric_tensor - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 matrix; Matrix_2x2 inv_matrix; @@ -146,35 +146,35 @@ 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); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); - IndexR const r_start(1); // avoid singular point. - IndexP const p_start(0); + IdxR const r_start(1); // avoid singular point. + IdxP const p_start(0); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); - ddc::DiscreteDomain domain_r(r_start, r_size); - ddc::DiscreteDomain domain_p(p_start, p_size); - ddc::DiscreteDomain grid(domain_r, domain_p); + ddc::DiscreteDomain idx_range_r(r_start, r_size); + ddc::DiscreteDomain idx_range_p(p_start, p_size); + ddc::DiscreteDomain grid(idx_range_r, idx_range_p); - FieldRP coords(grid); - ddc::for_each(grid, [&](IndexRP const irp) { + FieldMemRP coords(grid); + ddc::for_each(grid, [&](IdxRP const irp) { coords(irp) = CoordRP( - r_min + dr * ddc::select(irp).uid(), - p_min + dp * ddc::select(irp).uid()); + r_min + dr * ddc::select(irp).uid(), + p_min + dp * ddc::select(irp).uid()); }); // Test for each coordinates if the inverse_metric_tensor is the inverse of the metric_tensor - ddc::for_each(grid, [&](IndexRP const irp) { + ddc::for_each(grid, [&](IdxRP const irp) { Matrix_2x2 matrix; Matrix_2x2 inv_matrix; diff --git a/vendor/sll/tests/polar_bsplines.cpp b/vendor/sll/tests/polar_bsplines.cpp index 18f70f0e0..1808c30b3 100644 --- a/vendor/sll/tests/polar_bsplines.cpp +++ b/vendor/sll/tests/polar_bsplines.cpp @@ -20,32 +20,29 @@ struct PolarBsplineFixture, std::integral_constant>> : public testing::Test { - struct DimR + struct R { static constexpr bool PERIODIC = false; }; - struct DimP + struct P { static constexpr bool PERIODIC = true; }; - struct DimX + struct X { static constexpr bool PERIODIC = false; }; - struct DimY + struct Y { static constexpr bool PERIODIC = false; }; static constexpr std::size_t spline_degree = D; static constexpr int continuity = C; - struct BSplineR : ddc::NonUniformBSplines + struct BSplineR : ddc::NonUniformBSplines { }; struct BSplineP - : std::conditional_t< - Uniform, - ddc::UniformBSplines, - ddc::NonUniformBSplines> + : std::conditional_t, ddc::NonUniformBSplines> { }; @@ -58,10 +55,10 @@ struct PolarBsplineFixture; - struct IDimR : GrevillePointsR::interpolation_discrete_dimension_type + struct GridR : GrevillePointsR::interpolation_discrete_dimension_type { }; - struct IDimP : GrevillePointsP::interpolation_discrete_dimension_type + struct GridP : GrevillePointsP::interpolation_discrete_dimension_type { }; struct BSplines : PolarBSplines @@ -79,49 +76,49 @@ TYPED_TEST_SUITE(PolarBsplineFixture, Cases); TYPED_TEST(PolarBsplineFixture, PartitionOfUnity) { - using DimR = typename TestFixture::DimR; - using IDimR = typename TestFixture::IDimR; - using DVectR = ddc::DiscreteVector; - using DimP = typename TestFixture::DimP; - using IDimP = typename TestFixture::IDimP; - using DVectP = ddc::DiscreteVector; - using DimX = typename TestFixture::DimX; - using DimY = typename TestFixture::DimY; - using PolarCoord = ddc::Coordinate; + using R = typename TestFixture::R; + using GridR = typename TestFixture::GridR; + using IdxStepR = ddc::DiscreteVector; + using P = typename TestFixture::P; + using GridP = typename TestFixture::GridP; + using IdxStepP = ddc::DiscreteVector; + using X = typename TestFixture::X; + using Y = typename TestFixture::Y; + using PolarCoord = ddc::Coordinate; using BSplinesR = typename TestFixture::BSplineR; using BSplinesP = typename TestFixture::BSplineP; - using CircToCart = CircularToCartesian; + using CircToCart = CircularToCartesian; using SplineRPBuilder = ddc::SplineBuilder2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using SplineRPEvaluator = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::NullExtrapolationRule, ddc::NullExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; - using DiscreteMapping = DiscreteToCartesian; + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; + using DiscreteMapping = DiscreteToCartesian; using BSplines = typename TestFixture::BSplines; - using CoordR = ddc::Coordinate; - using CoordP = ddc::Coordinate; + using CoordR = ddc::Coordinate; + using CoordP = ddc::Coordinate

; using GrevillePointsR = typename TestFixture::GrevillePointsR; using GrevillePointsP = typename TestFixture::GrevillePointsP; @@ -133,7 +130,7 @@ TYPED_TEST(PolarBsplineFixture, PartitionOfUnity) // 1. Create BSplines { - DVectR constexpr npoints(ncells + 1); + IdxStepR constexpr npoints(ncells + 1); std::vector breaks(npoints); const double dr = (rN - r0) / ncells; for (int i(0); i < npoints; ++i) { @@ -144,7 +141,7 @@ TYPED_TEST(PolarBsplineFixture, PartitionOfUnity) if constexpr (BSplinesP::is_uniform()) { ddc::init_discrete_space(p0, pN, ncells); } else { - DVectP constexpr npoints(ncells + 1); + IdxStepP constexpr npoints(ncells + 1); std::vector breaks(npoints); const double dp = (pN - p0) / ncells; for (int i(0); i < npoints; ++i) { @@ -153,19 +150,19 @@ TYPED_TEST(PolarBsplineFixture, PartitionOfUnity) ddc::init_discrete_space(breaks); } - ddc::init_discrete_space(GrevillePointsR::template get_sampling()); - ddc::init_discrete_space(GrevillePointsP::template get_sampling()); - ddc::DiscreteDomain interpolation_domain_R( - GrevillePointsR::template get_domain()); - ddc::DiscreteDomain interpolation_domain_P( - GrevillePointsP::template get_domain()); - ddc::DiscreteDomain - interpolation_domain(interpolation_domain_R, interpolation_domain_P); + ddc::init_discrete_space(GrevillePointsR::template get_sampling()); + ddc::init_discrete_space(GrevillePointsP::template get_sampling()); + ddc::DiscreteDomain interpolation_idx_range_R( + GrevillePointsR::template get_domain()); + ddc::DiscreteDomain interpolation_idx_range_P( + GrevillePointsP::template get_domain()); + ddc::DiscreteDomain + interpolation_idx_range(interpolation_idx_range_R, interpolation_idx_range_P); - SplineRPBuilder builder_rp(interpolation_domain); + SplineRPBuilder builder_rp(interpolation_idx_range); ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; + ddc::PeriodicExtrapolationRule

p_extrapolation_rule; SplineRPEvaluator evaluator_rp( r_extrapolation_rule, r_extrapolation_rule, diff --git a/vendor/sll/tests/polar_splines.cpp b/vendor/sll/tests/polar_splines.cpp index 5713d911a..fa02514ce 100644 --- a/vendor/sll/tests/polar_splines.cpp +++ b/vendor/sll/tests/polar_splines.cpp @@ -12,19 +12,19 @@ #include -struct DimR +struct R { static constexpr bool PERIODIC = false; }; -struct DimP +struct P { static constexpr bool PERIODIC = true; }; -struct DimX +struct X { static constexpr bool PERIODIC = false; }; -struct DimY +struct Y { static constexpr bool PERIODIC = false; }; @@ -33,15 +33,15 @@ static constexpr std::size_t spline_r_degree = DEGREE_R; static constexpr std::size_t spline_p_degree = DEGREE_P; static constexpr int continuity = CONTINUITY; -struct BSplinesR : ddc::NonUniformBSplines +struct BSplinesR : ddc::NonUniformBSplines { }; #if defined(BSPLINES_TYPE_UNIFORM) -struct BSplinesP : ddc::UniformBSplines +struct BSplinesP : ddc::UniformBSplines { }; #elif defined(BSPLINES_TYPE_NON_UNIFORM) -struct BSplinesP : ddc::NonUniformBSplines +struct BSplinesP : ddc::NonUniformBSplines { }; #endif @@ -51,10 +51,10 @@ using GrevillePointsR = ddc:: using GrevillePointsP = ddc:: GrevilleInterpolationPoints; -struct IDimR : GrevillePointsR::interpolation_discrete_dimension_type +struct GridR : GrevillePointsR::interpolation_discrete_dimension_type { }; -struct IDimP : GrevillePointsP::interpolation_discrete_dimension_type +struct GridP : GrevillePointsP::interpolation_discrete_dimension_type { }; struct BSplines : PolarBSplines @@ -62,16 +62,16 @@ 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) { - using PolarCoord = ddc::Coordinate; - using CoordR = ddc::Coordinate; - using CoordP = ddc::Coordinate; + using PolarCoord = ddc::Coordinate; + using CoordR = ddc::Coordinate; + using CoordP = ddc::Coordinate

; using Spline = PolarSpline; using Evaluator = PolarSplineEvaluator; using BuilderRP = ddc::SplineBuilder2D< @@ -79,30 +79,30 @@ TEST(PolarSplineTest, ConstantEval) Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using EvaluatorRP = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::NullExtrapolationRule, ddc::NullExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; - using DiscreteMapping = DiscreteToCartesian; + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; + using DiscreteMapping = DiscreteToCartesian; CoordR constexpr r0(0.); CoordR constexpr rN(1.); @@ -112,7 +112,7 @@ TEST(PolarSplineTest, ConstantEval) // 1. Create BSplines { - ddc::DiscreteVector constexpr npoints_r(ncells + 1); + ddc::DiscreteVector constexpr npoints_r(ncells + 1); std::vector breaks_r(npoints_r); const double dr = (rN - r0) / ncells; for (int i(0); i < npoints_r; ++i) { @@ -122,7 +122,7 @@ TEST(PolarSplineTest, ConstantEval) #if defined(BSPLINES_TYPE_UNIFORM) ddc::init_discrete_space(p0, pN, ncells); #elif defined(BSPLINES_TYPE_NON_UNIFORM) - ddc::DiscreteVector constexpr npoints_p(ncells + 1); + ddc::DiscreteVector constexpr npoints_p(ncells + 1); std::vector breaks_p(npoints_p); const double dp = (pN - p0) / ncells; for (int i(0); i < npoints_r; ++i) { @@ -132,17 +132,17 @@ TEST(PolarSplineTest, ConstantEval) #endif } - ddc::init_discrete_space(GrevillePointsR::get_sampling()); - ddc::init_discrete_space(GrevillePointsP::get_sampling()); - ddc::DiscreteDomain interpolation_domain_R(GrevillePointsR::get_domain()); - ddc::DiscreteDomain interpolation_domain_P(GrevillePointsP::get_domain()); - ddc::DiscreteDomain - interpolation_domain(interpolation_domain_R, interpolation_domain_P); + ddc::init_discrete_space(GrevillePointsR::get_sampling()); + ddc::init_discrete_space(GrevillePointsP::get_sampling()); + ddc::DiscreteDomain interpolation_idx_range_R(GrevillePointsR::get_domain()); + ddc::DiscreteDomain interpolation_idx_range_P(GrevillePointsP::get_domain()); + ddc::DiscreteDomain + interpolation_idx_range(interpolation_idx_range_R, interpolation_idx_range_P); - BuilderRP builder_rp(interpolation_domain); + BuilderRP builder_rp(interpolation_idx_range); ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; + ddc::PeriodicExtrapolationRule

p_extrapolation_rule; EvaluatorRP evaluator_rp( r_extrapolation_rule, r_extrapolation_rule, diff --git a/vendor/sll/tests/polynomial_evaluator.hpp b/vendor/sll/tests/polynomial_evaluator.hpp index 6c25c386c..2df6a3a1d 100644 --- a/vendor/sll/tests/polynomial_evaluator.hpp +++ b/vendor/sll/tests/polynomial_evaluator.hpp @@ -1,5 +1,4 @@ #pragma once - #include #include #include @@ -9,13 +8,22 @@ #include +#include "ddc_aliases.hpp" + +/** + * A wrapper around a polynomial Evaluator that can be used in tests. + */ struct PolynomialEvaluator { - template + /** + * @brief An analytical evaluator to be used for exact comparisons in tests. + */ + template class Evaluator { public: - using Dim = DDim; + /// The grid dimension. + using Dim = Grid1D; private: std::array m_coeffs; @@ -23,45 +31,76 @@ struct PolynomialEvaluator double const m_xN; public: - template - Evaluator(Domain domain) + /** + * @brief A constructor taking the range over which the evaluator will be applied. + * The polynomial is initialised with random values between 0.0 and 1.0 + * + * @param[in] idx_range The range of the grid over which the evaluator will be applied. + */ + template + Evaluator(IdxRange idx_range) : m_degree(Degree) - , m_xN(std::max(std::abs(rmin(domain)), std::abs(rmax(domain)))) + , m_xN(std::max(std::abs(rmin(idx_range)), std::abs(rmax(idx_range)))) { for (int i(0); i < m_degree + 1; ++i) { m_coeffs[i] = double(rand() % 100) / 100.0; } } + /** + * Evaluate the equation at x + * @param[in] x The position where the equation is evaluated. + * @returns The result of the equation. + */ double operator()(double const x) const noexcept { return eval(x, 0); } - void operator()(ddc::ChunkSpan> chunk) const + /** + * Evaluate the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + */ + void operator()(ddc::ChunkSpan> chunk) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - for (ddc::DiscreteElement const i : domain) { + for (ddc::DiscreteElement const i : idx_range) { chunk(i) = eval(ddc::coordinate(i), 0); } } + /** + * Evaluate the derivative of the equation at x + * @param[in] x The position where the equation is evaluated. + * @param[in] derivative The order of the derivative. + * @returns The result of the equation. + */ double deriv(double const x, int const derivative) const noexcept { return eval(x, derivative); } - void deriv(ddc::ChunkSpan> chunk, int const derivative) + /** + * Evaluate the derivative of the equation at the positions on which chunk is defined. + * @param[out] chunk The result of the equation. + * @param[in] derivative The order of the derivative. + */ + void deriv(ddc::ChunkSpan> chunk, int const derivative) const { - auto const& domain = chunk.domain(); + auto const& idx_range = chunk.idx_range(); - for (ddc::DiscreteElement const i : domain) { + for (ddc::DiscreteElement const i : idx_range) { chunk(i) = eval(ddc::coordinate(i), derivative); } } + /** + * The maximum norm of this equation. Used for normalisation. + * @param[in] diff The derivative whose max norm should be returned. + * @returns The maximum value of the infinite norm. + */ double max_norm(int diff = 0) const { return std::abs(deriv(m_xN, diff)); diff --git a/vendor/sll/tests/pseudo_cartesian_matrix.cpp b/vendor/sll/tests/pseudo_cartesian_matrix.cpp index 94c68789d..48925982d 100644 --- a/vendor/sll/tests/pseudo_cartesian_matrix.cpp +++ b/vendor/sll/tests/pseudo_cartesian_matrix.cpp @@ -17,32 +17,32 @@ template class PseudoCartesianJacobianMatrixTest { public: - struct DimX + struct X { }; - struct DimY + struct Y { }; - struct DimR + struct R { static bool constexpr PERIODIC = false; }; - struct DimP + struct P { static bool constexpr PERIODIC = true; }; - using CoordR = ddc::Coordinate; - using CoordP = ddc::Coordinate; - using CoordRP = ddc::Coordinate; + using CoordR = ddc::Coordinate; + using CoordP = ddc::Coordinate

; + using CoordRP = ddc::Coordinate; static int constexpr BSDegree = 3; - struct BSplinesR : ddc::NonUniformBSplines + struct BSplinesR : ddc::NonUniformBSplines { }; - struct BSplinesP : ddc::NonUniformBSplines + struct BSplinesP : ddc::NonUniformBSplines { }; @@ -57,32 +57,32 @@ class PseudoCartesianJacobianMatrixTest ddc::BoundCond::PERIODIC>; - struct IDimR : InterpPointsR::interpolation_discrete_dimension_type + struct GridR : InterpPointsR::interpolation_discrete_dimension_type { }; - struct IDimP : InterpPointsP::interpolation_discrete_dimension_type + struct GridP : InterpPointsP::interpolation_discrete_dimension_type { }; - using BSDomainR = ddc::DiscreteDomain; - using BSDomainP = ddc::DiscreteDomain; - using BSDomainRP = ddc::DiscreteDomain; + using BSIdxRangeR = ddc::DiscreteDomain; + using BSIdxRangeP = ddc::DiscreteDomain; + using BSIdxRangeRP = ddc::DiscreteDomain; - using IndexR = ddc::DiscreteElement; - using IndexP = ddc::DiscreteElement; - using IndexRP = ddc::DiscreteElement; + using IdxR = ddc::DiscreteElement; + using IdxP = ddc::DiscreteElement; + using IdxRP = ddc::DiscreteElement; - using IVectR = ddc::DiscreteVector; - using IVectP = ddc::DiscreteVector; - using IVectRP = ddc::DiscreteVector; + using IdxStepR = ddc::DiscreteVector; + using IdxStepP = ddc::DiscreteVector; + using IdxStepRP = ddc::DiscreteVector; - using IDomainR = ddc::DiscreteDomain; - using IDomainP = ddc::DiscreteDomain; - using IDomainRP = ddc::DiscreteDomain; + using IdxRangeR = ddc::DiscreteDomain; + using IdxRangeP = ddc::DiscreteDomain; + using IdxRangeRP = ddc::DiscreteDomain; using SplineRPBuilder = ddc::SplineBuilder2D< @@ -90,37 +90,37 @@ class PseudoCartesianJacobianMatrixTest Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using SplineRPEvaluator = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::NullExtrapolationRule, ddc::NullExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; - using DiscreteMapping = DiscreteToCartesian; + using DiscreteMapping = DiscreteToCartesian; - using spline_domain = ddc::DiscreteDomain; + using spline_idx_range = ddc::DiscreteDomain; template - using FieldRP = ddc::Chunk; + using FieldMemRP = ddc::Chunk; using Matrix_2x2 = std::array, 2>; @@ -139,7 +139,7 @@ class PseudoCartesianJacobianMatrixTest * mapping of a circular mapping and a discrete mapping of a Czarny * mapping. * - * Create a Nr x Nt discrete domain with Nr = N, Nt = 2*N and N a + * Create a Nr x Nt discrete index range with Nr = N, Nt = 2*N and N a * templated parameter. * Then compute the infinity norm of the difference between the * the pseudo Cartesian Jacobian matrix of the analytical mapping @@ -151,11 +151,11 @@ class PseudoCartesianJacobianMatrixTest // --- Define the grid. --------------------------------------------------------------------------- CoordR const r_min(0.0); CoordR const r_max(1.0); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); @@ -178,17 +178,17 @@ class PseudoCartesianJacobianMatrixTest ddc::init_discrete_space(r_knots); ddc::init_discrete_space(p_knots); - ddc::init_discrete_space(InterpPointsR::template get_sampling()); - ddc::init_discrete_space(InterpPointsP::template get_sampling()); + ddc::init_discrete_space(InterpPointsR::template get_sampling()); + ddc::init_discrete_space(InterpPointsP::template get_sampling()); - IDomainR interpolation_domain_R(InterpPointsR::template get_domain()); - IDomainP interpolation_domain_P(InterpPointsP::template get_domain()); - IDomainRP grid(interpolation_domain_R, interpolation_domain_P); + IdxRangeR interpolation_idx_range_R(InterpPointsR::template get_domain()); + IdxRangeP interpolation_idx_range_P(InterpPointsP::template get_domain()); + IdxRangeRP grid(interpolation_idx_range_R, interpolation_idx_range_P); // --- Define the operators. ---------------------------------------------------------------------- SplineRPBuilder const builder(grid); ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; + ddc::PeriodicExtrapolationRule

p_extrapolation_rule; SplineRPEvaluator spline_evaluator( r_extrapolation_rule, r_extrapolation_rule, @@ -202,7 +202,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; DiscreteMapping const discrete_mapping_circ = DiscreteMapping:: analytical_to_discrete(analytical_mapping_circ, builder, spline_evaluator); @@ -216,7 +216,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); DiscreteMapping const discrete_mapping_czar = DiscreteMapping:: analytical_to_discrete(analytical_mapping_czar, builder, spline_evaluator); diff --git a/vendor/sll/tests/refined_discrete_mapping.cpp b/vendor/sll/tests/refined_discrete_mapping.cpp index f13d33116..b850fb398 100644 --- a/vendor/sll/tests/refined_discrete_mapping.cpp +++ b/vendor/sll/tests/refined_discrete_mapping.cpp @@ -16,34 +16,34 @@ namespace { -struct RDimX +struct X { }; -struct RDimY +struct Y { }; -struct RDimR +struct R { static bool constexpr PERIODIC = false; }; -struct RDimP +struct P { static bool constexpr PERIODIC = true; }; -using CoordR = ddc::Coordinate; -using CoordP = ddc::Coordinate; -using CoordRP = ddc::Coordinate; +using CoordR = ddc::Coordinate; +using CoordP = ddc::Coordinate

; +using CoordRP = ddc::Coordinate; -using CoordXY = ddc::Coordinate; +using CoordXY = ddc::Coordinate; int constexpr BSDegree = 3; -struct BSplinesR : ddc::NonUniformBSplines +struct BSplinesR : ddc::NonUniformBSplines { }; -struct BSplinesP : ddc::NonUniformBSplines +struct BSplinesP : ddc::NonUniformBSplines { }; @@ -53,10 +53,10 @@ using SplineInterpPointsR = ddc:: using SplineInterpPointsP = ddc:: GrevilleInterpolationPoints; -struct IDimR : SplineInterpPointsR::interpolation_discrete_dimension_type +struct GridR : SplineInterpPointsR::interpolation_discrete_dimension_type { }; -struct IDimP : SplineInterpPointsP::interpolation_discrete_dimension_type +struct GridP : SplineInterpPointsP::interpolation_discrete_dimension_type { }; @@ -65,66 +65,66 @@ using SplineRPBuilder = ddc::SplineBuilder2D< Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, + GridR, + GridP, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE, ddc::BoundCond::PERIODIC, ddc::BoundCond::PERIODIC, ddc::SplineSolver::LAPACK, - IDimR, - IDimP>; + GridR, + GridP>; using SplineRPEvaluator = ddc::SplineEvaluator2D< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesR, BSplinesP, - IDimR, - IDimP, - ddc::ConstantExtrapolationRule, - ddc::ConstantExtrapolationRule, - ddc::PeriodicExtrapolationRule, - ddc::PeriodicExtrapolationRule, - IDimR, - IDimP>; - -using BSDomainR = ddc::DiscreteDomain; -using BSDomainP = ddc::DiscreteDomain; -using BSDomainRP = ddc::DiscreteDomain; - -using IDomainR = ddc::DiscreteDomain; -using IDomainP = ddc::DiscreteDomain; -using IDomainRP = ddc::DiscreteDomain; - -using IndexR = ddc::DiscreteElement; -using IndexP = ddc::DiscreteElement; -using IndexRP = ddc::DiscreteElement; - -using IVectR = ddc::DiscreteVector; -using IVectP = ddc::DiscreteVector; -using IVectRP = ddc::DiscreteVector; + GridR, + GridP, + ddc::ConstantExtrapolationRule, + ddc::ConstantExtrapolationRule, + ddc::PeriodicExtrapolationRule

, + ddc::PeriodicExtrapolationRule

, + GridR, + GridP>; + +using BSIdxRangeR = ddc::DiscreteDomain; +using BSIdxRangeP = ddc::DiscreteDomain; +using BSIdxRangeRP = ddc::DiscreteDomain; + +using IdxRangeR = ddc::DiscreteDomain; +using IdxRangeP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; + +using IdxR = ddc::DiscreteElement; +using IdxP = ddc::DiscreteElement; +using IdxRP = ddc::DiscreteElement; + +using IdxStepR = ddc::DiscreteVector; +using IdxStepP = ddc::DiscreteVector; +using IdxStepRP = ddc::DiscreteVector; template -using FieldR = ddc::Chunk; +using FieldMemR = ddc::Chunk; template -using FieldP = ddc::Chunk; +using FieldMemP = ddc::Chunk; template -using FieldRP = ddc::Chunk; +using FieldMemRP = ddc::Chunk; -using IDomainRP = ddc::DiscreteDomain; +using IdxRangeRP = ddc::DiscreteDomain; -using CzarnyMapping = CzarnyToCartesian; -using CircularMapping = CircularToCartesian; +using CzarnyMapping = CzarnyToCartesian; +using CircularMapping = CircularToCartesian; template -using FieldRP = ddc::Chunk; +using FieldMemRP = ddc::Chunk; using Matrix_2x2 = std::array, 2>; @@ -138,29 +138,27 @@ using Matrix_2x2 = std::array, 2>; * The mapping we are testing. * @param[in] analytical_mappping * The mapping analytically defined. - * @param[in] domain - * The domain on which we test the values. + * @param[in] idx_range + * The index range on which we test the values. */ template double check_value_on_grid( DiscreteMapping const& mapping, Mapping const& analytical_mapping, - IDomainRP const& domain) + IdxRangeRP const& idx_range) { const double TOL = 1e-7; double max_err = 0.; - ddc::for_each(domain, [&](IndexRP const irp) { + ddc::for_each(idx_range, [&](IdxRP const irp) { const CoordRP coord(ddc::coordinate(irp)); const CoordXY discrete_coord = mapping(coord); const CoordXY analytical_coord = analytical_mapping(coord); - EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); - EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); + EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); + EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); - const double diff_x - = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); - const double diff_y - = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); + const double diff_x = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); + const double diff_y = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); max_err = max_err > diff_x ? max_err : diff_x; max_err = max_err > diff_y ? max_err : diff_y; @@ -181,22 +179,22 @@ double check_value_on_grid( * The mapping we are testing. * @param[in] analytical_mappping * The mapping analytically defined. - * @param[in] domain - * The domain on which we test the values. + * @param[in] idx_range + * The index range on which we test the values. */ template double check_value_not_on_grid( DiscreteMapping const& mapping, Mapping const& analytical_mapping, - IDomainRP const& domain) + IdxRangeRP const& idx_range) { std::srand(100); - FieldRP coords(domain); - IndexR ir_max(ddc::select(domain).back()); - IndexP ip_max(ddc::select(domain).back()); - ddc::for_each(domain, [&](IndexRP const irp) { - IndexR ir(ddc::select(irp)); + FieldMemRP coords(idx_range); + IdxR ir_max(ddc::select(idx_range).back()); + IdxP ip_max(ddc::select(idx_range).back()); + ddc::for_each(idx_range, [&](IdxRP const irp) { + IdxR ir(ddc::select(irp)); CoordR coordr_0 = ddc::coordinate(ir); CoordR coordr_1 = ddc::coordinate(ir + 1); double coord_r; @@ -207,7 +205,7 @@ double check_value_not_on_grid( coord_r = coordr_0; } - IndexP ip(ddc::select(irp)); + IdxP ip(ddc::select(irp)); CoordP coordp_0 = ddc::coordinate(ip); CoordP coordp_1 = ddc::coordinate(ip + 1); double coord_p; @@ -222,18 +220,16 @@ double check_value_not_on_grid( const double TOL = 5e-5; double max_err = 0.; - ddc::for_each(domain, [&](IndexRP const irp) { + ddc::for_each(idx_range, [&](IdxRP const irp) { const CoordRP coord(coords(irp)); const CoordXY discrete_coord = mapping(coord); const CoordXY analytical_coord = analytical_mapping(coord); - EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); - EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); + EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); + EXPECT_NEAR(ddc::get(discrete_coord), ddc::get(analytical_coord), TOL); - const double diff_x - = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); - const double diff_y - = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); + const double diff_x = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); + const double diff_y = double(ddc::get(discrete_coord) - ddc::get(analytical_coord)); max_err = max_err > diff_x ? max_err : diff_x; max_err = max_err > diff_y ? max_err : diff_y; @@ -249,15 +245,14 @@ double test_on_grid_and_not_on_grid( int const Nt, Mapping const& analytical_mapping, RefinedDiscreteToCartesian< - RDimX, - RDimY, + X, + Y, SplineRPBuilder, SplineRPEvaluator, refined_Nr, refined_Nt> const& refined_mapping, - DiscreteToCartesian const& - discrete_mapping, - IDomainRP const& grid) + DiscreteToCartesian const& discrete_mapping, + IdxRangeRP const& grid) { std::cout << std::endl << "DOMAIN Nr x Nt = " << Nr << " x " << Nt @@ -290,17 +285,17 @@ double test_on_grid_and_not_on_grid( * The mapping we are testing. * @param[in] analytical_mappping * The mapping analytically defined. - * @param[in] domain - * The domain on which we test the values. + * @param[in] idx_range + * The index range on which we test the values. */ template double check_Jacobian_on_grid( DiscreteMapping const& mapping, Mapping const& analytical_mapping, - IDomainRP const& domain) + IdxRangeRP const& idx_range) { double max_err = 0.; - ddc::for_each(domain, [&](IndexRP const irp) { + ddc::for_each(idx_range, [&](IdxRP const irp) { const CoordRP coord(ddc::coordinate(irp)); Matrix_2x2 discrete_Jacobian; @@ -336,22 +331,22 @@ double check_Jacobian_on_grid( * The mapping we are testing. * @param[in] analytical_mappping * The mapping analytically defined. - * @param[in] domain - * The domain on which we test the values. + * @param[in] idx_range + * The index range on which we test the values. */ template double check_Jacobian_not_on_grid( DiscreteMapping const& mapping, Mapping const& analytical_mapping, - IDomainRP const& domain) + IdxRangeRP const& idx_range) { std::srand(100); - FieldRP coords(domain); - IndexR ir_max(ddc::select(domain).back()); - IndexP ip_max(ddc::select(domain).back()); - ddc::for_each(domain, [&](IndexRP const irp) { - IndexR ir(ddc::select(irp)); + FieldMemRP coords(idx_range); + IdxR ir_max(ddc::select(idx_range).back()); + IdxP ip_max(ddc::select(idx_range).back()); + ddc::for_each(idx_range, [&](IdxRP const irp) { + IdxR ir(ddc::select(irp)); CoordR coordr_0 = ddc::coordinate(ir); CoordR coordr_1 = ddc::coordinate(ir + 1); double coord_r; @@ -362,7 +357,7 @@ double check_Jacobian_not_on_grid( coord_r = coordr_0; } - IndexP ip(ddc::select(irp)); + IdxP ip(ddc::select(irp)); CoordP coordp_0 = ddc::coordinate(ip); CoordP coordp_1 = ddc::coordinate(ip + 1); double coord_p; @@ -377,7 +372,7 @@ double check_Jacobian_not_on_grid( double max_err = 0.; - ddc::for_each(domain, [&](IndexRP const irp) { + ddc::for_each(idx_range, [&](IdxRP const irp) { const CoordRP coord(coords(irp)); Matrix_2x2 discrete_Jacobian; @@ -407,15 +402,14 @@ double test_Jacobian( int const Nt, Mapping const& analytical_mapping, RefinedDiscreteToCartesian< - RDimX, - RDimY, + X, + Y, SplineRPBuilder, SplineRPEvaluator, refined_Nr, refined_Nt> const& refined_mapping, - DiscreteToCartesian const& - discrete_mapping, - IDomainRP const& grid) + DiscreteToCartesian const& discrete_mapping, + IdxRangeRP const& grid) { std::cout << std::endl << "DOMAIN Nr x Nt = " << Nr << " x " << Nt @@ -449,20 +443,21 @@ double test_Jacobian( * The mapping we are testing. * @param[in] analytical_mappping * The mapping analytically defined. - * @param[in] domain - * The domain on which we test the values. + * @param[in] idx_range + * The index range on which we test the values. */ template double check_pseudo_Cart( DiscreteMapping const& mapping, Mapping const& analytical_mapping, - IDomainRP const& domain) + IdxRangeRP const& idx_range) { double max_err = 0.; Matrix_2x2 discrete_pseudo_Cart; Matrix_2x2 analytical_pseudo_Cart; - mapping.to_pseudo_cartesian_jacobian_center_matrix(domain, discrete_pseudo_Cart); - analytical_mapping.to_pseudo_cartesian_jacobian_center_matrix(domain, analytical_pseudo_Cart); + mapping.to_pseudo_cartesian_jacobian_center_matrix(idx_range, discrete_pseudo_Cart); + analytical_mapping + .to_pseudo_cartesian_jacobian_center_matrix(idx_range, analytical_pseudo_Cart); const double diff_11 = double(discrete_pseudo_Cart[0][0] - analytical_pseudo_Cart[0][0]); const double diff_12 = double(discrete_pseudo_Cart[0][1] - analytical_pseudo_Cart[0][1]); @@ -485,15 +480,14 @@ double test_pseudo_Cart( int const Nt, Mapping const& analytical_mapping, RefinedDiscreteToCartesian< - RDimX, - RDimY, + X, + Y, SplineRPBuilder, SplineRPEvaluator, refined_Nr, refined_Nt> const& refined_mapping, - DiscreteToCartesian const& - discrete_mapping, - IDomainRP const& grid) + DiscreteToCartesian const& discrete_mapping, + IdxRangeRP const& grid) { std::cout << std::endl << "DOMAIN Nr x Nt = " << Nr << " x " << Nt @@ -520,17 +514,17 @@ TEST(RefinedDiscreteMapping, TestRefinedDiscreteMapping) const CzarnyMapping analytical_mapping(0.3, 1.4); //const CircularMapping analytical_mapping; - // Discrete domain --- + // Discrete index range --- int const Nr = 16; int const Nt = 32; CoordR const r_min(0.0); CoordR const r_max(1.0); - IVectR const r_size(Nr); + IdxStepR const r_size(Nr); CoordP const p_min(0.0); CoordP const p_max(2.0 * M_PI); - IVectP const p_size(Nt); + IdxStepP const p_size(Nt); double const dr((r_max - r_min) / r_size); double const dp((p_max - p_min) / p_size); @@ -552,48 +546,44 @@ TEST(RefinedDiscreteMapping, TestRefinedDiscreteMapping) ddc::init_discrete_space(r_knots); ddc::init_discrete_space(p_knots); - ddc::init_discrete_space(SplineInterpPointsR::get_sampling()); - ddc::init_discrete_space(SplineInterpPointsP::get_sampling()); + ddc::init_discrete_space(SplineInterpPointsR::get_sampling()); + ddc::init_discrete_space(SplineInterpPointsP::get_sampling()); - IDomainR interpolation_domain_R(SplineInterpPointsR::get_domain()); - IDomainP interpolation_domain_P(SplineInterpPointsP::get_domain()); - IDomainRP grid(interpolation_domain_R, interpolation_domain_P); + IdxRangeR interpolation_idx_range_R(SplineInterpPointsR::get_domain()); + IdxRangeP interpolation_idx_range_P(SplineInterpPointsP::get_domain()); + IdxRangeRP grid(interpolation_idx_range_R, interpolation_idx_range_P); // Operators --- SplineRPBuilder builder(grid); - ddc::ConstantExtrapolationRule boundary_condition_r_left(r_min); - ddc::ConstantExtrapolationRule boundary_condition_r_right(r_max); + ddc::ConstantExtrapolationRule boundary_condition_r_left(r_min); + ddc::ConstantExtrapolationRule boundary_condition_r_right(r_max); SplineRPEvaluator spline_evaluator( boundary_condition_r_left, boundary_condition_r_right, - ddc::PeriodicExtrapolationRule(), - ddc::PeriodicExtrapolationRule()); + ddc::PeriodicExtrapolationRule

(), + ddc::PeriodicExtrapolationRule

()); // Tests --- std::array results; DiscreteToCartesian discrete_mapping - = DiscreteToCartesian:: + = DiscreteToCartesian:: analytical_to_discrete(analytical_mapping, builder, spline_evaluator); RefinedDiscreteToCartesian refined_mapping_16x32 - = RefinedDiscreteToCartesian:: + = RefinedDiscreteToCartesian:: analytical_to_refined(analytical_mapping, grid); RefinedDiscreteToCartesian refined_mapping_32x64 - = RefinedDiscreteToCartesian:: + = RefinedDiscreteToCartesian:: analytical_to_refined(analytical_mapping, grid); - RefinedDiscreteToCartesian refined_mapping_64x128 = RefinedDiscreteToCartesian< - RDimX, - RDimY, - SplineRPBuilder, - SplineRPEvaluator, - 64, - 128>::analytical_to_refined(analytical_mapping, grid); + RefinedDiscreteToCartesian refined_mapping_64x128 + = RefinedDiscreteToCartesian:: + analytical_to_refined(analytical_mapping, grid); std::cout << std::endl