diff --git a/docs/Adding_docs.md b/docs/Adding_docs.md index 3509b0cc9..5febece7f 100644 --- a/docs/Adding_docs.md +++ b/docs/Adding_docs.md @@ -107,7 +107,7 @@ The `README.md` file should begin with a title: ```markdown # ``` -If the folder doesn't contain a `.private` file (indicating that the folder is not yet ready to be shared publically), the new page can be referenced from any enclosing page that may exist. For example if you create a new folder in `src/`, you should add the following line to the list of folders in `src/README.md`: +If the folder doesn't contain a `.private` file (indicating that the folder is not yet ready to be shared publicly), the new page can be referenced from any enclosing page that may exist. For example if you create a new folder in `src/`, you should add the following line to the list of folders in `src/README.md`: ```markdown - [my_new_folder](./my_new_folder/README.md) : Short description of contents. ``` @@ -150,4 +150,4 @@ $$ ``` The first syntax can be used when the expression fits in one line and doesn't use markdown special characters. The second syntax can be used in all contexts. -In C++ header files it should be blocked with Doxagen syntax, i.e. `@f$` instead of `$`, and `@f[` and `@f]` instead of `$$`. +In C++ header files it should be blocked with Doxygen syntax, i.e. `@f$` instead of `$`, and `@f[` and `@f]` instead of `$$`. diff --git a/simulations/geometryVparMu/collisions/CMakeLists.txt b/simulations/geometryVparMu/collisions/CMakeLists.txt index 971267822..fbd7a480d 100644 --- a/simulations/geometryVparMu/collisions/CMakeLists.txt +++ b/simulations/geometryVparMu/collisions/CMakeLists.txt @@ -8,6 +8,7 @@ target_link_libraries(test_collSpVparMu PDI::pdi gslx::collisions gslx::geometry_vparmu + gslx::geometry_collisions_vparmu gslx::initialization_vparmu gslx::io gslx::quadrature diff --git a/simulations/geometryVparMu/collisions/coll_1specie_maxwellian.yaml b/simulations/geometryVparMu/collisions/coll_1specie_maxwellian.yaml index 7bac316ed..3a98d23e9 100644 --- a/simulations/geometryVparMu/collisions/coll_1specie_maxwellian.yaml +++ b/simulations/geometryVparMu/collisions/coll_1specie_maxwellian.yaml @@ -13,9 +13,9 @@ SpeciesInfo: temperature_eq: 1. mean_velocity_eq: 0. -CollisionsInfo: +Collisions: nustar0_rpeak: 1. - collisions_interspecies: false + interspecies: false Algorithm: deltat: 0.01 # Attention: must be < 1.e-2 because plays the role of nustar here diff --git a/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentT.yaml b/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentT.yaml index 0dea8e9d0..131961aac 100644 --- a/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentT.yaml +++ b/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentT.yaml @@ -18,9 +18,9 @@ SpeciesInfo: temperature_eq: 1.1 mean_velocity_eq: 0. -CollisionsInfo: +Collisions: nustar0_rpeak: 1. - collisions_interspecies: true + interspecies: true Algorithm: deltat: 0.0001 # Attention: must be < 1.e-2 because plays the role of nustar here diff --git a/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentV.yaml b/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentV.yaml index 890c5f429..37f72fb42 100644 --- a/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentV.yaml +++ b/simulations/geometryVparMu/collisions/coll_2species_maxwellian_differentV.yaml @@ -18,9 +18,9 @@ SpeciesInfo: temperature_eq: 1. mean_velocity_eq: 0.1 -CollisionsInfo: +Collisions: nustar0_rpeak: 1. - collisions_interspecies: true + interspecies: true Algorithm: deltat: 0.0001 # Attention: must be < 1.e-2 because plays the role of nustar here diff --git a/simulations/geometryVparMu/collisions/params.yaml.hpp b/simulations/geometryVparMu/collisions/params.yaml.hpp index 995c0a0d9..77c4f4d9c 100644 --- a/simulations/geometryVparMu/collisions/params.yaml.hpp +++ b/simulations/geometryVparMu/collisions/params.yaml.hpp @@ -22,9 +22,9 @@ constexpr char const* const params_yaml = R"PDI_CFG(SplineMesh: temperature_eq: 1.1 mean_velocity_eq: 0. -CollisionsInfo: +Collisions: nustar0_rpeak: 1. - collisions_interspecies: true + interspecies: true Algorithm: deltat: 0.01 diff --git a/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp b/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp index 98a5b2685..914caf4a8 100644 --- a/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp +++ b/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp @@ -11,8 +11,8 @@ #include <paraconf.h> #include <pdi.h> -#include "CollisionSpVparMu.hpp" -#include "collisioninfo.hpp" +#include "collision_configuration.hpp" +#include "collision_operator.hpp" #include "ddc_alias_inline_functions.hpp" #include "geometry.hpp" #include "input.hpp" @@ -67,7 +67,7 @@ int main(int argc, char** argv) = MaxwellianEquilibrium::init_from_input(idxrange_kinsp, conf_collision); init_fequilibrium(get_field(allfequilibrium)); - // ---> Initialisation of the distribution function as a pertubed Maxwellian + // ---> Initialisation of the distribution function as a perturbed Maxwellian DFieldMemSpVparMu allfdistribu(idxrange_spvparmu); NoPerturbInitialization const init(get_const_field(allfequilibrium)); init(get_field(allfdistribu)); @@ -98,14 +98,14 @@ int main(int argc, char** argv) SplineMuBuilder const builder_mu(idxrange_mu); DFieldMemMu const coeff_intdmu(simpson_trapezoid_quadrature_coefficients_1d< Kokkos::DefaultExecutionSpace>(idxrange_mu, Extremity::BACK)); - CollisionInfo const collision_info(conf_collision); - CollisionSpVparMu<CollisionInfo, IdxRangeSpVparMu, GridVpar, GridMu, double> collision_operator( - collision_info, + CollisionConfiguration const collision_configuration( + conf_collision, idxrange_spvparmu, get_const_field(coeff_intdmu), get_const_field(coeff_intdvpar), B_norm, get_const_field(Bstar_s)); + CollisionOperator collision_operator(collision_configuration); // --------- TIME ITERATION --------- // ---> Reading of algorithm info from input YAML file diff --git a/src/collisions/CMakeLists.txt b/src/collisions/CMakeLists.txt index cef4c3b0d..f40989430 100644 --- a/src/collisions/CMakeLists.txt +++ b/src/collisions/CMakeLists.txt @@ -1,7 +1,7 @@ # SPDX-License-Identifier: MIT add_library("collisions" STATIC - koliop_interface.cpp + collision_operator.cpp ) target_include_directories("collisions" diff --git a/src/collisions/CollisionSpVparMu.hpp b/src/collisions/CollisionSpVparMu.hpp deleted file mode 100644 index 17e9c22b5..000000000 --- a/src/collisions/CollisionSpVparMu.hpp +++ /dev/null @@ -1,355 +0,0 @@ -// SPDX-License-Identifier: MIT -#pragma once -#include <ddc/ddc.hpp> - -#include "assert.hpp" -#include "collisions_dimensions.hpp" -#include "ddc_alias_inline_functions.hpp" -#include "ddc_aliases.hpp" -#include "ddc_helper.hpp" -#include "koliop_interface.hpp" -#include "species_info.hpp" - -/** - * @brief A class which computes the collision operator in (vpar,mu) - */ -template < - class CollInfo, - class IdxRangeFDistrib, - class GridVpar, - class GridMu, - class InputDFieldThetaR> -class CollisionSpVparMu /* : public IRightHandSide */ -{ -private: - using InputDFieldR = typename CollInfo::radial_chunk_type; - using fdistrib_grids = ddc::to_type_seq_t<IdxRangeFDistrib>; - using GridR = typename collisions_dimensions::ExtractRDim<InputDFieldR>::type; - using GridTheta = - typename collisions_dimensions::ExtractThetaDim<InputDFieldThetaR, GridR>::type; - using InputThetaRTags = - typename collisions_dimensions::ExtractRThetaTags<InputDFieldThetaR>::type; - using InputSpThetaRVparTags = ddc::type_seq_merge_t< - ddc::detail::TypeSeq<Species>, - ddc::type_seq_merge_t<InputThetaRTags, ddc::detail::TypeSeq<GridVpar>>>; - using InputIdxRangeSpThetaRVpar - = ddc::detail::convert_type_seq_to_discrete_domain_t<InputSpThetaRVparTags>; - - // Validate template types - static_assert(ddc::is_discrete_domain_v<IdxRangeFDistrib>); - static_assert(IdxRangeFDistrib::rank() >= 3 && IdxRangeFDistrib::rank() <= 6); - static_assert((std::is_same_v<InputDFieldR, double>) || ddc::is_borrowed_chunk_v<InputDFieldR>); - static_assert( - (std::is_same_v<InputDFieldThetaR, double>) - || (ddc::is_borrowed_chunk_v<InputDFieldThetaR>)); - // Ensure expected types appear in distribution index range - static_assert( - ddc::in_tags_v<Species, fdistrib_grids>, - "The distribution function must be defined over the species"); - static_assert( - ddc::in_tags_v<GridVpar, fdistrib_grids>, - "The distribution function must be defined over the Vpar direction"); - static_assert( - ddc::in_tags_v<GridMu, fdistrib_grids>, - "The distribution function must be defined over the Mu direction"); - // Ensure that uniform - static_assert( - ddc::is_uniform_point_sampling_v<GridVpar>, - "The grid must be uniform in the vpar direction."); - static_assert( - ddc::is_uniform_point_sampling_v<GridMu>, - "The grid must be uniform in the mu direction."); - - static_assert( - ddc::type_seq_rank_v<Species, fdistrib_grids> == 0, - "Species should be the first index in the distribution function index set"); - static_assert( - collisions_dimensions:: - order_of_last_grids<fdistrib_grids, GridTheta, GridR, GridVpar, GridMu>(), - "Misordered inex set for the distribution function. Koliop expects (Sp, Phi, Theta, R, " - "Vpar, Mu)"); - -public: - /// Type alias for a field on a grid of (species, theta, r, vpar) or a subset - using InputDFieldSpThetaRVpar = DConstField<InputIdxRangeSpThetaRVpar>; - - /// Type alias for the index range of the radial points - using IdxRangeR = IdxRange<GridR>; - /// Type alias for the index range of the poloidal plane - using IdxRangeThetaR = IdxRange<GridTheta, GridR>; - /// Type alias for the index range containing (species, theta, r, vpar) - using IdxRangeSpThetaRVpar = IdxRange<Species, GridTheta, GridR, GridVpar>; - /// Type alias for the index range of the magnetic moment. - using IdxRangeMu = IdxRange<GridMu>; - /// Type alias for the index range of the velocity parallel to the magnetic field. - using IdxRangeVpar = IdxRange<GridVpar>; - - /// Type alias for the index of the poloidal plane - using IdxThetaR = Idx<GridTheta, GridR>; - /// Type alias for the index containing (species, theta, r, vpar) - using IdxSpThetaRVpar = Idx<Species, GridTheta, GridR, GridVpar>; - - /// Type alias for a field memory block on a grid of radial values - using DFieldMemR = DFieldMem<IdxRangeR>; - /// Type alias for a field memory block on a grid of magnetic moments - using DFieldMemMu = DFieldMem<IdxRangeMu>; - /// Type alias for a field memory block on a grid of parallel velocities - using DFieldMemVpar = DFieldMem<IdxRangeVpar>; - /// Type alias for a field memory block on a grid on a poloidal plane - using DFieldMemThetaR = DFieldMem<IdxRangeThetaR>; - /// Type alias for a field memory block on a grid of species, poloidal plane and parallel velocities - using DFieldMemSpThetaRVpar = DFieldMem<IdxRangeSpThetaRVpar>; - /// Type alias for a field defined on a grid of radial values - using DFieldR = DField<IdxRangeR>; - /// Type alias for a field defined on a grid on a poloidal plane - using DFieldThetaR = DField<IdxRangeThetaR>; - /// Type alias for a field on a grid of species, poloidal plane and parallel velocities - using DFieldSpThetaRVpar = DField<IdxRangeSpThetaRVpar>; - /// Type alias for a constant field on GPU defined on a grid of magnetic moments. - using DConstFieldMu = DConstField< - IdxRangeMu, - Kokkos::DefaultExecutionSpace:: - memory_space>; // Equivalent to Field<double const, IdxRangeMu> - /// Type alias for a constant field on GPU defined on a grid of parallel velocities. - using DConstFieldVpar = DConstField< - IdxRangeVpar, - Kokkos::DefaultExecutionSpace:: - memory_space>; // Equivalent to Field<double const, IdxRangeMu> - - /// Type alias for the distribution function stored on GPU. - using FDistribField = DField<IdxRangeFDistrib>; - -private: - /** - * @brief Copy the information in src (received as an input to the class) into a 1D radial profile. - * - * @param[in] src The source from which the data is copied. - * @param[out] dst The destination into which the data is copied. - */ - void deepcopy_radial_profile(DFieldR dst, InputDFieldR src) - { - if constexpr (collisions_dimensions::is_spoofed_dim_v<GridR>) { - // If InputDFieldR is a double because R is not provided - ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), dst, src); - } else { - // If InputDFieldR and DFieldR are the same type - ddc::parallel_deepcopy(dst, src); - } - } - -public: - /** - * @brief Copy the information in src (received as an input to the class) into a 2D field defined - * over (r, theta). - * This function should be private but is public due to Kokkos restrictions. - * - * @param[in] src The source from which the data is copied. - * @param[out] dst The destination into which the data is copied. - */ - void deepcopy_poloidal_plane(DFieldThetaR dst, InputDFieldThetaR src) - { - if constexpr ((collisions_dimensions::is_spoofed_dim_v<GridR>)&&( - collisions_dimensions::is_spoofed_dim_v<GridTheta>)) { - // If InputDFieldThetaR is a double because R and Theta are not provided - ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), dst, src); - } else if constexpr ((!collisions_dimensions::is_spoofed_dim_v<GridR>)&&( - !collisions_dimensions::is_spoofed_dim_v<GridTheta>)) { - // If InputDFieldThetaR and DFieldThetaR are the same type - ddc::parallel_deepcopy(dst, src); - } else { - using NonSpoofDim = std:: - conditional_t<collisions_dimensions::is_spoofed_dim_v<GridR>, GridTheta, GridR>; - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - get_idx_range(dst), - KOKKOS_LAMBDA(Idx<GridTheta, GridR> idx) { - dst(idx) = src(ddc::select<NonSpoofDim>(idx)); - }); - } - } - - /** - * @brief Copy the information in src (received as an input to the class) into a 4D field defined - * over (species, theta, r, vpar). - * This function should be private but is public due to Kokkos restrictions. - * - * @param[in] src The source from which the data is copied. - * @param[out] dst The destination into which the data is copied. - */ - void deepcopy_Bstar(DFieldSpThetaRVpar dst, InputDFieldSpThetaRVpar src) - { - if constexpr (std::is_same_v<DFieldSpThetaRVpar, InputDFieldSpThetaRVpar>) { - ddc::parallel_deepcopy(dst, src); - } else { - using InputIdxSpThetaRVpar = typename InputIdxRangeSpThetaRVpar::discrete_element_type; - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - get_idx_range(dst), - KOKKOS_LAMBDA(IdxSpThetaRVpar idx) { - InputIdxSpThetaRVpar src_idx(idx); - dst(idx) = src(src_idx); - }); - } - } - -public: - /** - * @brief Create instance of CollisionSpVparMu class - * - * @param[in] collision_info - * class containing rg, safety_factor, nustar0 and coeff_AD - * @param[in] idxrange_fdistrib - * index range (species, 3D space, 2D velocity space) on which the collision operator acts - * ATTENTION it must contain the all index range in species and 2D velocity. - * @param[in] coeff_intdmu - * quadrature coefficient in mu direction - * @param[in] coeff_intdvpar - * quadrature coefficient in vpar direction - * @param[in] B_norm - * magnetic field norm in (r,theta) - * @param[in] Bstar_s - * Bstar value for each species in (r,theta,vpar) - */ - CollisionSpVparMu( - CollInfo const& collision_info, - IdxRangeFDistrib idxrange_fdistrib, - DConstFieldMu coeff_intdmu, - DConstFieldVpar coeff_intdvpar, - InputDFieldThetaR B_norm, - InputDFieldSpThetaRVpar Bstar_s) - : m_operator_handle {} - , m_hat_As {"m_hat_As", collisions_dimensions::get_expanded_idx_range<Species>(idxrange_fdistrib)} - , m_hat_Zs {"m_hat_Zs", collisions_dimensions::get_expanded_idx_range<Species>(idxrange_fdistrib)} - , m_mask_buffer_r {"m_mask_buffer_r", collisions_dimensions::get_expanded_idx_range<GridR>(idxrange_fdistrib)} - , m_mask_LIM {"m_mask_LIM", IdxRangeThetaR {collisions_dimensions::get_expanded_idx_range<GridTheta, GridR>(idxrange_fdistrib)}} - , m_B_norm {"m_B_norm", IdxRangeThetaR {collisions_dimensions::get_expanded_idx_range<GridTheta, GridR>(idxrange_fdistrib)}} - , m_Bstar_s {"m_Bstar_s", IdxRangeSpThetaRVpar {collisions_dimensions::get_expanded_idx_range<Species, GridTheta, GridR, GridVpar>(idxrange_fdistrib)}} - , m_coeff_AD {"m_coeff_AD", collisions_dimensions::get_expanded_idx_range<GridR>(idxrange_fdistrib)} - , m_mug {"m_mug", ddc::select<GridMu>(idxrange_fdistrib)} - , m_vparg {"m_vparg", ddc::select<GridVpar>(idxrange_fdistrib)} - { - IdxRangeVpar idxrange_vpar(idxrange_fdistrib); - if (idxrange_vpar.size() % 2 != 0) { - throw std::runtime_error("The number of points in the vpar direction must be a " - "multiple of 2. This ensures that there is no grid point at " - "vpar=0 (this would cause division by 0)."); - } - IdxRangeSp idxrange_sp(idxrange_fdistrib); - // --> Initialize the mass species - host_t<DConstFieldSp> hat_As_host - = ddc::host_discrete_space<Species>().masses()[idxrange_sp]; - ddc::parallel_deepcopy(get_field(m_hat_As), hat_As_host); - // --> Initialize the charge species - host_t<DConstFieldSp> hat_Zs_host - = ddc::host_discrete_space<Species>().charges()[idxrange_sp]; - ddc::parallel_deepcopy(get_field(m_hat_Zs), hat_Zs_host); - - // --> Initialize the other quantities needed in koliop - deepcopy_radial_profile(get_field(m_coeff_AD), collision_info.coeff_AD()); - deepcopy_poloidal_plane(get_field(m_B_norm), B_norm); - deepcopy_Bstar(get_field(m_Bstar_s), Bstar_s); - - // --> Initialization of the masks that have no sense here to 0. - ddc::parallel_fill(get_field(m_mask_buffer_r), 0.0); // Masked if >= 0.99 - ddc::parallel_fill(get_field(m_mask_LIM), 0.0); // Masked if >= 0.99 - - // --> Initialization of vpar and mu grids - ddcHelper::dump_coordinates(Kokkos::DefaultExecutionSpace(), get_field(m_mug)); - ddcHelper::dump_coordinates(Kokkos::DefaultExecutionSpace(), get_field(m_vparg)); - - // NOTE: We need to fence because DumpCoordinates is asynchronous on GPUs. - Kokkos::fence(); - - std::size_t const n_mu = ddc::select<GridMu>(idxrange_fdistrib).size(); - std::size_t const n_vpar = idxrange_vpar.size(); - std::size_t const n_r = get_idx_range<GridR>(m_mask_buffer_r).size(); - std::size_t const n_theta = get_idx_range<GridTheta>(m_B_norm).size(); - std::size_t const n_sp = idxrange_sp.size(); - std::size_t const n_batch - = idxrange_fdistrib.size() / (n_mu * n_vpar * n_r * n_theta * n_sp); - - m_operator_handle = koliop_interface::DoOperatorInitialization( - n_mu, - n_vpar, - n_r, - n_theta, - n_batch, - n_sp, - collision_info.collisions_interspecies(), - /* the_local_index range_r_offset */ 0 + n_r - 1, - m_hat_As.data_handle(), - m_hat_Zs.data_handle(), - m_mug.data_handle(), - m_vparg.data_handle(), - coeff_intdmu.data_handle(), - coeff_intdvpar.data_handle(), - m_coeff_AD.data_handle(), - m_mask_buffer_r.data_handle(), - m_mask_LIM.data_handle(), - m_B_norm.data_handle(), - m_Bstar_s.data_handle()); - } - - ~CollisionSpVparMu() - { - koliop_interface::DoOperatorDeinitialization( - static_cast<::koliop_Operator>(m_operator_handle)); - } - - /** - * @brief Apply the collision operator to the distribution functions of all species on all species - * - * @param[inout] all_f_distribution - * All distribution functions - * @param[in] deltat_coll - * Collision time step - */ - void operator()(FDistribField all_f_distribution, double deltat_coll) const - { - if (::koliop_Collision( - static_cast<::koliop_Operator>(m_operator_handle), - deltat_coll, - all_f_distribution.data_handle()) - != KOLIOP_STATUS_SUCCESS) { - GSLX_ASSERT(false); - } - - // NOTE: While gslx is not fully stream compatible, just fence at the end of - // the computation. - if (::koliop_Fence(static_cast<::koliop_Operator>(m_operator_handle)) - != KOLIOP_STATUS_SUCCESS) { - GSLX_ASSERT(false); - } - } - -protected: - /** - * Opaque type representing the operator (due to the C interface) - */ - ::koliop_Operator m_operator_handle; - // NOTE: Some of these arrays should come from a parent class or manager. - // They resides in this class while we wait for their implementation. - /// Combinatory (6x6) matrix computed only one times at initialisation. Rk: 6 = 2*(Npolmax-1) + 1 + 1 - koliop_interface::MDL<double[6][6]> m_comb_mat; - /// Normalized masses for all species - DFieldMemSp m_hat_As; - /// Normalized charges for all species - DFieldMemSp m_hat_Zs; - /// Radial AD coefficients - DFieldMemR m_coeff_AD; - /// Mask used to avoid to apply collision in certain region - // [TODO]: This mask should maybe be deleted in C++ version - DFieldMemR m_mask_buffer_r; - /// Limiter mask in (r,theta) - DFieldMemThetaR m_mask_LIM; - /// B norm in (r,theta) - // [TODO] Attention this must be 3D for generalization to 3D geometry--> transfer it in a 1D array ? - DFieldMemThetaR m_B_norm; - /// Bstar(species,r,theta,vpar) - // [TODO] Must be 5D for full 3D geometry - DFieldMemSpThetaRVpar m_Bstar_s; - /// grid in mu direction - DFieldMemMu m_mug; - /// grid in vpar direction - DFieldMemVpar m_vparg; -}; diff --git a/src/collisions/README.md b/src/collisions/README.md index af44692e7..b112c64dc 100644 --- a/src/collisions/README.md +++ b/src/collisions/README.md @@ -1,5 +1,9 @@ -# Collisions +# Collisions -Collision operator in (vpar,mu) applied to all species at the same time. +Collision operator in (vpar,mu) applied to all against all species at the same time. As such, the operator is O(N^2) in number of species. -Rk: translated from Fortran version in C++ in koliop submodule +It was rewritten from Fortran into C++. The operator's name is KOkkos coLIsion OPerator (KOLIOP) and is provided as a submodule (see `/vendor`). + +While the operator is in C++, the interface with Fortran and other C++ code is made through a C layer. This is also known as the hourglass pattern (conceptually, the narrow neck of the device represent the lowering to C). + +To integrate Koliop into gyselalibxx, we wrap its functionalities into a DDC aware operator present in `collision_operator.hpp`. In gyselalibxx, operator are expected to support multiple if not all kind of geometries. But Koliop expect some data in layout right [sp, phi, theta, r, vpar, mu] instead of the [sp, phi, r, theta, vpar, mu] layout that is going to be favored in gyselalibxx. We have some machinery that setup input configuration data depending on the geometry. These are in `collision_configuration_sprvparmu.hpp`, `collision_configuration_spvparmu.hpp`. diff --git a/src/collisions/collision_common_configuration.hpp b/src/collisions/collision_common_configuration.hpp new file mode 100644 index 000000000..833e05558 --- /dev/null +++ b/src/collisions/collision_common_configuration.hpp @@ -0,0 +1,394 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include <paraconf.h> + +#include "ddc_alias_inline_functions.hpp" +#include "ddc_aliases.hpp" +#include "paraconfpp.hpp" +#include "species_info.hpp" + +namespace detail { + +/** + * @brief Not all simulations expose r,theta,phi (tor1, tor2, tor3) so we + * "spoof" the presence of the dimensions. This is the class from which fake + * dimensions inherit. These fake dimensions are inserted in order to match + * Koliop's interface. + * These fake dimensions are then assigned to an: + * IdxRange<Grid1D>(Idx<Grid1D> {0}, IdxStep<Grid1D> {1}); + * See create_scalar_index_range(); + */ +struct InternalSpoofGrid +{ +}; + +/** + * @brief Fake radial dimension to be used if there is no radial dimension in + * the simulation. + */ +struct InternalSpoofGridR : InternalSpoofGrid +{ +}; + +/** + * @brief Fake poloidal dimension to be used if there is no poloidal dimension + * in the simulation. + */ +struct InternalSpoofGridTheta : InternalSpoofGrid +{ +}; + +/** + * @brief Fake toroidal dimension to be used if there is no toroidal dimension + * in the simulation. + */ +struct InternalSpoofGridPhi : InternalSpoofGrid +{ +}; + +/** + * @brief A storage class for the fields the collision operator requires. + */ +template < + class IdxRangeDistributionFunction, + class GridSp, + class GridPhi, + class GridTheta, + class GridR, + class GridVpar, + class GridMu> +class CollisionConfigurationData +{ +public: + /** + * @brief Distribution function index range. + */ + using IdxRangeDistributionFunctionType = IdxRangeDistributionFunction; + + /** + * @brief Mu grid. + */ + using GridMuType = GridMu; + /** + * @brief Vpar grid. + */ + using GridVparType = GridVpar; + /** + * @brief R grid. + */ + using GridRType = GridR; + /** + * @brief Theta grid. + */ + using GridThetaType = GridTheta; + /** + * @brief Phi grid. + */ + using GridPhiType = GridPhi; + /** + * @brief Sp grid. + */ + using GridSpType = GridSp; + + /** + * @brief Sp,Theta,R,Vpar index range. + */ + using IdxRangeSpThetaRVparType = IdxRange<GridSpType, GridThetaType, GridRType, GridVparType>; + /** + * @brief Theta,R index range. + */ + using IdxRangeThetaRType = IdxRange<GridThetaType, GridRType>; + + /** + * @brief Mu index range. + */ + using IdxRangeMuType = IdxRange<GridMuType>; + /** + * @brief Vpar index range. + */ + using IdxRangeVparType = IdxRange<GridVparType>; + /** + * @brief R index range. + */ + using IdxRangeRType = IdxRange<GridRType>; + /** + * @brief Species index range. + */ + using IdxRangeSpType = IdxRange<GridSpType>; + + static_assert( + ddc::is_uniform_point_sampling_v<GridVparType>, + "The grid must be uniform in the vpar direction."); + static_assert( + ddc::is_uniform_point_sampling_v<GridVparType>, + "The grid must be uniform in the mu direction."); + +public: + /** + * @brief Construct a new Collision Configuration Data object/ + * + * @param[in] yaml_input_file simulation input file. + * @param[in] index_range_fdistribution the index range that will represent the + * distribution fuction given to the operator. + * @param[in] coeff_intdmu quadrature coefficient. + * @param[in] coeff_intdvpar quadrature coefficient. + */ + CollisionConfigurationData( + PC_tree_t const& yaml_input_file, + IdxRangeDistributionFunctionType index_range_fdistribution, + DConstField<IdxRangeMuType> coeff_intdmu, + DConstField<IdxRangeVparType> coeff_intdvpar) + : m_hat_As_allocation {extract_md_index_range<GridSpType>(index_range_fdistribution)} + , m_hat_Zs_allocation {extract_md_index_range<GridSpType>(index_range_fdistribution)} + , m_mask_buffer_r_allocation {extract_md_index_range<GridRType>(index_range_fdistribution)} + , m_mask_LIM_allocation {extract_md_index_range<GridThetaType, GridRType>( + index_range_fdistribution)} + , m_B_norm_allocation {extract_md_index_range<GridThetaType, GridRType>( + index_range_fdistribution)} + , m_Bstar_s_allocation {extract_md_index_range< + GridSpType, + GridThetaType, + GridRType, + GridVparType>(index_range_fdistribution)} + , m_mug_allocation {extract_md_index_range<GridMuType>(index_range_fdistribution)} + , m_vparg_allocation {extract_md_index_range<GridVparType>(index_range_fdistribution)} + , m_coeff_AD_allocation {extract_md_index_range<GridRType>(index_range_fdistribution)} + , m_coeff_intdmu {coeff_intdmu} + , m_coeff_intdvpar {coeff_intdvpar} + , m_collisions_interspecies {::PCpp_bool(yaml_input_file, ".Collisions.interspecies")} + , m_mu_extent {extract_md_index_range<GridMuType>(index_range_fdistribution).size()} + , m_vpar_extent {extract_md_index_range<GridVparType>(index_range_fdistribution).size()} + , m_r_extent {extract_md_index_range<GridRType>(index_range_fdistribution).size()} + , m_theta_extent {extract_md_index_range<GridThetaType>(index_range_fdistribution).size()} + , m_phi_extent {extract_md_index_range<GridPhiType>(index_range_fdistribution).size()} + , m_sp_extent {extract_md_index_range<GridSpType>(index_range_fdistribution).size()} + { + if (IdxRangeVparType {index_range_fdistribution}.size() % 2 != 0) { + throw std::runtime_error("The number of points in the vpar direction must be a " + "multiple of 2. This ensures that there is no grid point at " + "vpar=0 (this would cause division by 0)."); + } + + /* Define default array values. The user of CollisionConfigurationData + * is free to modify these. + */ + + { + IdxRangeSpType const index_range_sp + = extract_md_index_range<GridSpType>(index_range_fdistribution); + + /* Initialise the mass species. + */ + host_t<DConstFieldSp> hat_As_host + = ddc::host_discrete_space<GridSpType>().masses()[index_range_sp]; + ddc::parallel_deepcopy(get_field(m_hat_As_allocation), hat_As_host); + + /* Initialise the species charge. + */ + host_t<DConstFieldSp> hat_Zs_host + = ddc::host_discrete_space<GridSpType>().charges()[index_range_sp]; + ddc::parallel_deepcopy(get_field(m_hat_Zs_allocation), hat_Zs_host); + } + + { + /* Disabling Initialisation of the masks that have no sense here to + * 0. + * Masked if >= 0.99 + */ + ddc::parallel_fill(get_field(m_mask_buffer_r_allocation), 0.0); + ddc::parallel_fill(get_field(m_mask_LIM_allocation), 0.0); + } + + { + /* Initialize grid coordinates. + */ + ddcHelper:: + dump_coordinates(Kokkos::DefaultExecutionSpace(), get_field(m_mug_allocation)); + ddcHelper::dump_coordinates( + Kokkos::DefaultExecutionSpace(), + get_field(m_vparg_allocation)); + } + + Kokkos::fence(); + } + + /** + * @brief We have some fields whose dimension vary depending on the + * simulation. Bstar_s in a typical 5D simulation has dimension [Sp, theta, + * R, Vpar]. Now some simulation do not have all these dimension, say only + * [Sp, Vpar, Mu]. We have to satisfy both situations. TransposeDeepCopy + * takes two fields, expecting the source field to be copied into the + * destination field. It may be that the source has less dimension but never + * more. If the source has less dimensions, the data is replicated. For + * instance from [Vpar, Sp] into [Sp, R, Vpar], the [Vpar, Sp] plane with be + * replicated along Rand the transpose done appropriately. + * + * @param[out] destination a field. + * @param[in] source a field. + */ + template <class IdxRange0, class IdxRange1, class IdxRangeSlice = IdxRange1> + static void transpose_replicate_deep_copy( + DField<IdxRange0> destination, + DConstField<IdxRange1> source) + { + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + get_idx_range(destination), + KOKKOS_LAMBDA(typename IdxRange0::discrete_element_type destination_index) { + const typename IdxRangeSlice::discrete_element_type source_index( + destination_index); + destination(destination_index) = source(source_index); + }); + } + +protected: + /** + * @brief Check if a dimension is spoofed but is not present in the actual + * simulation. + */ + template <class Grid1D> + static constexpr bool is_spoofed_dimension_v = std::is_base_of_v<InternalSpoofGrid, Grid1D>; + + /** + * @brief Create a scalar index range object. + * + * @tparam Grid1D + * @return IdxRange<Grid1D> + */ + template <class Grid1D> + static IdxRange<Grid1D> create_scalar_index_range() + { + return IdxRange<Grid1D>(Idx<Grid1D> {0}, IdxStep<Grid1D> {1}); + } + + /** + * @brief Get the index range for a specific grid from a multi-D index range. + * If the dimension is spoofed and does not appear in the multi-D index range + * then an index range which only iterates over the index(0) is returned. + * + * @tparam Grid The tag for the specific grid. + * @param[in] index_range The multi-D index range. + * @returns The index range for the specific grid. + */ + template <class Grid1D, class IdxRange0> + static IdxRange<Grid1D> extract_1d_index_range(IdxRange0 index_range) + { + if constexpr (is_spoofed_dimension_v<Grid1D>) { + return create_scalar_index_range<Grid1D>(); + } else { + return ddc::select<Grid1D>(index_range); + } + } + + /** + * @brief Get the index range for specific grid dimensions from a multi-D index + * range. + * + * @tparam Grid1D The tags for the specific grid dimensions. + * @param[in] index_range The multi-D index range. + * @returns The index range for the specific grid. + */ + template <class... Grid1D, class IdxRange0> + static IdxRange<Grid1D...> extract_md_index_range(IdxRange0 index_range) + { + return IdxRange<Grid1D...>(extract_1d_index_range<Grid1D>(index_range)...); + } + + /* Not protected nor private because there is no internal state to + * protect/condition to satisfy. + */ +public: + /** + * @brief Normalised masses for all species. + */ + DFieldMem<IdxRangeSpType> m_hat_As_allocation; + + /** + * @brief Normalised charges for all species. + */ + DFieldMem<IdxRangeSpType> m_hat_Zs_allocation; + + /** + * @brief Mask used to avoid applying collision in certain region. + * [TODO]: This mask is not exposed to the user in the C++ version. + */ + DFieldMem<IdxRangeRType> m_mask_buffer_r_allocation; + + /** + * @brief Limiter mask in (theta,r). + */ + DFieldMem<IdxRangeThetaRType> m_mask_LIM_allocation; + + /** + * @brief B norm in (theta,r). + * [TODO] Attention this must be 3D for generalisation to 3D geometry. + */ + DFieldMem<IdxRangeThetaRType> m_B_norm_allocation; + + /** + * @brief Bstar_s(species,theta,r,vpar). + * [TODO] Must be 5D for full 3D geometry. + */ + DFieldMem<IdxRangeSpThetaRVparType> m_Bstar_s_allocation; + + /** + * @brief Grid in mu direction. + */ + DFieldMem<IdxRangeMuType> m_mug_allocation; + + /** + * @brief Grid in vpar direction. + */ + DFieldMem<IdxRangeVparType> m_vparg_allocation; + + /** + * @brief Radial coefficients of AD, almost equivalent to the collision + * frequency. + */ + DFieldMem<IdxRangeRType> m_coeff_AD_allocation; + + /** + * @brief Quadrature coefficients in mu direction. + */ + DConstField<IdxRangeMuType> m_coeff_intdmu; + + /** + * @brief Quadrature coefficients in vpar direction. + */ + DConstField<IdxRangeVparType> m_coeff_intdvpar; + + /** + * @brief Boolean that is equal to true if inter-species collisions are + * taken into account. + * Read in the YAML input file. + */ + bool m_collisions_interspecies; + + /** + * @brief Mu dimension extent. + */ + std::size_t m_mu_extent; + /** + * @brief Vpar dimension extent. + */ + std::size_t m_vpar_extent; + /** + * @brief R dimension extent. + */ + std::size_t m_r_extent; + /** + * @brief Theta dimension extent. + */ + std::size_t m_theta_extent; + /** + * @brief Phi dimension extent. + */ + std::size_t m_phi_extent; + /** + * @brief Species dimension extent. + */ + std::size_t m_sp_extent; +}; +} // namespace detail diff --git a/src/collisions/collision_operator.cpp b/src/collisions/collision_operator.cpp new file mode 100644 index 000000000..e51bb6c1b --- /dev/null +++ b/src/collisions/collision_operator.cpp @@ -0,0 +1,86 @@ +#include <KOLIOP/koliop.h> + +#include "assert.hpp" +#include "collision_operator.hpp" + +// NOTE: Offset of in the r dimension. Used in the Fortran code. +static inline constexpr std::size_t local_idx_range_r_offset = 0; +// NOTE: Offset of in the theta dimension. Used in the Fortran code. +static inline constexpr std::size_t local_idx_range_theta_offset = 0; + +// NOTE: For now, these are constant. When more physics get implemented, pass +// them to the Collisions operator's constructor. +// NOTE: Enable masking. +static inline constexpr std::int8_t is_sol_enabled = false; +static inline constexpr std::int8_t is_limiter_mask_enabled = true; + +namespace detail { + +::koliop_Operator do_operator_initialisation( + std::size_t mu_extent, + std::size_t vpar_extent, + std::size_t tor1_extent, + std::size_t tor2_extent, + std::size_t tor3_extent, + std::size_t species_extent, + std::int8_t collision_interspecies, + std::int64_t ir_sol_separatrix, + double const* hat_As, + double const* hat_Zs, + double const* mug, + double const* vparg, + double const* coeff_intdmu, + double const* coeff_intdvpar, + double const* coeff_AD_r, + double const* mask_buffer_r, + double const* mask_limiter, + double const* B_norm, + double const* Bstar_s) +{ + ::koliop_Operator operator_handle; + + std::size_t const r_extent = tor2_extent; + std::size_t const theta_extent = tor1_extent; + + if (::koliop_Create( + &operator_handle, + mu_extent, + vpar_extent, + tor1_extent, + tor2_extent, + tor3_extent, + species_extent, + r_extent, // NOTE: No MPI index range decomposition. + theta_extent, // NOTE: No MPI index range decomposition. + local_idx_range_r_offset, + local_idx_range_theta_offset, + collision_interspecies, + is_sol_enabled, + is_limiter_mask_enabled, + ir_sol_separatrix, + hat_As, + hat_Zs, + mug, + vparg, + coeff_intdmu, + coeff_intdvpar, + coeff_AD_r, + mask_buffer_r, + mask_limiter, + B_norm, + Bstar_s) + != KOLIOP_STATUS_SUCCESS) { + GSLX_ASSERT(false); + } + + return operator_handle; +} + +void do_operator_deinitialisation(::koliop_Operator operator_handle) +{ + if (::koliop_Release(operator_handle) != KOLIOP_STATUS_SUCCESS) { + GSLX_ASSERT(false); + } +} + +} // namespace detail diff --git a/src/collisions/collision_operator.hpp b/src/collisions/collision_operator.hpp new file mode 100644 index 000000000..39c659471 --- /dev/null +++ b/src/collisions/collision_operator.hpp @@ -0,0 +1,169 @@ +// SPDX-License-Identifier: MIT +#pragma once + +#include <ddc/ddc.hpp> + +#include "assert.hpp" +#include "ddc_alias_inline_functions.hpp" +#include "ddc_aliases.hpp" +#include "ddc_helper.hpp" +#include "species_info.hpp" + +// SPDX-License-Identifier: MIT +#pragma once + +#include <ddc/ddc.hpp> + +#include <KOLIOP/koliop.h> + +namespace detail { + +/** + * @brief Initialise the collision operator defined in koliop. + * + * @param[in] mu_extent The number of points in the magnetic moment dimension. + * @param[in] vpar_extent The number of points in the parallel velocity dimension. + * @param[in] tor1_extent The number of points in the tor1 dimension. + * @param[in] tor2_extent The number of points in the tor2 dimension. + * @param[in] tor3_extent The number of points in the tor3 dimension. + * @param[in] species_extent The number of species. + * @param[in] collision_interspecies Boolean that is equal to true if inter-species collisions are taken into account. + * @param[in] ir_sol_separatrix The index in the radial dimension where the SOL separatrix is found. + * @param[in] hat_As The normalised masses for all species. + * @param[in] hat_Zs The normalised charges for all species. + * @param[in] mug The coordinates of the grid of magnetic moments. + * @param[in] vparg The coordinates of the grid of parallel velocities. + * @param[in] coeff_intdmu The quadrature coefficients to integrate on the grid of magnetic moments. + * @param[in] coeff_intdvpar The quadrature coefficients to integrate on the grid of parallel velocities. + * @param[in] coeff_AD_r + * @param[in] mask_buffer_r The mask used to avoid applying collisions in certain radial region. + * @param[in] mask_limiter The limiter mask defined over the polar slice. + * @param[in] B_norm The norm of the magnetic field defined over the polar slice. + * @param[in] Bstar_s Bstar defined on the coordinates (species,r,theta,vpar). + * + * @returns The handle which identifies the operator. + */ +::koliop_Operator do_operator_initialisation( + std::size_t mu_extent, + std::size_t vpar_extent, + std::size_t tor1_extent, + std::size_t tor2_extent, + std::size_t tor3_extent, + std::size_t species_extent, + std::int8_t collision_interspecies, + std::int64_t ir_sol_separatrix, + double const* hat_As, + double const* hat_Zs, + double const* mug, + double const* vparg, + double const* coeff_intdmu, + double const* coeff_intdvpar, + double const* coeff_AD_r, + double const* mask_buffer_r, + double const* mask_limiter, + double const* B_norm, + double const* Bstar_s); + +/** + * @brief Destructor for koliop. + * + * @param[in] operator_handle The handle which identifies the operator. + */ +void do_operator_deinitialisation(::koliop_Operator operator_handle); + +}; // namespace detail + +/** + * @brief A class which computes the collision operator in (Sp,vpar,mu). + */ +template <class CollisionConfiguration> +class CollisionOperator /* : public IRightHandSide */ +{ +public: + /** + * @brief The type of collision configuration. Namely, what kind of geometry + * the operator will ... operate on. + */ + using CollisionConfigurationType = CollisionConfiguration; + /** + * @brief The distribution function that we should expect from the operator + * user/caller. + */ + using IdxRangeDistributionFunctionType = + typename CollisionConfigurationType::IdxRangeDistributionFunctionType; + +public: + /** + * @brief Construct a new Collision Operator object + * + * @param[in] collision_configuration This parametrise the operator. It contains the + * data required for its initialization. + */ + explicit CollisionOperator(CollisionConfigurationType const& collision_configuration) + : m_operator_handle {detail::do_operator_initialisation( + collision_configuration.configuration().m_mu_extent, + collision_configuration.configuration().m_vpar_extent, + collision_configuration.configuration().m_r_extent, + collision_configuration.configuration().m_theta_extent, + collision_configuration.configuration().m_phi_extent, + collision_configuration.configuration().m_sp_extent, + collision_configuration.configuration().m_collisions_interspecies, + /* One less than the number of point in the r direction. */ + 0 /* Null offset */ + collision_configuration.configuration().m_r_extent - 1, + collision_configuration.configuration().m_hat_As_allocation.data_handle(), + collision_configuration.configuration().m_hat_Zs_allocation.data_handle(), + collision_configuration.configuration().m_mug_allocation.data_handle(), + collision_configuration.configuration().m_vparg_allocation.data_handle(), + collision_configuration.configuration().m_coeff_intdmu.data_handle(), + collision_configuration.configuration().m_coeff_intdvpar.data_handle(), + collision_configuration.configuration().m_coeff_AD_allocation.data_handle(), + collision_configuration.configuration().m_mask_buffer_r_allocation.data_handle(), + collision_configuration.configuration().m_mask_LIM_allocation.data_handle(), + collision_configuration.configuration().m_B_norm_allocation.data_handle(), + collision_configuration.configuration().m_Bstar_s_allocation.data_handle())} + { + // EMPTY + } + + /** + * @brief Destroy the Collision operator object. + */ + ~CollisionOperator() + { + detail::do_operator_deinitialisation(m_operator_handle); + } + + /** + * @brief Apply the collision operator to the distribution functions of all + * species on all species. + * + * @param[inout] all_f_distribution All the distribution function, depending + * on the CollisionConfigurationType, the existence of the theta, r, phi + * dimension may vary. At most, we have (species, phi, r, theta, vpar, mu) + * in layout right. + * @param[in] deltat_coll Collision time step. + */ + void operator()(DField<IdxRangeDistributionFunctionType> all_f_distribution, double deltat_coll) + const + { + if (::koliop_Collision( + static_cast<::koliop_Operator>(m_operator_handle), + deltat_coll, + all_f_distribution.data_handle()) + != KOLIOP_STATUS_SUCCESS) { + GSLX_ASSERT(false); + } + + // NOTE: Koliop is asynchronous, fence to ensure the operator ended. + if (::koliop_Fence(static_cast<::koliop_Operator>(m_operator_handle)) + != KOLIOP_STATUS_SUCCESS) { + GSLX_ASSERT(false); + } + } + +protected: + /** + * @brief Opaque C type representing the operator. + */ + ::koliop_Operator m_operator_handle; +}; diff --git a/src/collisions/collisioninfo.hpp b/src/collisions/collisioninfo.hpp deleted file mode 100644 index 65651073b..000000000 --- a/src/collisions/collisioninfo.hpp +++ /dev/null @@ -1,71 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once - -#include <paraconf.h> - -#include "paraconfpp.hpp" -#include "species_info.hpp" - -/// Class to collect information to initialize the collision operator -class CollisionInfo -{ - // value of nustar0 = to nustar0_rpeak read in the YAML input file - double const m_nustar0; - - // boolean that is equal to true if inter-species collisions are taken into account - // read in the YAML input file - std::int8_t m_collisions_interspecies; - - // AD coefficient almost equivalent to the collision frequency - double m_coeff_AD; - -public: - /// radial_chunk_type used to treat the 0D case for radial profile - using radial_chunk_type = double; - - /** - * @brief The constructor for the CollisionFrequency class. - * @param[in] yaml_input_file YAML input file containing CollisionsInfo - */ - explicit CollisionInfo(PC_tree_t const& yaml_input_file) - : m_nustar0 {PCpp_double(yaml_input_file, ".CollisionsInfo.nustar0_rpeak")} - , m_collisions_interspecies { - PCpp_bool(yaml_input_file, ".CollisionsInfo.collisions_interspecies")} - { - // nustar fixed to 1. => normalisation of time is equal to 1./nustar - if (m_nustar0 != 1.0) { - throw std::invalid_argument("nustar0 must be equal to 1"); - } - - /** - * Computation of the radial profile of ADi that is required for koliop interface - * @f[ AD(rpeak) = \sqrt(2)*eps_rpeak^(3/2)/(q_rpeak R0)*nustar0_rpeak @f] - * eps_rpeak = eps(rpeak) with eps(r) = r/R0 is the aspect ratio. - * - * In this specific case, where AD is a scalar all values are forced to 1, so - * coeff_AD_atrpeak = sqrt(2) - **/ - m_coeff_AD = Kokkos::sqrt(2.); - }; - - ~CollisionInfo() = default; - - /** - * @brief A method for accessing the collisions_interspecies member variable of the class. - * @return A boolean containing the collisions_interspecies value. - */ - std::int8_t collisions_interspecies() const - { - return m_collisions_interspecies; - } - - /** - * @brief A method for accessing the coeff_AD variable of the class. - * @return A double containing the coeff_AD value. - */ - double coeff_AD() const - { - return m_coeff_AD; - } -}; diff --git a/src/collisions/collisioninfo_radial.hpp b/src/collisions/collisioninfo_radial.hpp deleted file mode 100644 index 48b4162bd..000000000 --- a/src/collisions/collisioninfo_radial.hpp +++ /dev/null @@ -1,130 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once -#include <paraconf.h> - -#include "ddc_alias_inline_functions.hpp" -#include "ddc_aliases.hpp" -#include "paraconfpp.hpp" -#include "species_info.hpp" - -/// Class to collect information to initialize the collision operator -template <class GridR> -class CollisionInfoRadial -{ -private: - using IdxR = Idx<GridR>; - using IdxRangeR = IdxRange<GridR>; - -public: - /// Type alias for a field memory block on a grid of radial values - using DFieldMemR = DFieldMem<IdxRangeR>; - /// Type alias for a field defined on a grid of radial values - using DFieldR = DField<IdxRangeR>; - /// Type alias for a field defined on a grid of radial values - using DConstFieldR = DConstField<IdxRangeR>; - /// radial_chunk_type used to treat the 0D case for radial profile - using radial_chunk_type = DConstFieldR; - -private: - /// value of nustar0_rpeak read in the YAML input file - double m_nustar0_rpeak; - - /// boolean that is equal to true if inter-species collisions are taken into account - // read in the YAML input file - std::int8_t m_collisions_interspecies; - - // [TODO] Values that are temporarily need for the interface with koliop - // but that will be deleted as soon as no more required - double const m_rpeak; - double const m_q_rpeak; - double const m_R0; - - DConstFieldR m_rg; //VG ? - DConstFieldR m_safety_factor; //VG? - - /// radial coefficients of AD - DFieldMemR m_coeff_AD; - -public: - /** - * Computation of the radial profile of AD - * @f[ AD(rpeak) = \sqrt(2)*eps_rpeak^(3/2)/(q_rpeak R0)*nustar0_rpeak @f] - * where R0, q_rpeak, nustar0_rpeak are given as input data and - * eps_rpeak = eps(rpeak) with eps(r) = r/R0 is the aspect ratio. - * - * @param[in] idxrange_r radial index range - * - * This function should be private but is public due to Kokkos restrictions. - */ - void compute_coeff_AD(IdxRangeR idxrange_r) - { - double const rpeak = m_rpeak; - double const q_rpeak = m_q_rpeak; - double const R0 = m_R0; - double const nustar0_rpeak = m_nustar0_rpeak; - - // Compute coeff_AD_rpeak - DFieldR coeff_AD = get_field(m_coeff_AD); - double const eps_rpeak = rpeak / R0; - double const coeff_AD_rpeak - = std::sqrt(2.) * eps_rpeak * std::sqrt(eps_rpeak) * nustar0_rpeak / (q_rpeak * R0); - - // [TODO] coeff_AD should be a scalar and not a radial profile - // but to change that koliop interface must be changed - // so coeff_AD(r) is fixed equal to coeff_AD_rpeak - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - idxrange_r, - KOKKOS_LAMBDA(Idx<GridR> idx_r) { coeff_AD(idx_r) = coeff_AD_rpeak; }); - } - -public: - /** - * @brief The constructor for the CollisionInfoRadial class. - * @param[in] rpeak the radial reference position r = rpeak - * @param[in] q_rpeak the safety factor value at r = rpeak - * @param[in] R0 the major radius at the magnetic axis - * @param[in] nustar0_rpeak the value of nustar at r = rpeak - * @param[in] collisions_interspecies logical value equal to true if interspecies collisions are taken into account - * @param[in] idxrange_r radial index range - */ - CollisionInfoRadial( - double const rpeak, - double const q_rpeak, - double const R0, - double const nustar0_rpeak, - std::int8_t const collisions_interspecies, - IdxRangeR idxrange_r) - : m_rpeak(rpeak) - , m_q_rpeak(q_rpeak) - , m_R0(R0) - , m_nustar0_rpeak(nustar0_rpeak) - , m_collisions_interspecies(collisions_interspecies) - , m_coeff_AD(idxrange_r) - { - compute_coeff_AD(idxrange_r); - }; - - ~CollisionInfoRadial() = default; - - - /** - * @brief A method for accessing the collisions_interspecies member variable of the class. - * @return A boolean containing the collisions_interspecies value. - */ - std::int8_t collisions_interspecies() const - { - return m_collisions_interspecies; - } - - - /** - * @brief A method for accessing coeff_AD (the radial profile of AD coeff) variable of the class. - * @return A field containing the radial profile of the AD coefficients. - */ - DConstFieldR coeff_AD() const - { - return get_const_field(m_coeff_AD); - } -}; diff --git a/src/collisions/collisions_dimensions.hpp b/src/collisions/collisions_dimensions.hpp deleted file mode 100644 index ec13c2702..000000000 --- a/src/collisions/collisions_dimensions.hpp +++ /dev/null @@ -1,260 +0,0 @@ -// SPDX-License-Identifier: MIT -#pragma once -#include "ddc_alias_inline_functions.hpp" -#include "ddc_aliases.hpp" - -/** - * @brief A namespace to collect classes which are necessary to create Fields with the - * correct number of dimensions to be compatible with Koliop. - */ -namespace collisions_dimensions { -/** - * Create dimensions to act as the radial/poloidal/toroidal tag for Koliop even if these dims don't - * exist in the simulation - */ - -/// Class from which fake dimensions inherit. These fake dimensions are inserted in order to match Koliop interface -struct InternalSpoofGrid -{ -}; - -/// Fake radial dimension to be used if there is no radial dimension in the simulation -struct InternalSpoofGridR : InternalSpoofGrid -{ -}; - -/// Fake poloidal dimension to be used if there is no poloidal dimension in the simulation -struct InternalSpoofGridTheta : InternalSpoofGrid -{ -}; - -/// Check if a dimension is spoofed but is not present in the actual simulation -template <class Grid1D> -inline constexpr bool is_spoofed_dim_v = std::is_base_of_v<InternalSpoofGrid, Grid1D>; - -/** - * Class to get the type of the radial dimension from a field containing a radial profile. - * @tparam Field The type of the field containing the radial profile. - */ -template <class Field> -struct ExtractRDim -{ - static_assert(!std::is_same_v<Field, Field>, "Unrecognised radial profile type"); -}; - -/** - * Class to get the type of the poloidal dimension from a field containing a profile on the poloidal plane. - * @tparam Field The type of the field containing the profile on the poloidal plane. - * @tparam GridR The tag for the discrete radial dimension. - */ -template <class Field, class GridR> -struct ExtractThetaDim -{ - static_assert(!std::is_same_v<Field, Field>, "Unrecognised poloidal profile type"); -}; - -/** - * @brief Get the index range for a specific grid from a multi-D index range. - * If the dimension is spoofed and does not appear in the multi-D index range then an index - * range which only iterates over the index(0) is returned. - * - * @tparam Grid The tag for the specific grid. - * @param idx_range The multi-D index range. - * @returns The index range for the specific grid. - */ -template <class Grid1D, class IdxRangeFDistrib> -inline IdxRange<Grid1D> get_1d_idx_range(IdxRangeFDistrib idx_range) -{ - if constexpr (is_spoofed_dim_v<Grid1D>) { - return IdxRange<Grid1D>(Idx<Grid1D> {0}, IdxStep<Grid1D> {1}); - } else { - return ddc::select<Grid1D>(idx_range); - } -} - -/** - * @brief Get the index range for specific grid dimensions from a multi-D index range. - * - * @tparam Grid1D The tags for the specific grid dimensions. - * @param idx_range The multi-D index range. - * @returns The index range for the specific grid. - */ -template <class... Grid1D, class IdxRangeFDistrib> -inline IdxRange<Grid1D...> get_expanded_idx_range(IdxRangeFDistrib idx_range) -{ - return IdxRange<Grid1D...>(get_1d_idx_range<Grid1D>(idx_range)...); -} - -/// If radial profile is stored in a double then the grid tag must be spoofed. -template <> -struct ExtractRDim<double> -{ - using type = InternalSpoofGridR; -}; - -/// If radial profile is stored in a 1D chunk then the grid tag is extracted. -template <class GridR, class Layout> -struct ExtractRDim< - ConstField<double, IdxRange<GridR>, Kokkos::DefaultExecutionSpace::memory_space, Layout>> -{ - using type = GridR; -}; - -/// If radial profile is stored in a chunk then information about why the chunk has the wrong type is provided -template <class ElementType, class IdxRange, class MemSpace, class Layout> -struct ExtractRDim<Field<ElementType, IdxRange, MemSpace, Layout>> -{ - static_assert( - std::is_same_v<ElementType, const double>, - "The radial profile should be a double and should not be modifiable. Please use a " - "constant field."); - static_assert( - std::is_same_v<MemSpace, Kokkos::DefaultExecutionSpace::memory_space>, - "The radial profile should be provided on the GPU"); - static_assert( - (ddc::type_seq_size_v<ddc::to_type_seq_t<IdxRange>>) > 1, - "The radial profile should not be defined on more than 1 dimensions."); -}; - -/// If the profile on the poloidal plane is stored in a double then the grid tag must be spoofed. -template <> -struct ExtractThetaDim<double, InternalSpoofGridR> -{ - using type = InternalSpoofGridTheta; -}; - -/// If the profile on the poloidal plane is stored in a chunk on radial values then the grid tag must be spoofed. -template <class GridR, class Layout> -struct ExtractThetaDim< - ConstField<double, IdxRange<GridR>, Kokkos::DefaultExecutionSpace::memory_space, Layout>, - GridR> -{ - using type = InternalSpoofGridTheta; -}; - -/// If the profile on the poloidal plane is stored in a chunk on poloidal values then the grid tag is extracted. -template <class GridR, class GridTheta, class Layout> -struct ExtractThetaDim< - ConstField< - double, - IdxRange<GridTheta>, - Kokkos::DefaultExecutionSpace::memory_space, - Layout>, - GridR> -{ - using type = GridTheta; -}; - -/// If the profile on the poloidal plane is stored in a 2D chunk then the grid tag is extracted. -template <class GridR, class GridTheta, class Layout> -struct ExtractThetaDim< - ConstField< - double, - IdxRange<GridTheta, GridR>, - Kokkos::DefaultExecutionSpace::memory_space, - Layout>, - GridR> -{ - using type = GridTheta; -}; - -/// If poloidal profile is stored in a chunk then information about why the chunk has the wrong type is provided -template <class GridR, class ElementType, class IdxRange, class MemSpace, class Layout> -struct ExtractThetaDim<Field<ElementType, IdxRange, MemSpace, Layout>, GridR> -{ - static_assert( - std::is_same_v<ElementType, const double>, - "The poloidal profile should be a double and should not be modifiable. Please use a " - "constant field."); - static_assert( - std::is_same_v<MemSpace, Kokkos::DefaultExecutionSpace::memory_space>, - "The poloidal profile should be provided on the GPU"); - static_assert( - (ddc::type_seq_size_v<ddc::to_type_seq_t<IdxRange>>) > 2, - "The poloidal profile should not be defined on more than 2 dimensions."); -}; - -/** - * Class to get a DDC TypeSeq containing the input radial and theta tags. - * This is useful to construct other related input types. - * - * @tparam The type of the field containing the profile on the poloidal plane. - */ -template <class Field> -struct ExtractRThetaTags -{ - static_assert(!std::is_same_v<Field, Field>, "Unrecognised poloidal profile type"); -}; - -/// If the profile on the poloidal plane is stored in a double then there are no input tags -template <> -struct ExtractRThetaTags<double> -{ - using type = ddc::detail::TypeSeq<>; -}; - -/// If the profile on the poloidal plane is stored in a ConstField then the input tags are extracted from the index range -template <class IdxRange> -struct ExtractRThetaTags<DConstField<IdxRange>> -{ - using type = ddc::to_type_seq_t<IdxRange>; -}; - -/** - * A helper function to check if a grid is found at the correct location in a set of - * ordered grids. - * - * @tparam Grid1D The grid we are checking. - * @tparam OrderedGrids A ddc::TypeSeq containing an ordered set of grids. - * - * @param idx The index where the dimension is expected to be found - * - * @returns True if the grid is at the expected location or if the grid is spoofed. - * False otherwise. - */ -template <class Grid1D, class OrderedGrids> -constexpr bool check_dimension_location(int& idx) -{ - if constexpr (!is_spoofed_dim_v<Grid1D>) { - bool success = ddc::type_seq_rank_v<Grid1D, OrderedGrids> == idx; - idx += 1; - return success; - } else { - return true; - } -} - -/** - * A helper function to check if the final grids in an ordered set of grids match - * the order of the provided expected grids. - * - * @tparam OrderedGrids A ddc::TypeSeq containing an ordered set of grids. - * @tparam ExpectedGrids The expected grids in the order that they are expected. - * - * @returns True if the grids are ordered correctly, false otherwise. - */ -template <class OrderedGrids, class... ExpectedGrids> -constexpr bool order_of_last_grids() -{ - int n_real_dims = ((is_spoofed_dim_v<ExpectedGrids> ? 0 : 1) + ...); - int idx(ddc::type_seq_size_v<OrderedGrids> - n_real_dims); - return ((check_dimension_location<ExpectedGrids, OrderedGrids>(idx)) && ...); -} - -/** - * A helper function to check if the first grids in an ordered set of grids match - * the order of the provided expected grids. - * - * @tparam OrderedGrids A ddc::TypeSeq containing an ordered set of grids. - * @tparam ExpectedGrids The expected grids in the order that they are expected. - * - * @returns True if the grids are ordered correctly, false otherwise. - */ -template <class OrderedGrids, class... ExpectedGrids> -constexpr bool order_of_first_grids() -{ - int idx(0); - return ((check_dimension_location<ExpectedGrids, OrderedGrids>(idx)) && ...); -} - -} // namespace collisions_dimensions diff --git a/src/collisions/koliop_interface.cpp b/src/collisions/koliop_interface.cpp deleted file mode 100644 index 2441bee9f..000000000 --- a/src/collisions/koliop_interface.cpp +++ /dev/null @@ -1,128 +0,0 @@ -#include <atomic> -#include <cassert> -#include <iostream> - -#include <KOLIOP/koliop.h> - -#include "assert.hpp" -#include "koliop_interface.hpp" - -// NOTE: For now, these are constant. When more physics get implemented, pass -// them via the Collisions' constructor. -static inline constexpr std::size_t Npolmax = 3; -static inline constexpr std::size_t index_max = 2 * (Npolmax - 1); -// pSpecies_id represents the species processed by a rank. -// It is fixed to 0 but is no longer used since the 2024/03/01 modification in the Fortran version. -// After this date and on both the C++ and Fortran versions, the value is meaningless as a rank contains all -// the species for a given global simulation space sub volume. -// [TODO] Delete this input (both for Fortran and C++ versions) -static inline constexpr std::int64_t pSpecies_id = 0; -static inline constexpr std::int64_t mu_id - = 0; // [TODO] Same than for pSpecies_id: must be deleted from the list of input -static inline constexpr std::size_t the_local_idx_range_r_offset = 0; -static inline constexpr std::size_t the_local_idx_range_theta_offset = 0; -static inline constexpr std::int8_t SOL = false; -static inline constexpr std::int8_t LIM = false; - -namespace koliop_interface { - -::koliop_Operator DoOperatorInitialization( - std::size_t the_mu_extent, - std::size_t the_vpar_extent, - std::size_t the_r_extent, - std::size_t the_theta_extent, - std::size_t the_phi_extent, - std::size_t the_species_extent, - std::int8_t collision_interspecies, - std::int64_t ir_SOL_separatrix, - const double* hat_As, - const double* hat_Zs, - const double* mug, - const double* vparg, - const double* coeff_intdmu, - const double* coeff_intdvpar, - const double* coeff_AD_r, - const double* mask_buffer_r, - const double* mask_LIM, - const double* B_norm, - const double* Bstar_s) -{ - ::koliop_Operator the_operator_handle; - - if (::koliop_Create( - &the_operator_handle, - the_mu_extent, - the_vpar_extent, - the_r_extent, - the_theta_extent, - the_phi_extent, - the_species_extent, - the_r_extent, // NOTE: No MPI index range decomposition. - the_theta_extent, // NOTE: No MPI index range decomposition. - the_local_idx_range_r_offset, - the_local_idx_range_theta_offset, - collision_interspecies, - SOL, - LIM, - ir_SOL_separatrix, - hat_As, - hat_Zs, - mug, - vparg, - coeff_intdmu, - coeff_intdvpar, - coeff_AD_r, - mask_buffer_r, - mask_LIM, - B_norm, - Bstar_s) - != KOLIOP_STATUS_SUCCESS) { - GSLX_ASSERT(false); - } - - return the_operator_handle; -} - -void DoOperatorDeinitialization(::koliop_Operator the_operator_handle) -{ - if (::koliop_Release(the_operator_handle) != KOLIOP_STATUS_SUCCESS) { - GSLX_ASSERT(false); - } -} - -/** - * UHL is a helper type which avoids repetition of a lengthy type name. - * In the future this type may be removed once the objects which use these - * types are defined elsewhere with DDC objects. - * UHL stands for: - * U : Unmanaged - This indicates that the type does not own its data and will not deallocate it. - * H : Host - The memory is saved on CPU. - * L : Left memory layout - The multi-indices are strided left-to-right (meaning indexing matches - * the standard indexing in Fortran not C) - */ -template <typename Extent> -using UHL = Kokkos::View< - Extent, - Kokkos::LayoutLeft, - Kokkos::HostSpace, - Kokkos::MemoryTraits<Kokkos::Unmanaged | Kokkos::Restrict>>; - -void DoCombMatComputation(MDL<double[6][6]>& comb_mat) -{ - // NOTE: In layout left. - // Precompute: n! /(k! * (n-k)!) - - // clang-format off - static constexpr double kCombMat[6 * 6] - = {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, - 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, - 1.0, 2.0, 1.0, 0.0, 0.0, 0.0, - 1.0, 3.0, 3.0, 1.0, 0.0, 0.0, - 1.0, 4.0, 6.0, 4.0, 1.0, 0.0, - 1.0, 5.0, 10.0, 10.0, 5.0, 1.0}; - // clang-format on - - Kokkos::deep_copy(comb_mat, UHL<const double**>(kCombMat, 6, 6)); -} - -} // namespace koliop_interface diff --git a/src/collisions/koliop_interface.hpp b/src/collisions/koliop_interface.hpp deleted file mode 100644 index e30628ac3..000000000 --- a/src/collisions/koliop_interface.hpp +++ /dev/null @@ -1,89 +0,0 @@ -// SPDX-License-Identifier: MIT -#pragma once - -#include <ddc/ddc.hpp> - -#include <KOLIOP/koliop.h> - -/** - * @brief A namespace containing functions which call the koliop library - */ -namespace koliop_interface { -/** - * MDL is a helper type which avoids repetition of a lengthy type name. - * In the future this type may be removed once the objects which use these - * types are defined elsewhere with DDC objects. - * MDL stands for: - * M : Managed - This indicates that the type owns its data and is responsible for deallocating it. - * D : Device - The memory is saved on GPU (when compiled for GPU). - * L : Left memory layout - The multi-indices are strided left-to-right (meaning indexing matches - * the standard indexing in Fortran not C) - */ -template <typename Extent> -using MDL = Kokkos::View< - Extent, - Kokkos::LayoutLeft, - Kokkos::DefaultExecutionSpace::memory_space, - Kokkos::MemoryTraits<Kokkos::Restrict>>; - -/** - * @brief Initialise the collision operator defined in koliop. - * - * @param[in] the_mu_extent The number of points in the magnetic moment dimension. - * @param[in] the_vpar_extent The number of points in the parallel velocity dimension. - * @param[in] the_r_extent The number of points in the radial dimension. - * @param[in] the_theta_extent The number of points in the poloidal dimension. - * @param[in] the_phi_extent The number of points in the toroidal dimension. - * @param[in] the_species_extent The number of species. - * @param[in] collision_interspecies Boolean that is equal to true if inter-species collisions are taken into account. - * @param[in] ir_SOL_separatrix The index in the radial dimension where the SOL separatrix is found. - * @param[in] hat_As The normalized masses for all species. - * @param[in] hat_Zs The normalized charges for all species. - * @param[in] mug The coordinates of the grid of magnetic moments. - * @param[in] vparg The coordinates of the grid of parallel velocities. - * @param[in] coeff_intdmu The quadrature coefficients to integrate on the grid of magnetic moments. - * @param[in] coeff_intdvpar The quadrature coefficients to integrate on the grid of parallel velocities. - * @param[in] coeff_AD_r - * @param[in] mask_buffer_r The mask used to avoid applying collisions in certain radial region. - * @param[in] mask_LIM The limiter mask defined over the polar slice. - * @param[in] B_norm The norm of the magnetic field defined over the polar slice. - * @param[in] Bstar_s Bstar defined on the coordinates (species,r,theta,vpar). - * - * @returns The handle which identifies the operator. - */ -::koliop_Operator DoOperatorInitialization( - std::size_t the_mu_extent, - std::size_t the_vpar_extent, - std::size_t the_r_extent, - std::size_t the_theta_extent, - std::size_t the_phi_extent, - std::size_t the_species_extent, - std::int8_t collision_interspecies, - std::int64_t ir_SOL_separatrix, - const double* hat_As, - const double* hat_Zs, - const double* mug, - const double* vparg, - const double* coeff_intdmu, - const double* coeff_intdvpar, - const double* coeff_AD_r, - const double* mask_buffer_r, - const double* mask_LIM, - const double* B_norm, - const double* Bstar_s); - -/** - * @brief Calculate the combinatory (6x6) matrix. - * - * @param[out] comb_mat The combinatory (6x6) matrix. - */ -void DoCombMatComputation(MDL<double[6][6]>& comb_mat); - - -/** - * @brief Destructor for koliop. - * - * @param the_operator_handle The handle which identifies the operator. - */ -void DoOperatorDeinitialization(::koliop_Operator the_operator_handle); -}; // namespace koliop_interface diff --git a/src/data_types/deriv_details.hpp b/src/data_types/deriv_details.hpp index 64ab13209..2c0302945 100644 --- a/src/data_types/deriv_details.hpp +++ b/src/data_types/deriv_details.hpp @@ -139,7 +139,7 @@ KOKKOS_FUNCTION Container<HeadQueryTag, QueryTags...> select_default( * describing derivatives from a type sequence of derivative dimensions. */ template <class DerivSeq> -class DefaultDerivElem; +struct DefaultDerivElem; template <class... Tags> struct DefaultDerivElem<ddc::detail::TypeSeq<Tags...>> diff --git a/src/geometryVparMu/CMakeLists.txt b/src/geometryVparMu/CMakeLists.txt index dd8515e91..f0f54f48e 100644 --- a/src/geometryVparMu/CMakeLists.txt +++ b/src/geometryVparMu/CMakeLists.txt @@ -1,4 +1,5 @@ # SPDX-License-Identifier: MIT +add_subdirectory(collisions) add_subdirectory(geometry) add_subdirectory(initialization) diff --git a/src/geometryVparMu/README.md b/src/geometryVparMu/README.md index d8d8b311c..39a27c11c 100644 --- a/src/geometryVparMu/README.md +++ b/src/geometryVparMu/README.md @@ -2,9 +2,11 @@ The `geometryVparMu` folder contains all the code describing methods which are specific to a 2D operator in velocity space (vpar,mu). It is broken up into the following sub-folders: -- [geometry](./geometry/README.md) - All the dimension tags used for (vpar,mu) velocity space geometry description. +- [collisions](./collisions/README.md) : Configuration of the collision operator for this geometry. -- [initialization](./initialization/README.md) - Initializes the distribution function (species,vpar,mu). +- [geometry](./geometry/README.md) - All the dimension tags used for (vpar,mu) velocity space geometry description. + +- [initialization](./initialization/README.md) - Initializes the distribution function (species,vpar,mu). diff --git a/src/geometryVparMu/collisions/CMakeLists.txt b/src/geometryVparMu/collisions/CMakeLists.txt new file mode 100644 index 000000000..627b208ee --- /dev/null +++ b/src/geometryVparMu/collisions/CMakeLists.txt @@ -0,0 +1,15 @@ +# SPDX-License-Identifier: MIT + +add_library("geometry_collisions_vparmu" INTERFACE) + +target_include_directories("geometry_collisions_vparmu" + INTERFACE + "${CMAKE_CURRENT_SOURCE_DIR}" +) + +target_link_libraries("geometry_collisions_vparmu" + INTERFACE + gslx::collisions +) + +add_library("gslx::geometry_collisions_vparmu" ALIAS "geometry_collisions_vparmu") diff --git a/src/geometryVparMu/collisions/README.md b/src/geometryVparMu/collisions/README.md new file mode 100644 index 000000000..1e93d7e16 --- /dev/null +++ b/src/geometryVparMu/collisions/README.md @@ -0,0 +1,3 @@ +# CollisionConfiguration : + +The collision folder contains a configuration file `collision_configuration.hpp` which provides the quantities the operator needs. This files maps this geometry to the full 5D the Koliop operator works on. diff --git a/src/geometryVparMu/collisions/collision_configuration.hpp b/src/geometryVparMu/collisions/collision_configuration.hpp new file mode 100644 index 000000000..64283aea6 --- /dev/null +++ b/src/geometryVparMu/collisions/collision_configuration.hpp @@ -0,0 +1,192 @@ +// SPDX-License-Identifier: MIT + +#pragma once + +#include <paraconf.h> + +#include "collision_common_configuration.hpp" +#include "paraconfpp.hpp" + +/** + * @brief Class to collect information to initialise the collision operator for + * a SpVparMu geometry. + * + * NOTE: Thanks to C++17 template constructor argument type deduction, we do not + * need to give the 3 template arguments. Giving a correct + * index_range_fdistribution is sufficient. + */ +template <class GridSp, class GridVpar, class GridMu> +class CollisionConfiguration +{ +public: + /** + * @brief Mu grid. + */ + using GridMuType = GridMu; + /** + * @brief Vpar grid. + */ + using GridVparType = GridVpar; + /** + * @brief Theta grid. + */ + using GridThetaType = detail::InternalSpoofGridTheta; + /** + * @brief R grid. + */ + using GridRType = detail::InternalSpoofGridR; + /** + * @brief Phi grid. + */ + using GridPhiType = detail::InternalSpoofGridPhi; + /** + * @brief Sp grid. + */ + using GridSpType = GridSp; + + /** + * @brief Mu index range. + */ + using IdxRangeMuType = IdxRange<GridMuType>; + /** + * @brief Vpar index range. + */ + using IdxRangeVparType = IdxRange<GridVparType>; + /** + * @brief Sp,Vpar index range. + */ + using IdxRangeSpVparType = IdxRange<GridSpType, GridVparType>; + /** + * @brief Distribution function index range. + */ + using IdxRangeDistributionFunctionType = IdxRange<GridSpType, GridVparType, GridMuType>; + + /** + * @brief Container for the operator input data. + */ + using CollisionConfigurationDataType = detail::CollisionConfigurationData< + IdxRangeDistributionFunctionType, + GridSpType, + GridPhiType, + GridThetaType, + GridRType, + GridVparType, + GridMuType>; + +public: + /** + * @brief The constructor for the CollisionConfiguration class. + * + * @param yaml_input_file yaml namelist file object. + * @param index_range_fdistribution index range of the whole distribution + * function. + * @param coeff_intdmu quadrature coefficients. + * @param coeff_intdvpar quadrature coefficients. + * @param[in] B_norm The norm of the magnetic field defined over the polar + * slice. + * @param[in] Bstar_s Bstar defined on the coordinates + * (species,r,theta,vpar). + */ + CollisionConfiguration( + PC_tree_t const& yaml_input_file, + IdxRangeDistributionFunctionType index_range_fdistribution, + DConstField<IdxRangeMuType> coeff_intdmu, + DConstField<IdxRangeVparType> coeff_intdvpar, + double B_norm, + DConstField<IdxRangeSpVparType> Bstar_s) + : m_operator_quantities {yaml_input_file, index_range_fdistribution, coeff_intdmu, coeff_intdvpar} + , m_nustar0 {PCpp_double(yaml_input_file, ".Collisions.nustar0_rpeak")} + { + /* nustar fixed to 1. => normalisation of time is equal to 1./nustar + */ + if (m_nustar0 != 1.0) { + throw std::invalid_argument("nustar0 must be equal to 1"); + } + + do_configuration_data_initialisation(B_norm, Bstar_s); + } + + /** + * @brief Can be used to obtain the quantities the operator needs for its + * initialization, modification is not authorized. + * TODO(Etienne M): AFAIK, the DFieldMem behave like pointers, we may need + * accessors to enforce constness propagation. + * + * @return const CollisionConfigurationDataType& + */ + const CollisionConfigurationDataType& configuration() const + { + return m_operator_quantities; + } + +protected: + /** + * @brief Does the initialization specific to this geometry. + * + * @param[in] input_B_norm + * @param[in] input_Bstar_s + */ + void do_configuration_data_initialisation( + double input_B_norm, + DConstField<IdxRangeSpVparType> input_Bstar_s) + { + { + /* input_B_norm is a scalar when there is no Theta nor R dimensions. + */ + ddc::parallel_fill( + Kokkos::DefaultExecutionSpace(), + get_field(m_operator_quantities.m_B_norm_allocation), + input_B_norm); + } + + { + /* Koliop's Bstar_s is needed in a format layout different from the + * rest of the code. It is needed in layout right + * [sp, theta, r, vpar] where as the input is in [sp, vpar]. + * We simply iterate over [sp, theta, r, vpar] and copy the input + * [sp, vpar] for every theta, r. + */ + CollisionConfigurationDataType::transpose_replicate_deep_copy( + get_field(m_operator_quantities.m_Bstar_s_allocation), + input_Bstar_s); + } + + { + /* Computation of the radial profile of ADi that is required for + * koliop interface + * @f[ AD(rpeak) = \sqrt(2)*eps_rpeak^(3/2)/(q_rpeak R0)*nustar0_rpeak @f] + * eps_rpeak = eps(rpeak) with eps(r) = r/R0 is the aspect ratio. + * + * In this specific case, where AD is a scalar all values are forced + * to 1, so coeff_AD_atrpeak = sqrt(2) + */ + double const coeff_AD = Kokkos::sqrt(2.0); + + /* [TODO] According to Peter D., the places where coeff_AD is needed + * can resynthesise the quantities from a scalar and some other + * fields that koliop is already aware of. This change has not + * landed in Gysela Fortran yet. For koliop to make do with a scalar + * and not a radial profile, it would require koliop to modify its + * interface and the Cpar koliop operator computation to be modified + * slightly. For now, the uniform coeff_AD_rpeak is applied to + * produce a koliop compatible radial profile. + */ + ddc::parallel_fill( + Kokkos::DefaultExecutionSpace(), + get_field(m_operator_quantities.m_coeff_AD_allocation), + coeff_AD); + } + + Kokkos::fence(); + } + + /** + * @brief Container for the arrays the quantities the operator needs. + */ + CollisionConfigurationDataType m_operator_quantities; + + /** + * @brief Value of nustar0 = to nustar0_rpeak read in the YAML input file. + */ + double m_nustar0; +}; diff --git a/src/geometryVparMu/geometry/geometry.hpp b/src/geometryVparMu/geometry/geometry.hpp index a01b475b5..799ea7eb0 100644 --- a/src/geometryVparMu/geometry/geometry.hpp +++ b/src/geometryVparMu/geometry/geometry.hpp @@ -127,7 +127,7 @@ using IdxStepMu = IdxStep<GridMu>; using IdxStepVparMu = IdxStep<GridVpar, GridMu>; using IdxStepSpVparMu = IdxStep<Species, GridVpar, GridMu>; -// IdxRange = to describe the wole index range (or a sub-index range) +// IdxRange = to describe the whole index range (or a sub-index range) using IdxRangeVpar = IdxRange<GridVpar>; using IdxRangeMu = IdxRange<GridMu>; using IdxRangeVparMu = IdxRange<GridVpar, GridMu>; diff --git a/toolchains/mi250.hipcc.adastra.spack/prepare.sh b/toolchains/mi250.hipcc.adastra.spack/prepare.sh index 788a31af8..a03ed19ee 100755 --- a/toolchains/mi250.hipcc.adastra.spack/prepare.sh +++ b/toolchains/mi250.hipcc.adastra.spack/prepare.sh @@ -48,6 +48,8 @@ eigen%gcc~ipo build_system=cmake build_type=Release generator=ninja arch=linux-r " # openblas@0.3.26%gcc@12.1.generic~bignuma~consistent_fpcsr+dynamic_dispatch+fortran~ilp64+locking+pic+shared build_system=makefile symbol_suffix=none threads=none arch=linux-rhel8-zen3 +which spack + # If we start preparing a new environment, ensure we wont get name clashes by # uninstalling previous products. # NOTE: This may fail if the product are not already installed. IMO this is a diff --git a/vendor/eigen b/vendor/eigen index ad13df7ea..2e76277bd 160000 --- a/vendor/eigen +++ b/vendor/eigen @@ -1 +1 @@ -Subproject commit ad13df7ea444f9b78f475d4dd6fc754953676430 +Subproject commit 2e76277bd049f7bec36b0f908c69734a42c5234f diff --git a/vendor/kokkos b/vendor/kokkos index 18cb82e57..4764a9a79 160000 --- a/vendor/kokkos +++ b/vendor/kokkos @@ -1 +1 @@ -Subproject commit 18cb82e5764a68d6f0633d81a50982f3804692e0 +Subproject commit 4764a9a79644b707a07a71ea712f84885f76beed