From 94ff94ac51731718cd41528cb41ac5fb3415d204 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 11 May 2026 23:16:31 +0200 Subject: [PATCH 1/8] Residual skibidi skibidi --- include/Residual/ResidualGive/applyAGive.inl | 260 +++++++----------- include/Residual/ResidualGive/residualGive.h | 4 +- .../Residual/ResidualGive/residualGive.inl | 112 +------- include/Residual/ResidualTake/applyATake.inl | 9 +- include/Residual/ResidualTake/residualTake.h | 6 +- .../Residual/ResidualTake/residualTake.inl | 10 +- include/Residual/residual.h | 3 +- 7 files changed, 121 insertions(+), 283 deletions(-) diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 2802e0ef..a97fad28 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -4,30 +4,33 @@ namespace residual_give { template -static inline void node_apply_a_give(int i_r, int i_theta, double r, double theta, const PolarGrid& grid, - const LevelCacheType& level_cache, bool DirBC_Interior, Vector& result, - ConstVector& x) +static inline void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid, const LevelCacheType& level_cache, + bool DirBC_Interior, Vector& result, ConstVector& x) { /* ---------------------------------------- */ /* Compute or retrieve stencil coefficients */ /* ---------------------------------------- */ - const int center = grid.index(i_r, i_theta); + const int center = grid.index(i_r, i_theta); + const double r = grid.radius(i_r); + const double theta = grid.theta(i_theta); + double coeff_beta, arr, att, art, detDF; level_cache.obtainValues(i_r, i_theta, center, r, theta, coeff_beta, arr, att, art, detDF); /* -------------------- */ /* Node in the interior */ /* -------------------- */ - if (i_r > 1 && i_r < grid.nr() - 2) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; + if (i_r > 0 && i_r < grid.nr() - 1) { + const double h1 = grid.radialSpacing(i_r - 1); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); + + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -38,23 +41,28 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ - coeff4 * att * x[top] /* Top */ + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[center]) /* Center: (Left, Right, Bottom, Top) */; + x[center]); /* Center: (Left, Right, Bottom, Top) */ + /* Fill result(i-1,j) */ - result[left] += (-coeff1 * arr * x[center] /* Right */ - + coeff1 * arr * x[left] /* Center: (Right) */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ + if (i_r > 1 || (i_r == 1 && !DirBC_Interior)) { /* Don't give to the inner dirichlet boundary! */ + result[left] += (-coeff1 * arr * x[center] /* Right */ + + coeff1 * arr * x[left] /* Center: (Right) */ + - 0.25 * art * x[top] /* Top Right */ + + 0.25 * art * x[bottom]); /* Bottom Right */ + } /* Fill result(i+1,j) */ - result[right] += (-coeff2 * arr * x[center] /* Left */ - + coeff2 * arr * x[right] /* Center: (Left) */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ + if (i_r < grid.nr() - 2) { /* Don't give to the outer dirichlet boundary! */ + result[right] += (-coeff2 * arr * x[center] /* Left */ + + coeff2 * arr * x[right] /* Center: (Left) */ + + 0.25 * art * x[top] /* Top Left */ + - 0.25 * art * x[bottom]); /* Bottom Left */ + } /* Fill result(i,j-1) */ result[bottom] += (-coeff3 * att * x[center] /* Top */ + coeff3 * att * x[bottom] /* Center: (Top) */ @@ -104,15 +112,16 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet /* h1 gets replaced with 2 * R0. */ /* (i_r-1,i_theta) gets replaced with (i_r, i_theta + (grid.ntheta()/2)). */ /* Some more adjustments from the changing the 9-point stencil to the artifical 7-point stencil. */ - double h1 = 2.0 * grid.radius(0); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + const double h1 = 2.0 * grid.radius(0); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff3 = 0.5 * (h1 + h2) / k1; + const double coeff4 = 0.5 * (h1 + h2) / k2; + const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -124,7 +133,7 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet const int top = grid.index(i_r, i_theta_P1); /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ + result[center] += (coeff5 * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - coeff1 * arr * x[left] /* Left */ - coeff2 * arr * x[right] /* Right */ - coeff3 * att * x[bottom] /* Bottom */ @@ -154,112 +163,6 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet /* - 0.25 * art * x[left]; // Bottom Left: REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */ } } - /* ------------------------------- */ - /* Node next to the inner boundary */ - /* ------------------------------- */ - else if (i_r == 1) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; - - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - const int left = grid.index(i_r - 1, i_theta); - const int right = grid.index(i_r + 1, i_theta); - const int bottom = grid.index(i_r, i_theta_M1); - const int top = grid.index(i_r, i_theta_P1); - - /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - - coeff1 * arr * x[left] /* Left */ - - coeff2 * arr * x[right] /* Right */ - - coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top] /* Top */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[center]); /* Center: (Left, Right, Bottom, Top) */ - /* Fill result(i-1,j) */ - if (!DirBC_Interior) { /* Don't give to the inner dirichlet boundary! */ - result[left] += (-coeff1 * arr * x[center] /* Right */ - + coeff1 * arr * x[left] /* Center: (Right) */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ - } - /* Fill result(i+1,j) */ - result[right] += (-coeff2 * arr * x[center] /* Left */ - + coeff2 * arr * x[right] /* Center: (Left) */ - + 0.25 * art * x[top] /* Top Left */ - - 0.25 * art * x[bottom]); /* Bottom Left */ - /* Fill result(i,j-1) */ - result[bottom] += (-coeff3 * att * x[center] /* Top */ - + coeff3 * att * x[bottom] /* Center: (Top) */ - - 0.25 * art * x[right] /* Top Right */ - + 0.25 * art * x[left]); /* Top Left */ - /* Fill result(i,j+1) */ - result[top] += (-coeff4 * att * x[center] /* Bottom */ - + coeff4 * att * x[top] /* Center: (Bottom) */ - + 0.25 * art * x[right] /* Bottom Right */ - - 0.25 * art * x[left]); /* Bottom Left */ - } - /* ------------------------------- */ - /* Node next to the outer boundary */ - /* ------------------------------- */ - else if (i_r == grid.nr() - 2) { - double h1 = grid.radialSpacing(i_r - 1); - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); - - double coeff1 = 0.5 * (k1 + k2) / h1; - double coeff2 = 0.5 * (k1 + k2) / h2; - double coeff3 = 0.5 * (h1 + h2) / k1; - double coeff4 = 0.5 * (h1 + h2) / k2; - - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); - - const int left = grid.index(i_r - 1, i_theta); - const int right = grid.index(i_r + 1, i_theta); - const int bottom = grid.index(i_r, i_theta_M1); - const int top = grid.index(i_r, i_theta_P1); - - /* Fill result(i,j) */ - result[center] += (0.25 * (h1 + h2) * (k1 + k2) * coeff_beta * std::fabs(detDF) * x[center] /* beta_{i,j} */ - - coeff1 * arr * x[left] /* Left */ - - coeff2 * arr * x[right] /* Right */ - - coeff3 * att * x[bottom] /* Bottom */ - - coeff4 * att * x[top] /* Top */ - + ((coeff1 + coeff2) * arr + (coeff3 + coeff4) * att) * - x[center]); /* Center: (Left, Right, Bottom, Top) */ - /* Fill result(i-1,j) */ - result[left] += (-coeff1 * arr * x[center] /* Right */ - + coeff1 * arr * x[left] /* Center: (Right) */ - - 0.25 * art * x[top] /* Top Right */ - + 0.25 * art * x[bottom]); /* Bottom Right */ - /* Don't give to the outer dirichlet boundary! */ - /* Fill result(i+1,j) */ - /* result[right] += ( */ - /* - coeff2 * arr * x[center] // Left */ - /* + coeff2 * arr * x[right] // Center: (Left) */ - /* + 0.25 * art * x[top] // Top Left */ - /* - 0.25 * art * x[bottom]); // Bottom Left */ - /* Fill result(i,j-1) */ - result[bottom] += (-coeff3 * att * x[center] /* Top */ - + coeff3 * att * x[bottom] /* Center: (Top) */ - - 0.25 * art * x[right] /* Top Right */ - + 0.25 * art * x[left]); /* Top Left */ - /* Fill result(i,j+1) */ - result[top] += (-coeff4 * att * x[center] /* Bottom */ - + coeff4 * att * x[top] /* Center: (Bottom) */ - + 0.25 * art * x[right] /* Bottom Right */ - - 0.25 * art * x[left]); /* Bottom Left */ - } /* ----------------------------- */ /* Node on to the outer boundary */ /* ----------------------------- */ @@ -292,32 +195,69 @@ static inline void node_apply_a_give(int i_r, int i_theta, double r, double thet } // namespace residual_give template -void ResidualGive::applyCircleSection(const int i_r, Vector result, ConstVector x) const +void ResidualGive::applySystemOperator(Vector result, ConstVector x) const { - using residual_give::node_apply_a_give; - - const PolarGrid& grid = Residual::grid_; + assert(result.size() == x.size()); - const double r = grid.radius(i_r); - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - const double theta = grid.theta(i_theta); - node_apply_a_give(i_r, i_theta, r, theta, grid, Residual::level_cache_, - Residual::DirBC_Interior_, result, x); - } -} + const PolarGrid& grid = Residual::grid_; + const LevelCacheType& level_cache = Residual::level_cache_; + const bool DirBC_Interior = Residual::DirBC_Interior_; -template -void ResidualGive::applyRadialSection(const int i_theta, Vector result, - ConstVector x) const -{ using residual_give::node_apply_a_give; - const PolarGrid& grid = Residual::grid_; + assign(result, 0.0); + + /* ---------------- */ + /* Circular section */ + /* ---------------- */ + // We parallelize over i_r (step 3) to avoid data race conditions between adjacent circles. + // The i_theta loop is sequential inside the kernel. + const int num_circle_tasks = grid.numberSmootherCircles(); + + for (int start_circle = 0; start_circle < 3; ++start_circle) { + const int num_circular_tasks = (num_circle_tasks - start_circle + 2) / 3; + Kokkos::parallel_for( + "ResidualGive: ApplyA (Circular)", Kokkos::RangePolicy<>(0, num_circular_tasks), + KOKKOS_LAMBDA(const int circle_task) { + const int i_r = start_circle + circle_task * 3; + for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { + node_apply_a_give(i_r, i_theta, grid, level_cache, DirBC_Interior, result, x); + } + }); + Kokkos::fence(); + } - const double theta = grid.theta(i_theta); - for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { - const double r = grid.radius(i_r); - node_apply_a_give(i_r, i_theta, r, theta, grid, Residual::level_cache_, - Residual::DirBC_Interior_, result, x); + /* -------------- */ + /* Radial section */ + /* -------------- */ + // We parallelize over i_theta (step 3) to avoid data race conditions between adjacent radial lines. + // The i_r loop is sequential inside the kernel. + // Due to periodicity in the angular direction, handle up to 2 additional + // radial lines (i_theta = 0 and 1) before the parallel passes. + const int additional_radial_tasks = grid.ntheta() % 3; + const int num_radial_tasks = grid.ntheta() - additional_radial_tasks; + + for (int i_theta = 0; i_theta < additional_radial_tasks; i_theta++) { + Kokkos::parallel_for( + "ResidualGive: ApplyA (Radial, additional)", Kokkos::RangePolicy<>(0, 1), KOKKOS_LAMBDA(const int) { + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + node_apply_a_give(i_r, i_theta, grid, level_cache, DirBC_Interior, result, x); + } + }); + Kokkos::fence(); + } + + for (int start_radial = 0; start_radial < 3; ++start_radial) { + const int num_radial_batches = (num_radial_tasks - start_radial + 2) / 3; + Kokkos::parallel_for( + "ResidualGive: ApplyA (Radial)", Kokkos::RangePolicy<>(0, num_radial_batches), + KOKKOS_LAMBDA(const int radial_task) { + const int i_theta = start_radial + radial_task * 3; + for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { + node_apply_a_give(i_r, i_theta, grid, level_cache, DirBC_Interior, result, x); + } + }); + Kokkos::fence(); } } +// clang-format on \ No newline at end of file diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h index 37d1040d..4bbe254b 100644 --- a/include/Residual/ResidualGive/residualGive.h +++ b/include/Residual/ResidualGive/residualGive.h @@ -13,9 +13,8 @@ class ResidualGive : public Residual const int num_omp_threads); ~ResidualGive() override = default; - void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; - void applySystemOperator(Vector result, ConstVector x) const final; + void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; private: void applyCircleSection(const int i_r, Vector result, ConstVector x) const; @@ -24,4 +23,5 @@ class ResidualGive : public Residual #include "residualGive.inl" #include "applyAGive.inl" + } // namespace gmgpolar diff --git a/include/Residual/ResidualGive/residualGive.inl b/include/Residual/ResidualGive/residualGive.inl index bc50f620..122efdf1 100644 --- a/include/Residual/ResidualGive/residualGive.inl +++ b/include/Residual/ResidualGive/residualGive.inl @@ -7,105 +7,6 @@ ResidualGive::ResidualGive(const PolarGrid& grid, const LevelCac { } -/* ------------ */ -/* result = A*x */ - -// clang-format off -template -void ResidualGive::applySystemOperator(Vector result, ConstVector x) const -{ - assert(result.size() == x.size()); - - assign(result, 0.0); - - const PolarGrid& grid = Residual::grid_; - - const int num_omp_threads = Residual::num_omp_threads_; - - /* Single-threaded execution */ - if (num_omp_threads == 1) { - for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) { - applyCircleSection(i_r, result, x); - } - for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { - applyRadialSection(i_theta, result, x); - } - } - /* Multi-threaded execution */ - else { - const int num_circle_tasks = grid.numberSmootherCircles(); - const int additional_radial_tasks = grid.ntheta() % 3; - const int num_radial_tasks = grid.ntheta() - additional_radial_tasks; - - #pragma omp parallel num_threads(num_omp_threads) - { - /* Circle Section 0 */ - #pragma omp for - for (int circle_task = 0; circle_task < num_circle_tasks; circle_task += 3) { - int i_r = grid.numberSmootherCircles() - circle_task - 1; - applyCircleSection(i_r, result, x); - } - /* Circle Section 1 */ - #pragma omp for - for (int circle_task = 1; circle_task < num_circle_tasks; circle_task += 3) { - int i_r = grid.numberSmootherCircles() - circle_task - 1; - applyCircleSection(i_r, result, x); - } - /* Circle Section 2 */ - #pragma omp for nowait - for (int circle_task = 2; circle_task < num_circle_tasks; circle_task += 3) { - int i_r = grid.numberSmootherCircles() - circle_task - 1; - applyCircleSection(i_r, result, x); - } - - /* Radial Section 0 */ - #pragma omp for - for (int radial_task = 0; radial_task < num_radial_tasks; radial_task += 3) { - if (radial_task > 0) { - int i_theta = radial_task + additional_radial_tasks; - applyRadialSection(i_theta, result, x); - } - else { - if (additional_radial_tasks == 0) { - applyRadialSection(0, result, x); - } - else if (additional_radial_tasks >= 1) { - applyRadialSection(0, result, x); - applyRadialSection(1, result, x); - } - } - } - /* Radial Section 1 */ - #pragma omp for - for (int radial_task = 1; radial_task < num_radial_tasks; radial_task += 3) { - if (radial_task > 1) { - int i_theta = radial_task + additional_radial_tasks; - applyRadialSection(i_theta, result, x); - } - else { - if (additional_radial_tasks == 0) { - applyRadialSection(1, result, x); - } - else if (additional_radial_tasks == 1) { - applyRadialSection(2, result, x); - } - else if (additional_radial_tasks == 2) { - applyRadialSection(2, result, x); - applyRadialSection(3, result, x); - } - } - } - /* Radial Section 2 */ - #pragma omp for - for (int radial_task = 2; radial_task < num_radial_tasks; radial_task += 3) { - int i_theta = radial_task + additional_radial_tasks; - applyRadialSection(i_theta, result, x); - } - } - } -} -// clang-format on - /* ------------------ */ /* result = rhs - A*x */ template @@ -117,10 +18,9 @@ void ResidualGive::computeResidual(Vector result, ConstV applySystemOperator(result, x); // Subtract A*x from rhs to get the residual. - const int n = result.size(); - const int num_omp_threads = Residual::num_omp_threads_; -#pragma omp parallel for num_threads(num_omp_threads) - for (int i = 0; i < n; i++) { - result[i] = rhs[i] - result[i]; - } -} + const int n = result.size(); + + Kokkos::parallel_for( + "Residual Give: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), + KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); +} \ No newline at end of file diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index ab50c20e..da5e1e18 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -3,10 +3,11 @@ namespace residual_take { -static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, - Vector& result, ConstVector& x, ConstVector& arr, - ConstVector& att, ConstVector& art, ConstVector& detDF, - ConstVector& coeff_beta) +static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_theta, const PolarGrid& grid, + bool DirBC_Interior, Vector& result, + ConstVector& x, ConstVector& arr, + ConstVector& att, ConstVector& art, + ConstVector& detDF, ConstVector& coeff_beta) { const int center = grid.index(i_r, i_theta); diff --git a/include/Residual/ResidualTake/residualTake.h b/include/Residual/ResidualTake/residualTake.h index 2c7cb173..d24b1880 100644 --- a/include/Residual/ResidualTake/residualTake.h +++ b/include/Residual/ResidualTake/residualTake.h @@ -13,11 +13,11 @@ class ResidualTake : public Residual const int num_omp_threads); ~ResidualTake() override = default; - void computeResidual(Vector result, ConstVector rhs, ConstVector x) const override; - - void applySystemOperator(Vector result, ConstVector x) const override; + void applySystemOperator(Vector result, ConstVector x) const final; + void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; }; #include "residualTake.inl" #include "applyATake.inl" + } // namespace gmgpolar diff --git a/include/Residual/ResidualTake/residualTake.inl b/include/Residual/ResidualTake/residualTake.inl index ab73cbda..40e6dfa9 100644 --- a/include/Residual/ResidualTake/residualTake.inl +++ b/include/Residual/ResidualTake/residualTake.inl @@ -15,14 +15,12 @@ void ResidualTake::computeResidual(Vector result, ConstV { assert(result.size() == x.size()); - const int num_omp_threads = Residual::num_omp_threads_; - applySystemOperator(result, x); // Subtract A*x from rhs to get the residual. const int n = result.size(); -#pragma omp parallel for num_threads(num_omp_threads) - for (int i = 0; i < n; i++) { - result[i] = rhs[i] - result[i]; - } + + Kokkos::parallel_for( + "Residual Take: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), + KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); } diff --git a/include/Residual/residual.h b/include/Residual/residual.h index f1525bd1..0b110ab7 100644 --- a/include/Residual/residual.h +++ b/include/Residual/residual.h @@ -26,10 +26,9 @@ class Residual } virtual ~Residual() = default; + virtual void applySystemOperator(Vector result, ConstVector x) const = 0; virtual void computeResidual(Vector result, ConstVector rhs, ConstVector x) const = 0; - virtual void applySystemOperator(Vector result, ConstVector x) const = 0; - protected: /* ------------------- */ /* Constructor members */ From ad070f1f8049a135c57aa11460db20779714204f Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 11 May 2026 23:21:08 +0200 Subject: [PATCH 2/8] skibidi --- include/Residual/ResidualGive/applyAGive.inl | 28 ++++++++++---------- include/Residual/ResidualTake/applyATake.inl | 9 +++---- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index a97fad28..565f7288 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -10,12 +10,12 @@ static inline void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid /* ---------------------------------------- */ /* Compute or retrieve stencil coefficients */ /* ---------------------------------------- */ - const int center = grid.index(i_r, i_theta); - const double r = grid.radius(i_r); - const double theta = grid.theta(i_theta); + const int center = grid.index(i_r, i_theta); + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); double coeff_beta, arr, att, art, detDF; - level_cache.obtainValues(i_r, i_theta, center, r, theta, coeff_beta, arr, att, art, detDF); + level_cache.obtainValues(i_r, i_theta, center, radius, theta, coeff_beta, arr, att, art, detDF); /* -------------------- */ /* Node in the interior */ @@ -83,14 +83,14 @@ static inline void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid /* ------------------------------------------------ */ if (DirBC_Interior) { /* Fill result(i,j) */ - result[center] += x[center]; + result[center] = x[center]; /* Give value to the interior nodes! */ - double h2 = grid.radialSpacing(i_r); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + const double h2 = grid.radialSpacing(i_r); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); - double coeff2 = 0.5 * (k1 + k2) / h2; + const double coeff2 = 0.5 * (k1 + k2) / h2; const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); @@ -168,14 +168,14 @@ static inline void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid /* ----------------------------- */ else if (i_r == grid.nr() - 1) { /* Fill result of (i,j) */ - result[center] += x[center]; + result[center] = x[center]; /* Give value to the interior nodes! */ - double h1 = grid.radialSpacing(i_r - 1); - double k1 = grid.angularSpacing(i_theta - 1); - double k2 = grid.angularSpacing(i_theta); + const double h1 = grid.radialSpacing(i_r - 1); + const double k1 = grid.angularSpacing(i_theta - 1); + const double k2 = grid.angularSpacing(i_theta); - double coeff1 = 0.5 * (k1 + k2) / h1; + const double coeff1 = 0.5 * (k1 + k2) / h1; const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index da5e1e18..b5c0a5f8 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -28,12 +28,11 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_ const double coeff4 = 0.5 * (h1 + h2) / k2; const double coeff5 = 0.25 * (h1 + h2) * (k1 + k2); - const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); - const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1); + const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1); + const int i_theta_Across = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2); - // Across origin: (i_r - 1, i_theta) gets replaced with (i_r, i_theta + (grid.ntheta() / 2)). - const int left = - (i_r > 0) ? grid.index(i_r - 1, i_theta) : grid.index(i_r, grid.wrapThetaIndex(i_theta + grid.ntheta() / 2)); + const int left = (i_r == 0) ? grid.index(i_r, i_theta_Across) : grid.index(i_r - 1, i_theta); const int right = grid.index(i_r + 1, i_theta); const int bottom = grid.index(i_r, i_theta_M1); const int top = grid.index(i_r, i_theta_P1); From 0d0db2940b6191b16c0f95e6062f560f6e80ea4c Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 11 May 2026 23:39:54 +0200 Subject: [PATCH 3/8] Remove --- include/Residual/ResidualGive/residualGive.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/include/Residual/ResidualGive/residualGive.h b/include/Residual/ResidualGive/residualGive.h index 4bbe254b..4aec6dbb 100644 --- a/include/Residual/ResidualGive/residualGive.h +++ b/include/Residual/ResidualGive/residualGive.h @@ -15,10 +15,6 @@ class ResidualGive : public Residual void applySystemOperator(Vector result, ConstVector x) const final; void computeResidual(Vector result, ConstVector rhs, ConstVector x) const final; - -private: - void applyCircleSection(const int i_r, Vector result, ConstVector x) const; - void applyRadialSection(const int i_theta, Vector result, ConstVector x) const; }; #include "residualGive.inl" From 4f9b114d0ad0d2a2802a8c352dc80bb7e62b8b58 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Tue, 12 May 2026 17:34:01 +0200 Subject: [PATCH 4/8] Update residualGive.inl --- include/Residual/ResidualGive/residualGive.inl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/Residual/ResidualGive/residualGive.inl b/include/Residual/ResidualGive/residualGive.inl index 122efdf1..0f1648d2 100644 --- a/include/Residual/ResidualGive/residualGive.inl +++ b/include/Residual/ResidualGive/residualGive.inl @@ -23,4 +23,6 @@ void ResidualGive::computeResidual(Vector result, ConstV Kokkos::parallel_for( "Residual Give: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); + + Kokkos::fence(); } \ No newline at end of file From 2e1d7208bfd74e23989d7d5c67c9f2d3257f2c70 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Tue, 12 May 2026 17:34:36 +0200 Subject: [PATCH 5/8] Update residualTake.inl --- include/Residual/ResidualTake/residualTake.inl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/Residual/ResidualTake/residualTake.inl b/include/Residual/ResidualTake/residualTake.inl index 40e6dfa9..dfe32e08 100644 --- a/include/Residual/ResidualTake/residualTake.inl +++ b/include/Residual/ResidualTake/residualTake.inl @@ -23,4 +23,6 @@ void ResidualTake::computeResidual(Vector result, ConstV Kokkos::parallel_for( "Residual Take: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); + + Kokkos::fence(); } From 1f30b1856a05979a64846498cc64fac010a62180 Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Tue, 12 May 2026 21:13:11 +0200 Subject: [PATCH 6/8] Adjust i_theta calculation for radial tasks --- include/Residual/ResidualGive/applyAGive.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 565f7288..3692f493 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -252,7 +252,7 @@ void ResidualGive::applySystemOperator(Vector result, Co Kokkos::parallel_for( "ResidualGive: ApplyA (Radial)", Kokkos::RangePolicy<>(0, num_radial_batches), KOKKOS_LAMBDA(const int radial_task) { - const int i_theta = start_radial + radial_task * 3; + const int i_theta = additional_radial_tasks + start_radial + radial_task * 3; for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { node_apply_a_give(i_r, i_theta, grid, level_cache, DirBC_Interior, result, x); } From 4f46c01ac2b0d97224e8f964c1b975bc085d4af0 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Wed, 13 May 2026 15:26:31 +0200 Subject: [PATCH 7/8] Fix cuda compilation (#257) --- include/Definitions/geometry_helper.h | 6 ++- .../poissonCoefficients.h | 7 +-- .../sonnendruckerCoefficients.h | 7 +-- .../sonnendruckerGyroCoefficients.h | 7 +-- .../zoniCoefficients.h | 7 +-- .../zoniGyroCoefficients.h | 7 +-- .../zoniShiftedCoefficients.h | 7 +-- .../zoniShiftedGyroCoefficients.h | 7 +-- .../DomainGeometry/circularGeometry.h | 13 +++--- .../DomainGeometry/circularGeometry.inl | 12 ++--- .../DomainGeometry/culhamGeometry.h | 46 ++++++++++--------- .../DomainGeometry/culhamGeometry.inl | 46 ++++++++++--------- .../DomainGeometry/czarnyGeometry.h | 13 +++--- .../DomainGeometry/czarnyGeometry.inl | 12 ++--- .../DomainGeometry/shafranovGeometry.h | 13 +++--- .../DomainGeometry/shafranovGeometry.inl | 12 ++--- include/Level/level.h | 9 ++-- include/Residual/ResidualGive/applyAGive.inl | 17 ++++--- .../Residual/ResidualGive/residualGive.inl | 4 +- .../Residual/ResidualTake/residualTake.inl | 2 +- .../DensityProfileCoefficients/CMakeLists.txt | 3 +- .../poissonCoefficients.cpp | 6 +-- .../sonnendruckerCoefficients.cpp | 6 +-- .../sonnendruckerGyroCoefficients.cpp | 6 +-- .../zoniCoefficients.cpp | 6 +-- .../zoniGyroCoefficients.cpp | 6 +-- .../zoniShiftedCoefficients.cpp | 6 +-- .../zoniShiftedGyroCoefficients.cpp | 6 +-- .../DomainGeometry/CMakeLists.txt | 3 +- 29 files changed, 163 insertions(+), 139 deletions(-) diff --git a/include/Definitions/geometry_helper.h b/include/Definitions/geometry_helper.h index 29821c9e..b3c645a7 100644 --- a/include/Definitions/geometry_helper.h +++ b/include/Definitions/geometry_helper.h @@ -4,10 +4,12 @@ #include "../InputFunctions/domainGeometry.h" #include +#include template -inline void compute_jacobian_elements(const DomainGeometry& domain_geometry, double r, double theta, double coeff_alpha, - double& arr, double& att, double& art, double& detDF) +KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& domain_geometry, double r, double theta, + double coeff_alpha, double& arr, double& att, double& art, + double& detDF) { /* Calculate the elements of the Jacobian matrix for the transformation mapping */ /* The Jacobian matrix is: */ diff --git a/include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h index 4cbf3d60..6f5d2e97 100644 --- a/include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class PoissonCoefficients PoissonCoefficients() = default; explicit PoissonCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.h index 62bc1a42..2b8953b8 100644 --- a/include/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class SonnendruckerCoefficients SonnendruckerCoefficients() = default; explicit SonnendruckerCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.h index 083d8361..6207372e 100644 --- a/include/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class SonnendruckerGyroCoefficients SonnendruckerGyroCoefficients() = default; explicit SonnendruckerGyroCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/zoniCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/zoniCoefficients.h index db14a255..d4845af1 100644 --- a/include/InputFunctions/DensityProfileCoefficients/zoniCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/zoniCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class ZoniCoefficients ZoniCoefficients() = default; explicit ZoniCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.h index d06a5cd0..9f48a6ee 100644 --- a/include/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class ZoniGyroCoefficients ZoniGyroCoefficients() = default; explicit ZoniGyroCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.h index 802dee8f..3a590a23 100644 --- a/include/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class ZoniShiftedCoefficients ZoniShiftedCoefficients() = default; explicit ZoniShiftedCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.h b/include/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.h index 80c25373..9f18ea0d 100644 --- a/include/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.h +++ b/include/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../densityProfileCoefficients.h" @@ -13,10 +14,10 @@ class ZoniShiftedGyroCoefficients ZoniShiftedGyroCoefficients() = default; explicit ZoniShiftedGyroCoefficients(double Rmax, double alpha); - double alpha(double r, double theta) const; - double beta(double r, double theta) const; + KOKKOS_FUNCTION double alpha(double r, double theta) const; + KOKKOS_FUNCTION double beta(double r, double theta) const; - double getAlphaJump() const; + KOKKOS_FUNCTION double getAlphaJump() const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DomainGeometry/circularGeometry.h b/include/InputFunctions/DomainGeometry/circularGeometry.h index 37815cde..081d86c4 100644 --- a/include/InputFunctions/DomainGeometry/circularGeometry.h +++ b/include/InputFunctions/DomainGeometry/circularGeometry.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../domainGeometry.h" @@ -13,12 +14,12 @@ class CircularGeometry CircularGeometry() = default; explicit CircularGeometry(double Rmax); - double Fx(double r, double theta) const; - double Fy(double r, double theta) const; - double dFx_dr(double r, double theta) const; - double dFy_dr(double r, double theta) const; - double dFx_dt(double r, double theta) const; - double dFy_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fx(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fy(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dt(double r, double theta) const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DomainGeometry/circularGeometry.inl b/include/InputFunctions/DomainGeometry/circularGeometry.inl index d80c482d..129cf91d 100644 --- a/include/InputFunctions/DomainGeometry/circularGeometry.inl +++ b/include/InputFunctions/DomainGeometry/circularGeometry.inl @@ -3,37 +3,37 @@ #include "circularGeometry.h" // In earlier versions denoted by 'x' -inline double CircularGeometry::Fx(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::Fx(double r, double theta) const { return (r / Rmax) * std::cos(theta); } // In earlier versions denoted by 'y' -inline double CircularGeometry::Fy(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::Fy(double r, double theta) const { return (r / Rmax) * std::sin(theta); } // In earlier versions denoted by 'Jrr' -inline double CircularGeometry::dFx_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::dFx_dr(double r, double theta) const { return (std::cos(theta)) / Rmax; } // In earlier versions denoted by 'Jtr' -inline double CircularGeometry::dFy_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::dFy_dr(double r, double theta) const { return (std::sin(theta)) / Rmax; } // In earlier versions denoted by 'Jrt' -inline double CircularGeometry::dFx_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::dFx_dt(double r, double theta) const { return (-(r / Rmax)) * std::sin(theta); } // In earlier versions denoted by 'Jtt' -inline double CircularGeometry::dFy_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CircularGeometry::dFy_dt(double r, double theta) const { return (r / Rmax) * std::cos(theta); } diff --git a/include/InputFunctions/DomainGeometry/culhamGeometry.h b/include/InputFunctions/DomainGeometry/culhamGeometry.h index c4e9543e..48005f18 100644 --- a/include/InputFunctions/DomainGeometry/culhamGeometry.h +++ b/include/InputFunctions/DomainGeometry/culhamGeometry.h @@ -4,6 +4,7 @@ #include #include #include +#include #include "../domainGeometry.h" @@ -16,34 +17,35 @@ class CulhamGeometry CulhamGeometry(); explicit CulhamGeometry(double Rmax); - double Fx(double r, double theta) const; - double Fy(double r, double theta) const; - double dFx_dr(double r, double theta) const; - double dFy_dr(double r, double theta) const; - double dFx_dt(double r, double theta) const; - double dFy_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fx(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fy(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dt(double r, double theta) const; private: const double Rmax = 1.3; void initializeGeometry(); - double my_sum(std::array& f, int64_t start_idx, int64_t end_idx) const; - double q(double rr) const; - double dq(double rr) const; - double p(double rr) const; - double dp(double rr) const; - double dg(double rr, double g) const; - double double_deriv(double rr, double c, double g, double dg, double val, double d_val) const; - double g(double rr) const; - double Delta(double rr) const; - double Delta_prime(double rr) const; - double E(double rr) const; - double T(double rr) const; - double E_prime(double rr) const; - double T_prime(double rr) const; - double P(double rr) const; - double dP(double rr) const; + KOKKOS_INLINE_FUNCTION double my_sum(std::array& f, int64_t start_idx, int64_t end_idx) const; + KOKKOS_INLINE_FUNCTION double q(double rr) const; + KOKKOS_INLINE_FUNCTION double dq(double rr) const; + KOKKOS_INLINE_FUNCTION double p(double rr) const; + KOKKOS_INLINE_FUNCTION double dp(double rr) const; + KOKKOS_INLINE_FUNCTION double dg(double rr, double g) const; + KOKKOS_INLINE_FUNCTION double double_deriv(double rr, double c, double g, double dg, double val, + double d_val) const; + KOKKOS_INLINE_FUNCTION double g(double rr) const; + KOKKOS_INLINE_FUNCTION double Delta(double rr) const; + KOKKOS_INLINE_FUNCTION double Delta_prime(double rr) const; + KOKKOS_INLINE_FUNCTION double E(double rr) const; + KOKKOS_INLINE_FUNCTION double T(double rr) const; + KOKKOS_INLINE_FUNCTION double E_prime(double rr) const; + KOKKOS_INLINE_FUNCTION double T_prime(double rr) const; + KOKKOS_INLINE_FUNCTION double P(double rr) const; + KOKKOS_INLINE_FUNCTION double dP(double rr) const; double rr; double dr; diff --git a/include/InputFunctions/DomainGeometry/culhamGeometry.inl b/include/InputFunctions/DomainGeometry/culhamGeometry.inl index 1def1f72..6f629573 100644 --- a/include/InputFunctions/DomainGeometry/culhamGeometry.inl +++ b/include/InputFunctions/DomainGeometry/culhamGeometry.inl @@ -3,7 +3,7 @@ #include "culhamGeometry.h" // In earlier versions denoted by 'x' -inline double CulhamGeometry::Fx(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::Fx(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -13,7 +13,7 @@ inline double CulhamGeometry::Fx(double r, double theta) const } // In earlier versions denoted by 'y' -inline double CulhamGeometry::Fy(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::Fy(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -23,7 +23,7 @@ inline double CulhamGeometry::Fy(double r, double theta) const } // In earlier versions denoted by 'Jrr' -inline double CulhamGeometry::dFx_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dFx_dr(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -34,7 +34,7 @@ inline double CulhamGeometry::dFx_dr(double r, double theta) const } // In earlier versions denoted by 'Jtr' -inline double CulhamGeometry::dFy_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dFy_dr(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -45,7 +45,7 @@ inline double CulhamGeometry::dFy_dr(double r, double theta) const } // In earlier versions denoted by 'Jrt' -inline double CulhamGeometry::dFx_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dFx_dt(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -55,7 +55,7 @@ inline double CulhamGeometry::dFx_dt(double r, double theta) const } // In earlier versions denoted by 'Jtt' -inline double CulhamGeometry::dFy_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dFy_dt(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -64,7 +64,8 @@ inline double CulhamGeometry::dFy_dt(double r, double theta) const 2.0 * T((r / Rmax)) * cos_two_theta; } -inline double CulhamGeometry::my_sum(std::array& f, int64_t start_idx, int64_t end_idx) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::my_sum(std::array& f, int64_t start_idx, + int64_t end_idx) const { int64_t i; double result; @@ -75,27 +76,27 @@ inline double CulhamGeometry::my_sum(std::array& f, int64_t start_ return result; } -inline double CulhamGeometry::q(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::q(double rr) const { return 0.8 - 0.1 * (rr * rr); } -inline double CulhamGeometry::dq(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dq(double rr) const { return (-0.2) * rr; } -inline double CulhamGeometry::p(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::p(double rr) const { return 100000.0 - 90000.0 * (rr * rr); } -inline double CulhamGeometry::dp(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dp(double rr) const { return (-180000.0) * rr; } -inline double CulhamGeometry::dg(double rr, double g) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dg(double rr, double g) const { return ((-g) * (0.0625000000000001 * (rr * rr) / pow((1.0 - 0.125 * (rr * rr)), 2.0) + 2.0 / (4.0 - 0.5 * (rr * rr))) + @@ -103,7 +104,8 @@ inline double CulhamGeometry::dg(double rr, double g) const (rr / (4.0 - 0.5 * (rr * rr)) + (4.0 - 0.5 * (rr * rr)) / (g * rr)); } -inline double CulhamGeometry::double_deriv(double rr, double c, double g, double dg, double val, double d_val) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::double_deriv(double rr, double c, double g, double dg, double val, + double d_val) const { return c * val / (rr * rr) - d_val * (pow(rr, (double)((-1))) + (4.0 - 0.5 * (rr * rr)) * @@ -113,7 +115,7 @@ inline double CulhamGeometry::double_deriv(double rr, double c, double g, double (g * rr)); } -inline double CulhamGeometry::g(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::g(double rr) const { int64_t ri; double dr; @@ -131,7 +133,7 @@ inline double CulhamGeometry::g(double rr) const } } -inline double CulhamGeometry::Delta(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::Delta(double rr) const { int64_t ri; double dr; @@ -149,7 +151,7 @@ inline double CulhamGeometry::Delta(double rr) const } } -inline double CulhamGeometry::Delta_prime(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::Delta_prime(double rr) const { int64_t ri; double dr; @@ -167,7 +169,7 @@ inline double CulhamGeometry::Delta_prime(double rr) const } } -inline double CulhamGeometry::E(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::E(double rr) const { int64_t ri; double dr; @@ -185,7 +187,7 @@ inline double CulhamGeometry::E(double rr) const } } -inline double CulhamGeometry::T(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::T(double rr) const { int64_t ri; double dr; @@ -203,7 +205,7 @@ inline double CulhamGeometry::T(double rr) const } } -inline double CulhamGeometry::E_prime(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::E_prime(double rr) const { int64_t ri; double dr; @@ -221,7 +223,7 @@ inline double CulhamGeometry::E_prime(double rr) const } } -inline double CulhamGeometry::T_prime(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::T_prime(double rr) const { int64_t ri; double dr; @@ -239,7 +241,7 @@ inline double CulhamGeometry::T_prime(double rr) const } } -inline double CulhamGeometry::P(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::P(double rr) const { if (rr == 0) { return 0.0; @@ -249,7 +251,7 @@ inline double CulhamGeometry::P(double rr) const } } -inline double CulhamGeometry::dP(double rr) const +KOKKOS_INLINE_FUNCTION double CulhamGeometry::dP(double rr) const { if (rr == 0) { return 0.0; diff --git a/include/InputFunctions/DomainGeometry/czarnyGeometry.h b/include/InputFunctions/DomainGeometry/czarnyGeometry.h index 9cd3fd4f..e0e50003 100644 --- a/include/InputFunctions/DomainGeometry/czarnyGeometry.h +++ b/include/InputFunctions/DomainGeometry/czarnyGeometry.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../domainGeometry.h" @@ -15,12 +16,12 @@ class CzarnyGeometry explicit CzarnyGeometry(); explicit CzarnyGeometry(double Rmax, double inverse_aspect_ratio_epsilon, double ellipticity_e); - double Fx(double r, double theta) const; - double Fy(double r, double theta) const; - double dFx_dr(double r, double theta) const; - double dFy_dr(double r, double theta) const; - double dFx_dt(double r, double theta) const; - double dFy_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fx(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fy(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dt(double r, double theta) const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DomainGeometry/czarnyGeometry.inl b/include/InputFunctions/DomainGeometry/czarnyGeometry.inl index 5556fcc8..3e71b6b3 100644 --- a/include/InputFunctions/DomainGeometry/czarnyGeometry.inl +++ b/include/InputFunctions/DomainGeometry/czarnyGeometry.inl @@ -3,7 +3,7 @@ #include "czarnyGeometry.h" // In earlier versions denoted by 'x' -inline double CzarnyGeometry::Fx(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::Fx(double r, double theta) const { double cos_theta = std::cos(theta); return (1.0 - @@ -12,7 +12,7 @@ inline double CzarnyGeometry::Fx(double r, double theta) const } // In earlier versions denoted by 'y' -inline double CzarnyGeometry::Fy(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::Fy(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -22,7 +22,7 @@ inline double CzarnyGeometry::Fy(double r, double theta) const } // In earlier versions denoted by 'Jrr' -inline double CzarnyGeometry::dFx_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::dFx_dr(double r, double theta) const { double cos_theta = std::cos(theta); return -(cos_theta) / @@ -31,7 +31,7 @@ inline double CzarnyGeometry::dFx_dr(double r, double theta) const } // In earlier versions denoted by 'Jtr' -inline double CzarnyGeometry::dFy_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::dFy_dr(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -43,7 +43,7 @@ inline double CzarnyGeometry::dFy_dr(double r, double theta) const } // In earlier versions denoted by 'Jrt' -inline double CzarnyGeometry::dFx_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::dFx_dt(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); @@ -52,7 +52,7 @@ inline double CzarnyGeometry::dFx_dt(double r, double theta) const } // In earlier versions denoted by 'Jtt' -inline double CzarnyGeometry::dFy_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double CzarnyGeometry::dFy_dt(double r, double theta) const { double sin_theta = std::sin(theta); double cos_theta = std::cos(theta); diff --git a/include/InputFunctions/DomainGeometry/shafranovGeometry.h b/include/InputFunctions/DomainGeometry/shafranovGeometry.h index cf87a15c..cd2a594d 100644 --- a/include/InputFunctions/DomainGeometry/shafranovGeometry.h +++ b/include/InputFunctions/DomainGeometry/shafranovGeometry.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include "../domainGeometry.h" @@ -15,12 +16,12 @@ class ShafranovGeometry ShafranovGeometry() = default; explicit ShafranovGeometry(double Rmax, double elongation_kappa, double shift_delta); - double Fx(double r, double theta) const; - double Fy(double r, double theta) const; - double dFx_dr(double r, double theta) const; - double dFy_dr(double r, double theta) const; - double dFx_dt(double r, double theta) const; - double dFy_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fx(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double Fy(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dr(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFx_dt(double r, double theta) const; + KOKKOS_INLINE_FUNCTION double dFy_dt(double r, double theta) const; private: const double Rmax = 1.3; diff --git a/include/InputFunctions/DomainGeometry/shafranovGeometry.inl b/include/InputFunctions/DomainGeometry/shafranovGeometry.inl index 690c4dca..e388ccfa 100644 --- a/include/InputFunctions/DomainGeometry/shafranovGeometry.inl +++ b/include/InputFunctions/DomainGeometry/shafranovGeometry.inl @@ -3,42 +3,42 @@ #include "shafranovGeometry.h" // In earlier versions denoted by 'x' -inline double ShafranovGeometry::Fx(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::Fx(double r, double theta) const { double cos_theta = std::cos(theta); return (1.0 - elongation_kappa) * (r / Rmax) * cos_theta - shift_delta * (r / Rmax) * (r / Rmax); } // In earlier versions denoted by 'y' -inline double ShafranovGeometry::Fy(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::Fy(double r, double theta) const { double sin_theta = std::sin(theta); return (1.0 + elongation_kappa) * (r / Rmax) * sin_theta; } // In earlier versions denoted by 'Jrr' -inline double ShafranovGeometry::dFx_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::dFx_dr(double r, double theta) const { double cos_theta = std::cos(theta); return ((Rmax - elongation_kappa * Rmax) * cos_theta - 2.0 * shift_delta * r) / (Rmax * Rmax); } // In earlier versions denoted by 'Jtr' -inline double ShafranovGeometry::dFy_dr(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::dFy_dr(double r, double theta) const { double sin_theta = std::sin(theta); return (elongation_kappa + 1.0) * sin_theta / Rmax; } // In earlier versions denoted by 'Jrt' -inline double ShafranovGeometry::dFx_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::dFx_dt(double r, double theta) const { double sin_theta = std::sin(theta); return ((elongation_kappa - 1.0) * r * sin_theta) / Rmax; } // In earlier versions denoted by 'Jtt' -inline double ShafranovGeometry::dFy_dt(double r, double theta) const +KOKKOS_INLINE_FUNCTION double ShafranovGeometry::dFy_dt(double r, double theta) const { double cos_theta = std::cos(theta); return ((elongation_kappa + 1.0) * r * cos_theta) / Rmax; diff --git a/include/Level/level.h b/include/Level/level.h index a791a7f1..9c1dd32f 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -19,6 +19,8 @@ class ExtrapolatedSmoother; #include #include +#include + #include "../PolarGrid/polargrid.h" #include "../InputFunctions/boundaryConditions.h" @@ -153,8 +155,9 @@ class LevelCache ConstVector art() const; ConstVector detDF() const; - inline void obtainValues(const int i_r, const int i_theta, const int global_index, double r, double theta, - double& coeff_beta, double& arr, double& att, double& art, double& detDF) const + KOKKOS_INLINE_FUNCTION void obtainValues(const int i_r, const int i_theta, const int global_index, double r, + double theta, double& coeff_beta, double& arr, double& att, double& art, + double& detDF) const { coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index] : density_profile_coefficients_.beta(r, theta); @@ -191,4 +194,4 @@ class LevelCache #include "levelCache.inl" #include "level.inl" -} // namespace gmgpolar \ No newline at end of file +} // namespace gmgpolar diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 3692f493..08c26658 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -1,11 +1,11 @@ #pragma once - namespace residual_give { template -static inline void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid, const LevelCacheType& level_cache, - bool DirBC_Interior, Vector& result, ConstVector& x) +static KOKKOS_INLINE_FUNCTION void node_apply_a_give(int i_r, int i_theta, const PolarGrid& grid, + const LevelCacheType& level_cache, bool DirBC_Interior, + Vector& result, ConstVector& x) { /* ---------------------------------------- */ /* Compute or retrieve stencil coefficients */ @@ -217,7 +217,8 @@ void ResidualGive::applySystemOperator(Vector result, Co for (int start_circle = 0; start_circle < 3; ++start_circle) { const int num_circular_tasks = (num_circle_tasks - start_circle + 2) / 3; Kokkos::parallel_for( - "ResidualGive: ApplyA (Circular)", Kokkos::RangePolicy<>(0, num_circular_tasks), + "ResidualGive: ApplyA (Circular)", + Kokkos::RangePolicy(0, num_circular_tasks), KOKKOS_LAMBDA(const int circle_task) { const int i_r = start_circle + circle_task * 3; for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) { @@ -239,7 +240,8 @@ void ResidualGive::applySystemOperator(Vector result, Co for (int i_theta = 0; i_theta < additional_radial_tasks; i_theta++) { Kokkos::parallel_for( - "ResidualGive: ApplyA (Radial, additional)", Kokkos::RangePolicy<>(0, 1), KOKKOS_LAMBDA(const int) { + "ResidualGive: ApplyA (Radial, additional)", Kokkos::RangePolicy(0, 1), + KOKKOS_LAMBDA(const int) { for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { node_apply_a_give(i_r, i_theta, grid, level_cache, DirBC_Interior, result, x); } @@ -250,7 +252,8 @@ void ResidualGive::applySystemOperator(Vector result, Co for (int start_radial = 0; start_radial < 3; ++start_radial) { const int num_radial_batches = (num_radial_tasks - start_radial + 2) / 3; Kokkos::parallel_for( - "ResidualGive: ApplyA (Radial)", Kokkos::RangePolicy<>(0, num_radial_batches), + "ResidualGive: ApplyA (Radial)", + Kokkos::RangePolicy(0, num_radial_batches), KOKKOS_LAMBDA(const int radial_task) { const int i_theta = additional_radial_tasks + start_radial + radial_task * 3; for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) { @@ -260,4 +263,4 @@ void ResidualGive::applySystemOperator(Vector result, Co Kokkos::fence(); } } -// clang-format on \ No newline at end of file +// clang-format on diff --git a/include/Residual/ResidualGive/residualGive.inl b/include/Residual/ResidualGive/residualGive.inl index 0f1648d2..00313f26 100644 --- a/include/Residual/ResidualGive/residualGive.inl +++ b/include/Residual/ResidualGive/residualGive.inl @@ -21,8 +21,8 @@ void ResidualGive::computeResidual(Vector result, ConstV const int n = result.size(); Kokkos::parallel_for( - "Residual Give: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), + "Residual Give: Subtract A*x from rhs", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); Kokkos::fence(); -} \ No newline at end of file +} diff --git a/include/Residual/ResidualTake/residualTake.inl b/include/Residual/ResidualTake/residualTake.inl index dfe32e08..df6bc7d3 100644 --- a/include/Residual/ResidualTake/residualTake.inl +++ b/include/Residual/ResidualTake/residualTake.inl @@ -21,7 +21,7 @@ void ResidualTake::computeResidual(Vector result, ConstV const int n = result.size(); Kokkos::parallel_for( - "Residual Take: Subtract A*x from rhs", Kokkos::RangePolicy<>(0, n), + "Residual Take: Subtract A*x from rhs", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(const int i) { result[i] = rhs[i] - result[i]; }); Kokkos::fence(); diff --git a/src/InputFunctions/DensityProfileCoefficients/CMakeLists.txt b/src/InputFunctions/DensityProfileCoefficients/CMakeLists.txt index 505bb9f8..0eef7d69 100644 --- a/src/InputFunctions/DensityProfileCoefficients/CMakeLists.txt +++ b/src/InputFunctions/DensityProfileCoefficients/CMakeLists.txt @@ -12,4 +12,5 @@ add_library(InputFunctions_DensityProfileCoefficients STATIC ${DENSITY_PROFILE_C target_include_directories(InputFunctions_DensityProfileCoefficients PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/include/InputFunctions/DensityProfileCoefficients -) \ No newline at end of file +) +target_link_libraries(InputFunctions_DensityProfileCoefficients PUBLIC Kokkos::kokkos) diff --git a/src/InputFunctions/DensityProfileCoefficients/poissonCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/poissonCoefficients.cpp index 94affce1..ab66832c 100644 --- a/src/InputFunctions/DensityProfileCoefficients/poissonCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/poissonCoefficients.cpp @@ -7,17 +7,17 @@ PoissonCoefficients::PoissonCoefficients(double Rmax, double alpha_jump) { } -double PoissonCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double PoissonCoefficients::alpha(double r, double theta) const { return 1.0; } -double PoissonCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double PoissonCoefficients::beta(double r, double theta) const { return 0.0; } -double PoissonCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double PoissonCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.cpp index 2cc66f2d..a1b74f77 100644 --- a/src/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.cpp @@ -7,17 +7,17 @@ SonnendruckerCoefficients::SonnendruckerCoefficients(double Rmax, double alpha_j { } -double SonnendruckerCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double SonnendruckerCoefficients::alpha(double r, double theta) const { return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r / Rmax) - 11.1111111111111); } -double SonnendruckerCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double SonnendruckerCoefficients::beta(double r, double theta) const { return 0.0; } -double SonnendruckerCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double SonnendruckerCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.cpp index 768bdac5..6b701c1f 100644 --- a/src/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.cpp @@ -7,18 +7,18 @@ SonnendruckerGyroCoefficients::SonnendruckerGyroCoefficients(double _Rmax, doubl { } -double SonnendruckerGyroCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double SonnendruckerGyroCoefficients::alpha(double r, double theta) const { return 0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r / Rmax) - 11.1111111111111); } -double SonnendruckerGyroCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double SonnendruckerGyroCoefficients::beta(double r, double theta) const { return pow((0.452961672473868 - 0.348432055749129 * atan(14.4444444444444 * (r / Rmax) - 11.1111111111111)), (double)((-1))); } -double SonnendruckerGyroCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double SonnendruckerGyroCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/zoniCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/zoniCoefficients.cpp index 3065415f..83b795ff 100644 --- a/src/InputFunctions/DensityProfileCoefficients/zoniCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/zoniCoefficients.cpp @@ -7,17 +7,17 @@ ZoniCoefficients::ZoniCoefficients(double Rmax, double alpha_jump) { } -double ZoniCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double ZoniCoefficients::alpha(double r, double theta) const { return exp(-tanh(10.0 * (r / Rmax) - 5.0)); } -double ZoniCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double ZoniCoefficients::beta(double r, double theta) const { return 0.0; } -double ZoniCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double ZoniCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.cpp index aba20f60..52e12213 100644 --- a/src/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.cpp @@ -7,17 +7,17 @@ ZoniGyroCoefficients::ZoniGyroCoefficients(double Rmax, double alpha_jump) { } -double ZoniGyroCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double ZoniGyroCoefficients::alpha(double r, double theta) const { return exp(-tanh(10.0 * (r / Rmax) - 5.0)); } -double ZoniGyroCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double ZoniGyroCoefficients::beta(double r, double theta) const { return exp(tanh(10.0 * (r / Rmax) - 5.0)); } -double ZoniGyroCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double ZoniGyroCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.cpp index 7d7a0607..1e63e4b7 100644 --- a/src/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.cpp @@ -7,17 +7,17 @@ ZoniShiftedCoefficients::ZoniShiftedCoefficients(double Rmax, double alpha_jump) { } -double ZoniShiftedCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double ZoniShiftedCoefficients::alpha(double r, double theta) const { return exp(-tanh(20.0 * (r / Rmax) - 14.0)); } -double ZoniShiftedCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double ZoniShiftedCoefficients::beta(double r, double theta) const { return 0.0; } -double ZoniShiftedCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double ZoniShiftedCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.cpp b/src/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.cpp index 1f3555a5..8fdb1388 100644 --- a/src/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.cpp +++ b/src/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.cpp @@ -7,17 +7,17 @@ ZoniShiftedGyroCoefficients::ZoniShiftedGyroCoefficients(double Rmax, double alp { } -double ZoniShiftedGyroCoefficients::alpha(double r, double theta) const +KOKKOS_FUNCTION double ZoniShiftedGyroCoefficients::alpha(double r, double theta) const { return exp(-tanh(20.0 * (r / Rmax) - 14.0)); } -double ZoniShiftedGyroCoefficients::beta(double r, double theta) const +KOKKOS_FUNCTION double ZoniShiftedGyroCoefficients::beta(double r, double theta) const { return exp(tanh(20.0 * (r / Rmax) - 14.0)); } -double ZoniShiftedGyroCoefficients::getAlphaJump() const +KOKKOS_FUNCTION double ZoniShiftedGyroCoefficients::getAlphaJump() const { return alpha_jump; } diff --git a/src/InputFunctions/DomainGeometry/CMakeLists.txt b/src/InputFunctions/DomainGeometry/CMakeLists.txt index 65b2c8b3..d83d211b 100644 --- a/src/InputFunctions/DomainGeometry/CMakeLists.txt +++ b/src/InputFunctions/DomainGeometry/CMakeLists.txt @@ -9,4 +9,5 @@ add_library(InputFunctions_DomainGeometry STATIC ${DOMAIN_GEOMETRY_SOURCES}) target_include_directories(InputFunctions_DomainGeometry PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR}/include/InputFunctions/DomainGeometry -) \ No newline at end of file +) +target_link_libraries(InputFunctions_DomainGeometry PUBLIC Kokkos::kokkos) From db8033436d871cf529ed1a8f09c4e86168c8c60f Mon Sep 17 00:00:00 2001 From: Julian Litz <91479202+julianlitz@users.noreply.github.com> Date: Wed, 13 May 2026 15:38:31 +0200 Subject: [PATCH 8/8] Update applyAGive.inl --- include/Residual/ResidualGive/applyAGive.inl | 1 - 1 file changed, 1 deletion(-) diff --git a/include/Residual/ResidualGive/applyAGive.inl b/include/Residual/ResidualGive/applyAGive.inl index 08c26658..338ea40c 100644 --- a/include/Residual/ResidualGive/applyAGive.inl +++ b/include/Residual/ResidualGive/applyAGive.inl @@ -263,4 +263,3 @@ void ResidualGive::applySystemOperator(Vector result, Co Kokkos::fence(); } } -// clang-format on