Clean geometryXYVxVY simulations

Closes #251

See merge request gysela-developpers/gyselalibxx!530

--------------------------------------------

simulations/geometryXYVxVy/landau/landau4d_fft.cpp

@@ -85,34 +85,14 @@ int main(int argc, char** argv)
     SplineVxBuilder_1d const builder_vx_1d(mesh_vx);
     SplineVyBuilder_1d const builder_vy_1d(mesh_vy);
 
-    host_t<DFieldSp> density_eq(dom_kinsp);
-    host_t<DFieldSp> temperature_eq(dom_kinsp);
-    host_t<DFieldSp> mean_velocity_eq(dom_kinsp);
-    host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
-    host_t<FieldSp<int>> init_perturb_mode(dom_kinsp);
-
-    for (IndexSp const isp : dom_kinsp) {
-        PC_tree_t const conf_isp = PCpp_get(conf_voicexx, ".SpeciesInfo[%d]", isp.uid());
-
-        density_eq(isp) = PCpp_double(conf_isp, ".density_eq");
-        temperature_eq(isp) = PCpp_double(conf_isp, ".temperature_eq");
-        mean_velocity_eq(isp) = PCpp_double(conf_isp, ".mean_velocity_eq");
-        init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
-        init_perturb_mode(isp) = PCpp_double(conf_isp, ".perturb_mode");
-    }
-
     // Initialization of the distribution function
     DFieldSpVxVy allfequilibrium(meshSpVxVy);
-    MaxwellianEquilibrium const init_fequilibrium(
-            std::move(density_eq),
-            std::move(temperature_eq),
-            std::move(mean_velocity_eq));
+    MaxwellianEquilibrium const init_fequilibrium
+            = MaxwellianEquilibrium::init_from_input(dom_kinsp, conf_voicexx);
     init_fequilibrium(allfequilibrium);
     DFieldSpXYVxVy allfdistribu(meshSpXYVxVy);
-    SingleModePerturbInitialization const
-            init(allfequilibrium,
-                 init_perturb_mode.span_cview(),
-                 init_perturb_amplitude.span_cview());
+    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());
 

src/geometryXVx/initialization/singlemodeperturbinitialization.cpp

