Use Kokkos::Parallel in SourceTerm Setup

include/GMGPolar/gmgpolar.h

@@ -139,10 +139,14 @@ class GMGPolar : public IGMGPolar
     /* --------------- */
     /* Setup Functions */
     int chooseNumberOfLevels(const PolarGrid& finest_grid);
+    // Public due to cuda restrictions
+public:
     template <concepts::BoundaryConditions BoundaryConditions, concepts::SourceTerm SourceTerm>
     void build_rhs_f(const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f,
                      const BoundaryConditions& boundary_conditions, const SourceTerm& source_term);
     void discretize_rhs_f(const Level<DomainGeometry, DensityProfileCoefficients>& level, Vector<double> rhs_f);
+
+private:
     bool checkUniformRefinement(const PolarGrid& grid, double tolerance) const;
 
     /* --------------- */

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CircularGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,11 @@ class CartesianR2_Boundary_CircularGeometry
 {
 public:
     explicit CartesianR2_Boundary_CircularGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_Boundary_CircularGeometry(const CartesianR2_Boundary_CircularGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_CzarnyGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -13,9 +14,10 @@ class CartesianR2_Boundary_CzarnyGeometry
     explicit CartesianR2_Boundary_CzarnyGeometry();
     explicit CartesianR2_Boundary_CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon,
                                                  double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION CartesianR2_Boundary_CzarnyGeometry(const CartesianR2_Boundary_CzarnyGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax                         = 1.3;

include/InputFunctions/BoundaryConditions/cartesianR2_Boundary_ShafranovGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,11 @@ class CartesianR2_Boundary_ShafranovGeometry
 {
 public:
     explicit CartesianR2_Boundary_ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_Boundary_ShafranovGeometry(const CartesianR2_Boundary_ShafranovGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax             = 1.3;

include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_CircularGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,11 @@ class CartesianR6_Boundary_CircularGeometry
 {
 public:
     explicit CartesianR6_Boundary_CircularGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR6_Boundary_CircularGeometry(const CartesianR6_Boundary_CircularGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_CzarnyGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -13,9 +14,10 @@ class CartesianR6_Boundary_CzarnyGeometry
     explicit CartesianR6_Boundary_CzarnyGeometry();
     explicit CartesianR6_Boundary_CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon,
                                                  double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION CartesianR6_Boundary_CzarnyGeometry(const CartesianR6_Boundary_CzarnyGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax                         = 1.3;

include/InputFunctions/BoundaryConditions/cartesianR6_Boundary_ShafranovGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,11 @@ class CartesianR6_Boundary_ShafranovGeometry
 {
 public:
     explicit CartesianR6_Boundary_ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR6_Boundary_ShafranovGeometry(const CartesianR6_Boundary_ShafranovGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax             = 1.3;

include/InputFunctions/BoundaryConditions/polarR6_Boundary_CircularGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class PolarR6_Boundary_CircularGeometry
 {
 public:
     explicit PolarR6_Boundary_CircularGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION PolarR6_Boundary_CircularGeometry(const PolarR6_Boundary_CircularGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/polarR6_Boundary_CulhamGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class PolarR6_Boundary_CulhamGeometry
 {
 public:
     explicit PolarR6_Boundary_CulhamGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION PolarR6_Boundary_CulhamGeometry(const PolarR6_Boundary_CulhamGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/polarR6_Boundary_CzarnyGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -12,9 +13,10 @@ class PolarR6_Boundary_CzarnyGeometry
 public:
     explicit PolarR6_Boundary_CzarnyGeometry();
     explicit PolarR6_Boundary_CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon, double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION PolarR6_Boundary_CzarnyGeometry(const PolarR6_Boundary_CzarnyGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax                         = 1.3;

include/InputFunctions/BoundaryConditions/polarR6_Boundary_ShafranovGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class PolarR6_Boundary_ShafranovGeometry
 {
 public:
     explicit PolarR6_Boundary_ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta);
+    KOKKOS_DEFAULTED_FUNCTION PolarR6_Boundary_ShafranovGeometry(const PolarR6_Boundary_ShafranovGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax             = 1.3;

include/InputFunctions/BoundaryConditions/refined_Boundary_CircularGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class Refined_Boundary_CircularGeometry
 {
 public:
     explicit Refined_Boundary_CircularGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION Refined_Boundary_CircularGeometry(const Refined_Boundary_CircularGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/refined_Boundary_CulhamGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class Refined_Boundary_CulhamGeometry
 {
 public:
     explicit Refined_Boundary_CulhamGeometry(double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION Refined_Boundary_CulhamGeometry(const Refined_Boundary_CulhamGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax = 1.3;

include/InputFunctions/BoundaryConditions/refined_Boundary_CzarnyGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -12,9 +13,10 @@ class Refined_Boundary_CzarnyGeometry
 public:
     explicit Refined_Boundary_CzarnyGeometry();
     explicit Refined_Boundary_CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon, double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION Refined_Boundary_CzarnyGeometry(const Refined_Boundary_CzarnyGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax                         = 1.3;

include/InputFunctions/BoundaryConditions/refined_Boundary_ShafranovGeometry.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include "../boundaryConditions.h"
 
@@ -11,9 +12,10 @@ class Refined_Boundary_ShafranovGeometry
 {
 public:
     explicit Refined_Boundary_ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta);
+    KOKKOS_DEFAULTED_FUNCTION Refined_Boundary_ShafranovGeometry(const Refined_Boundary_ShafranovGeometry&) = default;
 
-    double u_D(double r, double theta) const;
-    double u_D_Interior(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D(double r, double theta) const;
+    KOKKOS_FUNCTION double u_D_Interior(double r, double theta) const;
 
 private:
     const double Rmax             = 1.3;

include/InputFunctions/SourceTerms/cartesianR2_Poisson_CircularGeometry.h

@@ -13,11 +13,13 @@ class CartesianR2_Poisson_CircularGeometry
 {
 public:
     explicit CartesianR2_Poisson_CircularGeometry(PolarGrid const& grid, double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_Poisson_CircularGeometry(const CartesianR2_Poisson_CircularGeometry&) = default;
 
-    double operator()(std::size_t i_r, std::size_t i_theta) const;
+    KOKKOS_FUNCTION double operator()(std::size_t i_r, std::size_t i_theta) const;
 
 private:
-    PolarGrid const& grid_;
+    PolarGrid grid_;
     const double Rmax = 1.3;
 };
 } // namespace gmgpolar

include/InputFunctions/SourceTerms/cartesianR2_Poisson_CzarnyGeometry.h

@@ -14,11 +14,12 @@ class CartesianR2_Poisson_CzarnyGeometry
 public:
     explicit CartesianR2_Poisson_CzarnyGeometry(PolarGrid const& grid, double Rmax, double inverse_aspect_ratio_epsilon,
                                                 double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION CartesianR2_Poisson_CzarnyGeometry(const CartesianR2_Poisson_CzarnyGeometry&) = default;
 
-    double operator()(std::size_t i_r, std::size_t i_theta) const;
+    KOKKOS_FUNCTION double operator()(std::size_t i_r, std::size_t i_theta) const;
 
 private:
-    PolarGrid const& grid_;
+    PolarGrid grid_;
     const double Rmax                         = 1.3;
     const double inverse_aspect_ratio_epsilon = 0.3;
     const double ellipticity_e                = 1.4;

include/InputFunctions/SourceTerms/cartesianR2_Poisson_ShafranovGeometry.h

@@ -14,11 +14,13 @@ class CartesianR2_Poisson_ShafranovGeometry
 public:
     explicit CartesianR2_Poisson_ShafranovGeometry(PolarGrid const& grid, double Rmax, double elongation_kappa,
                                                    double shift_delta);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_Poisson_ShafranovGeometry(const CartesianR2_Poisson_ShafranovGeometry&) = default;
 
-    double operator()(std::size_t i_r, std::size_t i_theta) const;
+    KOKKOS_FUNCTION double operator()(std::size_t i_r, std::size_t i_theta) const;
 
 private:
-    PolarGrid const& grid_;
+    PolarGrid grid_;
     const double Rmax             = 1.3;
     const double elongation_kappa = 0.3;
     const double shift_delta      = 0.2;

include/InputFunctions/SourceTerms/cartesianR2_SonnendruckerGyro_CircularGeometry.h

@@ -13,11 +13,13 @@ class CartesianR2_SonnendruckerGyro_CircularGeometry
 {
 public:
     explicit CartesianR2_SonnendruckerGyro_CircularGeometry(PolarGrid const& grid, double Rmax);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_SonnendruckerGyro_CircularGeometry(const CartesianR2_SonnendruckerGyro_CircularGeometry&) = default;
 
-    double operator()(std::size_t i_r, std::size_t i_theta) const;
+    KOKKOS_FUNCTION double operator()(std::size_t i_r, std::size_t i_theta) const;
 
 private:
-    PolarGrid const& grid_;
+    PolarGrid grid_;
     const double Rmax = 1.3;
 };
 } // namespace gmgpolar

include/InputFunctions/SourceTerms/cartesianR2_SonnendruckerGyro_CzarnyGeometry.h

@@ -14,11 +14,13 @@ class CartesianR2_SonnendruckerGyro_CzarnyGeometry
 public:
     explicit CartesianR2_SonnendruckerGyro_CzarnyGeometry(PolarGrid const& grid, double Rmax,
                                                           double inverse_aspect_ratio_epsilon, double ellipticity_e);
+    KOKKOS_DEFAULTED_FUNCTION
+    CartesianR2_SonnendruckerGyro_CzarnyGeometry(const CartesianR2_SonnendruckerGyro_CzarnyGeometry&) = default;
 
-    double operator()(std::size_t i_r, std::size_t i_theta) const;
+    KOKKOS_FUNCTION double operator()(std::size_t i_r, std::size_t i_theta) const;
 
 private:
-    PolarGrid const& grid_;
+    PolarGrid grid_;
     const double Rmax                         = 1.3;
     const double inverse_aspect_ratio_epsilon = 0.3;
     const double ellipticity_e                = 1.4;