@@ -25,10 +25,10 @@ DSpanSpXVx SingleModePerturbInitialization::operator()(DSpanSpXVx const allfdist
     // Initialization of the perturbation
     DFieldX perturbation_alloc(gridx);
     ddc::ChunkSpan fequilibrium_proxy = m_fequilibrium.span_view();
-    ddc::ChunkSpan perturbation = perturbation_alloc.span_view();
+    ddc::ChunkSpan perturbation_proxy = perturbation_alloc.span_view();
     ddc::for_each(gridsp, [&](IndexSp const isp) {
         perturbation_initialization(
-                perturbation,
+                perturbation_proxy,
                 m_init_perturb_mode(isp),
                 m_init_perturb_amplitude(isp));
         // Initialization of the distribution function --> fill values
@@ -38,7 +38,8 @@ DSpanSpXVx SingleModePerturbInitialization::operator()(DSpanSpXVx const allfdist
                 KOKKOS_LAMBDA(IndexXVx const ixvx) {
                     IndexX const ix = ddc::select<IDimX>(ixvx);
                     IndexVx const ivx = ddc::select<IDimVx>(ixvx);
-                    double fdistribu_val = fequilibrium_proxy(isp, ivx) * (1. + perturbation(ix));
+                    double fdistribu_val
+                            = fequilibrium_proxy(isp, ivx) * (1. + perturbation_proxy(ix));
                     if (fdistribu_val < 1.e-60) {
                         fdistribu_val = 1.e-60;
                     }

src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp

@@ -39,6 +39,30 @@ DSpanSpVxVy MaxwellianEquilibrium::operator()(DSpanSpVxVy const allfequilibrium)
     return allfequilibrium;
 }
 
+
+MaxwellianEquilibrium MaxwellianEquilibrium::init_from_input(
+        IDomainSp dom_kinsp,
+        PC_tree_t const& yaml_input_file)
+{
+    host_t<DFieldSp> density_eq(dom_kinsp);
+    host_t<DFieldSp> temperature_eq(dom_kinsp);
+    host_t<DFieldSp> mean_velocity_eq(dom_kinsp);
+
+    for (IndexSp const isp : dom_kinsp) {
+        PC_tree_t const conf_isp = PCpp_get(yaml_input_file, ".SpeciesInfo[%d]", isp.uid());
+
+        density_eq(isp) = PCpp_double(conf_isp, ".density_eq");
+        temperature_eq(isp) = PCpp_double(conf_isp, ".temperature_eq");
+        mean_velocity_eq(isp) = PCpp_double(conf_isp, ".mean_velocity_eq");
+    }
+
+    return MaxwellianEquilibrium(
+            std::move(density_eq),
+            std::move(temperature_eq),
+            std::move(mean_velocity_eq));
+}
+
+
 /*
  Computing the Maxwellian function as
   fM(vx,vy) = n/(2*PI*T)*exp(-(vx**2+vy**2)/(2*T))

src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp

@@ -2,10 +2,12 @@
 
 #pragma once
 
+#include <paraconf.h>
 #include <species_info.hpp>
 
 #include "geometry.hpp"
 #include "iequilibrium.hpp"
+#include "paraconfpp.hpp"
 
 /// Equilibrium operator as Maxwellian. This initializes all species.
 class MaxwellianEquilibrium : public IEquilibrium
@@ -33,6 +35,16 @@ class MaxwellianEquilibrium : public IEquilibrium
 
     ~MaxwellianEquilibrium() override = default;
 
+    /**
+     * @brief Read the density, temperature and mean velocity required to initialize the Maxwellian in a YAML input file.
+     * @param[in] dom_kinsp Discrete Domain for the kinetic species
+     * @param[in] yaml_input_file YAML input file
+     * @return an instance of Maxwellian distribution function.
+     */
+    static MaxwellianEquilibrium init_from_input(
+            IDomainSp dom_kinsp,
+            PC_tree_t const& yaml_input_file);
+
     /**
      * @brief Initializes allfequilibrium as a Maxwellian.
      * @param[out] allfequilibrium A Span containing a Maxwellian distribution function.

src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp

@@ -6,14 +6,13 @@
 #include "singlemodeperturbinitialization.hpp"
 
 
-
 SingleModePerturbInitialization::SingleModePerturbInitialization(
         DViewSpVxVy fequilibrium,
-        host_t<ViewSp<int>> const init_perturb_mode,
-        host_t<DViewSp> const init_perturb_amplitude)
+        host_t<IFieldSp> init_perturb_mode,
+        host_t<DFieldSp> init_perturb_amplitude)
     : m_fequilibrium(fequilibrium)
-    , m_init_perturb_mode(init_perturb_mode)
-    , m_init_perturb_amplitude(init_perturb_amplitude)
+    , m_init_perturb_mode(std::move(init_perturb_mode))
+    , m_init_perturb_amplitude(std::move(init_perturb_amplitude))
 {
 }
 
@@ -25,15 +24,15 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al
 
     // Initialization of the perturbation
     DFieldXY perturbation_alloc(gridxy);
-    DSpanXY perturbation = perturbation_alloc.span_view();
+    ddc::ChunkSpan fequilibrium_proxy = m_fequilibrium.span_view();
+    ddc::ChunkSpan perturbation_proxy = perturbation_alloc.span_view();
     ddc::for_each(gridsp, [&](IndexSp const isp) {
         perturbation_initialization(
-                perturbation,
+                perturbation_proxy,
                 m_init_perturb_mode(isp),
                 m_init_perturb_amplitude(isp));
 
         // Initialization of the distribution function --> fill values
-        DViewSpVxVy fequilibrium_proxy = m_fequilibrium;
         ddc::parallel_for_each(
                 Kokkos::DefaultExecutionSpace(),
                 gridxyvxvy,
@@ -43,7 +42,7 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al
                     IndexVx const ivx = ddc::select<IDimVx>(ixyvxvy);
                     IndexVy const ivy = ddc::select<IDimVy>(ixyvxvy);
                     double fdistribu_val
-                            = fequilibrium_proxy(isp, ivx, ivy) * (1. + perturbation(ix, iy));
+                            = fequilibrium_proxy(isp, ivx, ivy) * (1. + perturbation_proxy(ix, iy));
                     if (fdistribu_val < 1.e-60) {
                         fdistribu_val = 1.e-60;
                     }
@@ -53,16 +52,39 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al
     return allfdistribu;
 }
 
+
+SingleModePerturbInitialization SingleModePerturbInitialization::init_from_input(
+        DViewSpVxVy allfequilibrium,
+        IDomainSp dom_kinsp,
+        PC_tree_t const& yaml_input_file)
+{
+    host_t<IFieldSp> init_perturb_mode(dom_kinsp);
+    host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
+
+    for (IndexSp const isp : dom_kinsp) {
+        PC_tree_t const conf_isp = PCpp_get(yaml_input_file, ".SpeciesInfo[%d]", isp.uid());
+
+        init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
+        init_perturb_mode(isp) = static_cast<int>(PCpp_int(conf_isp, ".perturb_mode"));
+    }
+
+    return SingleModePerturbInitialization(
+            allfequilibrium,
+            std::move(init_perturb_mode),
+            std::move(init_perturb_amplitude));
+}
+
+
 void SingleModePerturbInitialization::perturbation_initialization(
         DSpanXY const perturbation,
-        int const mode,
+        int const perturb_mode,
         double const perturb_amplitude) const
 {
     IDomainXY const gridxy = perturbation.domain<IDimX, IDimY>();
-    double const kx
-            = mode * 2. * M_PI / ddcHelper::total_interval_length(ddc::select<IDimX>(gridxy));
-    double const ky
-            = mode * 2. * M_PI / ddcHelper::total_interval_length(ddc::select<IDimY>(gridxy));
+    double const kx = perturb_mode * 2. * M_PI
+                      / ddcHelper::total_interval_length(ddc::select<IDimX>(gridxy));
+    double const ky = perturb_mode * 2. * M_PI
+                      / ddcHelper::total_interval_length(ddc::select<IDimY>(gridxy));
 
     ddc::parallel_for_each(
             Kokkos::DefaultExecutionSpace(),

src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp

@@ -3,30 +3,34 @@
 #pragma once
 
 #include <geometry.hpp>
+#include <paraconf.h>
 #include <species_info.hpp>
 
 #include "iinitialization.hpp"
+#include "paraconfpp.hpp"
 
 /// Initialization operator with a sinusoidal perturbation of a Maxwellian. This initializes all species.
 class SingleModePerturbInitialization : public IInitialization
 {
     DViewSpVxVy m_fequilibrium;
 
-    host_t<ViewSp<int>> m_init_perturb_mode;
+    host_t<IFieldSp> m_init_perturb_mode;
 
-    host_t<DViewSp> m_init_perturb_amplitude;
+    host_t<DFieldSp> m_init_perturb_amplitude;
 
 public:
     /**
      * @brief Initialization of the perturbation.
      * @param[in, out] perturbation On input: an uninitialized array
      *                              On output: an array containing a values that has a 
      *                              sinusoidal variation with given amplitude and mode. 
-     * @param[in] mode The mode of the perturbation. 
+     * @param[in] perturb_mode The mode of the perturbation. 
      * @param[in] perturb_amplitude The amplitude of the perturbation. 
      */
-    void perturbation_initialization(DSpanXY perturbation, int mode, double perturb_amplitude)
-            const;
+    void perturbation_initialization(
+            DSpanXY perturbation,
+            int const perturb_mode,
+            double const perturb_amplitude) const;
 
     /**
      * @brief Creates an instance of the SingleModePerturbInitialization class.
@@ -36,8 +40,8 @@ class SingleModePerturbInitialization : public IInitialization
      */
     SingleModePerturbInitialization(
             DViewSpVxVy fequilibrium,
-            host_t<ViewSp<int>> init_perturb_mode,
-            host_t<DViewSp> init_perturb_amplitude);
+            host_t<IFieldSp> init_perturb_mode,
+            host_t<DFieldSp> init_perturb_amplitude);
 
     ~SingleModePerturbInitialization() override = default;
 
@@ -47,4 +51,17 @@ class SingleModePerturbInitialization : public IInitialization
      * @return The initialized distribution function.
      */
     DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu) const override;
-};
+
+    /**
+     * @brief Read init_perturb_mode and init_perturb amplitude in a YAML input file 
+     *      to initialize the perturbation. 
+     * @param[in] allfequilibrium equilibrium distribution function.
+     * @param[in] dom_kinsp Discrete Domain for the kinetic species.
+     * @param[in] yaml_input_file YAML input file.
+     * @return an instance of SingleModePerturbInitialization class.
+     */
+    static SingleModePerturbInitialization init_from_input(
+            DViewSpVxVy allfequilibrium,
+            IDomainSp dom_kinsp,
+            PC_tree_t const& yaml_input_file);
+};

tests/geometryXYVxVy/landau/landau.yaml

@@ -19,7 +19,7 @@ SpeciesInfo:
   temperature_eq: 1.
   mean_velocity_eq: 0.
   perturb_amplitude: 0.05
-  perturb_mode: 0.5
+  perturb_mode: 1
 
 Algorithm:
   deltat: 0.0625

tests/geometryXYVxVy/landau/landau_small.yaml

@@ -19,7 +19,7 @@ SpeciesInfo:
   temperature_eq: 1.
   mean_velocity_eq: 0.
   perturb_amplitude: 0.05
-  perturb_mode: 0.5
+  perturb_mode: 1
 
 Algorithm:
   deltat: 0.0625