diff --git a/include/ConfigParser/config_parser.h b/include/ConfigParser/config_parser.h index 6bb961665..a737c76f0 100644 --- a/include/ConfigParser/config_parser.h +++ b/include/ConfigParser/config_parser.h @@ -10,12 +10,14 @@ #include "../../include/PolarGrid/polargrid.h" #include "../../include/GMGPolar/test_cases.h" #include "../../include/GMGPolar/igmgpolar.h" -#include "../../include/GMGPolar/gmgpolar.h" #include "test_selection.h" namespace gmgpolar { +template +class GMGPolar; + class ConfigParser { public: diff --git a/include/Definitions/geometry_helper.h b/include/Definitions/geometry_helper.h index b3c645a7a..1d4499c69 100644 --- a/include/Definitions/geometry_helper.h +++ b/include/Definitions/geometry_helper.h @@ -7,9 +7,9 @@ #include template -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) +KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& domain_geometry, const double r, + const double theta, const 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: */ @@ -26,9 +26,9 @@ KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& doma /* which is represented by: */ /* [arr, 0.5*art] */ /* [0.5*atr, att] */ - arr = 0.5 * (Jtt * Jtt + Jrt * Jrt) * coeff_alpha / std::fabs(detDF); - att = 0.5 * (Jtr * Jtr + Jrr * Jrr) * coeff_alpha / std::fabs(detDF); - art = (-Jtt * Jtr - Jrt * Jrr) * coeff_alpha / std::fabs(detDF); + arr = 0.5 * (Jtt * Jtt + Jrt * Jrt) * coeff_alpha / Kokkos::fabs(detDF); + att = 0.5 * (Jtr * Jtr + Jrr * Jrr) * coeff_alpha / Kokkos::fabs(detDF); + art = (-Jtt * Jtr - Jrt * Jrr) * coeff_alpha / Kokkos::fabs(detDF); /* Note that the inverse Jacobian matrix DF^{-1} is: */ /* 1.0 / det(DF) * */ /* [Jtt, -Jrt] */ diff --git a/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl b/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl index 1d5a0845e..042e5b8ef 100644 --- a/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl +++ b/include/DirectSolver/DirectSolverTake/applySymmetryShift.inl @@ -12,9 +12,7 @@ void DirectSolverTake::applySymmetryShiftInnerBoundary(HostVecto assert(DirectSolver::DirBC_Interior_); - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - HostConstVector arr = level_cache.arr(); HostConstVector att = level_cache.att(); HostConstVector art = level_cache.art(); @@ -55,9 +53,7 @@ void DirectSolverTake::applySymmetryShiftOuterBoundary(HostVecto const PolarGrid& grid = DirectSolver::grid_; const LevelCacheType& level_cache = DirectSolver::level_cache_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - HostConstVector arr = level_cache.arr(); HostConstVector att = level_cache.att(); HostConstVector art = level_cache.art(); diff --git a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl index b12897863..1956c8aaf 100644 --- a/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl +++ b/include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl @@ -30,8 +30,7 @@ template static KOKKOS_INLINE_FUNCTION void nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, const SystemMatrix& solver_matrix, HostConstVector& arr, HostConstVector& att, - HostConstVector& art, HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& art, HostConstVector& detDF, const double coeff_beta) { int ptr, offset; int row, column; @@ -71,13 +70,12 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - const double center_value = - (coeff5 * coeff_beta[center_index] * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + const double center_value = (coeff5 * coeff_beta * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ @@ -197,13 +195,12 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - const double center_value = - (coeff5 * coeff_beta[center_index] * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + const double center_value = (coeff5 * coeff_beta * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */ @@ -289,13 +286,12 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - const double center_value = - (coeff5 * coeff_beta[center_index] * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + const double center_value = (coeff5 * coeff_beta * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ @@ -397,13 +393,12 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri const double bottom_value = -coeff3 * (att(center_index) + att(bottom_index)); /* Bottom */ const double top_value = -coeff4 * (att(center_index) + att(top_index)); /* Top */ - const double center_value = - (coeff5 * coeff_beta[center_index] * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ - - left_value /* Center: (Left) */ - - right_value /* Center: (Right) */ - - bottom_value /* Center: (Bottom) */ - - top_value /* Center: (Top) */ - ); + const double center_value = (coeff5 * coeff_beta * Kokkos::fabs(detDF(center_index)) /* beta_{i,j} */ + - left_value /* Center: (Left) */ + - right_value /* Center: (Right) */ + - bottom_value /* Center: (Bottom) */ + - top_value /* Center: (Top) */ + ); const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */ const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */ @@ -503,14 +498,11 @@ typename DirectSolverTake::SystemMatrix DirectSolverTake solver_matrix(n, n, nnz_per_row); #endif - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); /* We split the loops into two regions to better respect the */ /* access patterns of the smoother and improve cache locality. */ @@ -524,6 +516,7 @@ typename DirectSolverTake::SystemMatrix DirectSolverTake::SystemMatrix DirectSolverTake& x, HostConstVector& rhs, HostVector& result, HostConstVector& arr, HostConstVector& att, HostConstVector& art, - HostConstVector& detDF, HostConstVector& coeff_beta) + HostConstVector& detDF) { KOKKOS_ASSERT(i_r >= 0 && i_r <= grid.numberSmootherCircles()); @@ -189,8 +189,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, HostConstVector& x, HostConstVector& rhs, HostVector& result, HostConstVector& arr, const HostConstVector& att, - HostConstVector& art, const HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& art, const HostConstVector& detDF) { assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr()); @@ -472,14 +471,11 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(H const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); /* The outer most circle next to the radial section is defined to be black. */ const int start_black_circles = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0; @@ -494,8 +490,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackCircleSection(H // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -511,16 +506,12 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(H const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; - const int num_omp_threads = ExtrapolatedSmoother::num_omp_threads_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); /* The outer most circle next to the radial section is defined to be black. */ const int start_white_circles = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1; @@ -535,8 +526,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteCircleSection(H // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -552,16 +542,12 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(H const PolarGrid& grid = ExtrapolatedSmoother::grid_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; - const int num_omp_threads = ExtrapolatedSmoother::num_omp_threads_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); assert(grid.ntheta() % 2 == 0); const int start_black_radials = 0; @@ -576,8 +562,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoBlackRadialSection(H // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -595,14 +580,11 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(H const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; const int num_omp_threads = ExtrapolatedSmoother::num_omp_threads_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); assert(grid.ntheta() % 2 == 0); const int start_white_radials = 1; @@ -617,8 +599,7 @@ void ExtrapolatedSmootherTake::applyAscOrthoWhiteRadialSection(H // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl index 0a5dd32c3..f2a54f3be 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildInnerBoundaryAsc.inl @@ -31,7 +31,7 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, const InnerBoundaryMatrix& matrix, HostConstVector& arr, HostConstVector& att, HostConstVector& art, - HostConstVector& detDF, HostConstVector& coeff_beta) + HostConstVector& detDF, const double coeff_beta) { using extrapolated_smoother_take::getCircleAscIndex; using extrapolated_smoother_take::getStencil; @@ -116,10 +116,10 @@ nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, const int top = grid.index(i_r, i_theta_P1); const int right = grid.index(i_r + 1, i_theta); - const double center_value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + + const double center_value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); - const double left_value = -coeff1 * (arr[center] + arr[left]); + const double left_value = -coeff1 * (arr[center] + arr[left]); /* Fill matrix row of (i,j) */ row = center_index; @@ -195,18 +195,16 @@ ExtrapolatedSmootherTake::buildInteriorBoundarySolverMatrix() SparseMatrixCSR inner_boundary_solver_matrix(ntheta, ntheta, nnz_per_row); #endif - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); Kokkos::parallel_for( "ExtrapolatedSmootherTake: BuildInnerBoundaryMatrix", Kokkos::RangePolicy(0, ntheta), KOKKOS_LAMBDA(const int i_theta) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildInteriorBoundarySolverMatrix(i_theta, grid, DirBC_Interior, inner_boundary_solver_matrix, arr, att, art, detDF, coeff_beta); }); diff --git a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl index 814d88e5d..3264500f4 100644 --- a/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl +++ b/include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/buildTridiagonalAsc.inl @@ -18,8 +18,7 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesCircleSecti const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, const BatchedTridiagonalSolver& circle_tridiagonal_solver, const BatchedTridiagonalSolver& radial_tridiagonal_solver, HostConstVector& arr, - HostConstVector& att, HostConstVector& art, HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& att, HostConstVector& art, HostConstVector& detDF, const double coeff_beta) { using extrapolated_smoother_take::updateMatrixElement; @@ -89,9 +88,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesCircleSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Bottom */ @@ -130,9 +129,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesCircleSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); } else { @@ -168,8 +167,7 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, const BatchedTridiagonalSolver& circle_tridiagonal_solver, const BatchedTridiagonalSolver& radial_tridiagonal_solver, HostConstVector& arr, - HostConstVector& att, HostConstVector& art, HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& att, HostConstVector& art, HostConstVector& detDF, const double coeff_beta) { using extrapolated_smoother_take::updateMatrixElement; @@ -242,9 +240,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -287,9 +285,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = 0.25 * (h1 + h2) * (k1 + k2) * coeff_beta[center] * Kokkos::fabs(detDF[center]) + - coeff1 * (arr[center] + arr[left]) + coeff2 * (arr[center] + arr[right]) + - coeff3 * (att[center] + att[bottom]) + coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); } else { @@ -351,9 +349,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Right */ @@ -374,9 +372,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); } else { @@ -441,9 +439,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -469,9 +467,9 @@ static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatricesRadialSecti /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + - coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + - coeff4 * (att[center] + att[top]); + value = coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) + coeff1 * (arr[center] + arr[left]) + + coeff2 * (arr[center] + arr[right]) + coeff3 * (att[center] + att[bottom]) + + coeff4 * (att[center] + att[top]); updateMatrixElement(solver, batch, row, column, value); } } @@ -539,14 +537,11 @@ void ExtrapolatedSmootherTake::buildTridiagonalSolverMatrices() const bool DirBC_Interior = ExtrapolatedSmoother::DirBC_Interior_; const LevelCacheType& level_cache = ExtrapolatedSmoother::level_cache_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); const BatchedTridiagonalSolver& circle_tridiagonal_solver = circle_tridiagonal_solver_; const BatchedTridiagonalSolver& radial_tridiagonal_solver = radial_tridiagonal_solver_; @@ -563,6 +558,7 @@ void ExtrapolatedSmootherTake::buildTridiagonalSolverMatrices() ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_r, const int i_theta) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildTridiagonalSolverMatricesCircleSection(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver, radial_tridiagonal_solver, arr, att, art, detDF, coeff_beta); @@ -577,6 +573,7 @@ void ExtrapolatedSmootherTake::buildTridiagonalSolverMatrices() ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_theta, const int i_r) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildTridiagonalSolverMatricesRadialSection(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver, radial_tridiagonal_solver, arr, att, art, detDF, coeff_beta); diff --git a/include/GMGPolar/setup.h b/include/GMGPolar/setup.h index eae98303b..2ffa83efc 100644 --- a/include/GMGPolar/setup.h +++ b/include/GMGPolar/setup.h @@ -10,11 +10,9 @@ void GMGPolar::setup() auto start_setup_createLevels = std::chrono::high_resolution_clock::now(); - if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE) { - if (!cache_density_profile_coefficients_ || !cache_domain_geometry_) { - throw std::runtime_error("Error: Caching must be enabled for both density profile coefficients and domain " - "geometry in 'Take' implementation strategy."); - } + if (stencil_distribution_method_ == StencilDistributionMethod::CPU_TAKE && !cache_domain_geometry_) { + throw std::runtime_error( + "Error: Caching must be enabled for domain geometry in 'Take' implementation strategy."); } // -------------------------------- // diff --git a/include/Level/level.h b/include/Level/level.h index 3390704d6..1b114bdb7 100644 --- a/include/Level/level.h +++ b/include/Level/level.h @@ -176,6 +176,20 @@ class LevelCache } } + // Write on just for obtaining coeff beta + KOKKOS_INLINE_FUNCTION double obtainBeta(const int i_r, const int i_theta, const PolarGrid& grid) const + { + if (cache_density_profile_coefficients_) { + const int index = grid.index(i_r, i_theta); + return coeff_beta_[index]; + } + else { + const double radius = grid.radius(i_r); + const double theta = grid.theta(i_theta); + return density_profile_coefficients_.beta(radius, theta); + } + } + private: const DomainGeometry& domain_geometry_; const DensityProfileCoefficients& density_profile_coefficients_; diff --git a/include/Residual/ResidualTake/applyATake.inl b/include/Residual/ResidualTake/applyATake.inl index 827c896df..d0a856d9f 100644 --- a/include/Residual/ResidualTake/applyATake.inl +++ b/include/Residual/ResidualTake/applyATake.inl @@ -4,11 +4,10 @@ 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, HostVector& result, + const bool DirBC_Interior, HostVector& result, HostConstVector& x, HostConstVector& arr, HostConstVector& att, HostConstVector& art, - HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& detDF, const double coeff_beta) { const int center = grid.index(i_r, i_theta); @@ -40,7 +39,7 @@ static KOKKOS_INLINE_FUNCTION void node_apply_a_take(const int i_r, const int i_ double value = 0.0; - value += coeff5 * coeff_beta[center] * std::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ + value += coeff5 * coeff_beta * std::fabs(detDF[center]) * x[center]; /* beta_{i,j} */ value += coeff1 * (arr[center] + arr[left]) * (x[center] - x[left]); /* Center: (Left) - Left */ value += coeff2 * (arr[center] + arr[right]) * (x[center] - x[right]); /* Center: (Right) - Right */ @@ -72,17 +71,15 @@ void ResidualTake::applySystemOperator(HostVector result using residual_take::node_apply_a_take; - const PolarGrid& grid = Residual::grid_; - const bool DirBC_Interior = Residual::DirBC_Interior_; + const PolarGrid& grid = Residual::grid_; + const bool DirBC_Interior = Residual::DirBC_Interior_; + const LevelCacheType& level_cache = Residual::level_cache_; - assert(Residual::level_cache_.cacheDensityProfileCoefficients()); assert(Residual::level_cache_.cacheDomainGeometry()); - - HostConstVector arr = Residual::level_cache_.arr(); - HostConstVector att = Residual::level_cache_.att(); - HostConstVector art = Residual::level_cache_.art(); - HostConstVector detDF = Residual::level_cache_.detDF(); - HostConstVector coeff_beta = Residual::level_cache_.coeff_beta(); + HostConstVector arr = Residual::level_cache_.arr(); + HostConstVector att = Residual::level_cache_.att(); + HostConstVector art = Residual::level_cache_.art(); + HostConstVector detDF = Residual::level_cache_.detDF(); /* We split the loops into two regions to better respect the */ /* access patterns of the smoother and improve cache locality. */ @@ -96,6 +93,7 @@ void ResidualTake::applySystemOperator(HostVector result ), // Kokkos lambda function to execute for each point in the index space KOKKOS_CLASS_LAMBDA(const int i_r, const int i_theta) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); }); @@ -108,6 +106,7 @@ void ResidualTake::applySystemOperator(HostVector result ), // Kokkos lambda function to execute for each point in the index space KOKKOS_CLASS_LAMBDA(const int i_theta, const int i_r) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); node_apply_a_take(i_r, i_theta, grid, DirBC_Interior, result, x, arr, att, art, detDF, coeff_beta); }); diff --git a/include/Smoother/SmootherTake/applyAscOrtho.inl b/include/Smoother/SmootherTake/applyAscOrtho.inl index d70542cd5..53df03f08 100644 --- a/include/Smoother/SmootherTake/applyAscOrtho.inl +++ b/include/Smoother/SmootherTake/applyAscOrtho.inl @@ -7,7 +7,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, HostConstVector& x, HostConstVector& rhs, HostVector& result, HostConstVector& arr, HostConstVector& att, HostConstVector& art, - HostConstVector& detDF, HostConstVector& coeff_beta) + HostConstVector& detDF) { KOKKOS_ASSERT(i_r >= 0 && i_r < grid.numberSmootherCircles()); @@ -59,8 +59,7 @@ static KOKKOS_INLINE_FUNCTION void nodeApplyAscOrthoRadialTake(const int i_r, const int i_theta, const PolarGrid& grid, bool DirBC_Interior, HostConstVector& x, HostConstVector& rhs, HostVector& result, HostConstVector& arr, const HostConstVector& att, - HostConstVector& art, const HostConstVector& detDF, - HostConstVector& coeff_beta) + HostConstVector& art, const HostConstVector& detDF) { KOKKOS_ASSERT(i_r >= grid.numberSmootherCircles() && i_r < grid.nr()); @@ -127,14 +126,11 @@ void SmootherTake::applyAscOrthoBlackCircleSection(HostConstVect const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); /* The outer most circle next to the radial section is defined to be black. */ const int start_black_circles = (grid.numberSmootherCircles() % 2 == 0) ? 1 : 0; @@ -149,8 +145,7 @@ void SmootherTake::applyAscOrthoBlackCircleSection(HostConstVect // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { int i_r = start_black_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -166,14 +161,11 @@ void SmootherTake::applyAscOrthoWhiteCircleSection(HostConstVect const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); /* The outer most circle next to the radial section is defined to be black. */ const int start_white_circles = (grid.numberSmootherCircles() % 2 == 0) ? 0 : 1; @@ -188,8 +180,7 @@ void SmootherTake::applyAscOrthoWhiteCircleSection(HostConstVect // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int circle_task, const int i_theta) { const int i_r = start_white_circles + circle_task * 2; - nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoCircleTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -205,14 +196,12 @@ void SmootherTake::applyAscOrthoBlackRadialSection(HostConstVect const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); assert(grid.ntheta() % 2 == 0); const int start_black_radials = 0; @@ -227,8 +216,7 @@ void SmootherTake::applyAscOrthoBlackRadialSection(HostConstVect // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_black_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); @@ -244,14 +232,12 @@ void SmootherTake::applyAscOrthoWhiteRadialSection(HostConstVect const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); assert(grid.ntheta() % 2 == 0); const int start_white_radials = 1; @@ -266,8 +252,7 @@ void SmootherTake::applyAscOrthoWhiteRadialSection(HostConstVect // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int radial_task, const int i_r) { const int i_theta = start_white_radials + radial_task * 2; - nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF, - coeff_beta); + nodeApplyAscOrthoRadialTake(i_r, i_theta, grid, DirBC_Interior, x, rhs, temp, arr, att, art, detDF); }); Kokkos::fence(); diff --git a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl index a7ed77d75..b0b12be61 100644 --- a/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl +++ b/include/Smoother/SmootherTake/buildInnerBoundaryAsc.inl @@ -8,8 +8,8 @@ namespace smoother_take #ifdef GMGPOLAR_USE_MUMPS // When using the MUMPS solver, the matrix is assembled in COO format. static KOKKOS_INLINE_FUNCTION void -update_CSR_COO_MatrixElement(const SparseMatrixCOO& matrix, int ptr, int offset, int row, - int column, double value) +update_CSR_COO_MatrixElement(const SparseMatrixCOO& matrix, const int ptr, const int offset, + const int row, const int column, const double value) { matrix.set_row_index(ptr + offset, row); matrix.set_col_index(ptr + offset, column); @@ -18,8 +18,8 @@ update_CSR_COO_MatrixElement(const SparseMatrixCOO& m #else // When using the in-house solver, the matrix is stored in CSR format. static KOKKOS_INLINE_FUNCTION void -update_CSR_COO_MatrixElement(const SparseMatrixCSR& matrix, int ptr, int offset, int row, - int column, double value) +update_CSR_COO_MatrixElement(const SparseMatrixCSR& matrix, const int ptr, const int offset, + const int row, const int column, const double value) { matrix.set_row_nz_index(row, offset, column); matrix.set_row_nz_entry(row, offset, value); @@ -28,10 +28,10 @@ update_CSR_COO_MatrixElement(const SparseMatrixCSR& m template static KOKKOS_INLINE_FUNCTION void -nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, bool DirBC_Interior, +nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, const InnerBoundaryMatrix& matrix, HostConstVector& arr, HostConstVector& att, HostConstVector& art, - HostConstVector& detDF, HostConstVector& coeff_beta) + HostConstVector& detDF, const double coeff_beta) { using smoother_take::getCircleAscIndex; using smoother_take::getStencil; @@ -106,8 +106,8 @@ nodeBuildInteriorBoundarySolverMatrix(const int i_theta, const PolarGrid& grid, const double bottom_value = -coeff3 * (att[center] + att[bottom]); const double top_value = -coeff4 * (att[center] + att[top]); - const double center_value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + const double center_value = + coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) - (left_value + right_value + bottom_value + top_value); /* Fill matrix row of (i,j) */ row = center_index; @@ -149,7 +149,6 @@ SmootherTake::buildInteriorBoundarySolverMatrix() const PolarGrid& grid = Smoother::grid_; const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - const int num_omp_threads = Smoother::num_omp_threads_; const int i_r = 0; const int ntheta = grid.ntheta(); @@ -172,18 +171,16 @@ SmootherTake::buildInteriorBoundarySolverMatrix() SparseMatrixCSR inner_boundary_solver_matrix(ntheta, ntheta, nnz_per_row); #endif - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); Kokkos::parallel_for( "SmootherTake: BuildInnerBoundaryMatrix", Kokkos::RangePolicy(0, ntheta), KOKKOS_LAMBDA(const int i_theta) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildInteriorBoundarySolverMatrix(i_theta, grid, DirBC_Interior, inner_boundary_solver_matrix, arr, att, art, detDF, coeff_beta); }); diff --git a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl index 02fe8020d..e1fd65f68 100644 --- a/include/Smoother/SmootherTake/buildTridiagonalAsc.inl +++ b/include/Smoother/SmootherTake/buildTridiagonalAsc.inl @@ -3,8 +3,8 @@ namespace smoother_take { -static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const BatchedTridiagonalSolver& solver, int batch, - int row, int column, double value) +static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const BatchedTridiagonalSolver& solver, const int batch, + const int row, const int column, const double value) { if (row == column) solver.set_main_diagonal(batch, row, value); @@ -14,13 +14,11 @@ static KOKKOS_INLINE_FUNCTION void updateMatrixElement(const BatchedTridiagonalS solver.set_cyclic_corner(batch, value); } -static KOKKOS_INLINE_FUNCTION void -nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior, - const BatchedTridiagonalSolver& circle_tridiagonal_solver, - const BatchedTridiagonalSolver& radial_tridiagonal_solver, - HostConstVector& arr, HostConstVector& att, - HostConstVector& art, HostConstVector& detDF, - HostConstVector& coeff_beta) +static KOKKOS_INLINE_FUNCTION void nodeBuildTridiagonalSolverMatrices( + const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior, + const BatchedTridiagonalSolver& circle_tridiagonal_solver, + const BatchedTridiagonalSolver& radial_tridiagonal_solver, HostConstVector& arr, + HostConstVector& att, HostConstVector& art, HostConstVector& detDF, const double coeff_beta) { using smoother_take::updateMatrixElement; @@ -75,8 +73,8 @@ nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = + coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) - (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Bottom */ @@ -130,8 +128,8 @@ nodeBuildTridiagonalSolverMatrices(int i_r, int i_theta, const PolarGrid& grid, /* Center: (Left, Right, Bottom, Top) */ row = center_index; column = center_index; - value = coeff5 * coeff_beta[center] * Kokkos::fabs(detDF[center]) - - (left_value + right_value + bottom_value + top_value); + value = + coeff5 * coeff_beta * Kokkos::fabs(detDF[center]) - (left_value + right_value + bottom_value + top_value); updateMatrixElement(solver, batch, row, column, value); /* Left */ @@ -194,14 +192,11 @@ void SmootherTake::buildTridiagonalSolverMatrices() const LevelCacheType& level_cache = Smoother::level_cache_; const bool DirBC_Interior = Smoother::DirBC_Interior_; - assert(level_cache.cacheDensityProfileCoefficients()); assert(level_cache.cacheDomainGeometry()); - - HostConstVector arr = level_cache.arr(); - HostConstVector att = level_cache.att(); - HostConstVector art = level_cache.art(); - HostConstVector detDF = level_cache.detDF(); - HostConstVector coeff_beta = level_cache.coeff_beta(); + HostConstVector arr = level_cache.arr(); + HostConstVector att = level_cache.att(); + HostConstVector art = level_cache.art(); + HostConstVector detDF = level_cache.detDF(); const BatchedTridiagonalSolver& circle_tridiagonal_solver = circle_tridiagonal_solver_; const BatchedTridiagonalSolver& radial_tridiagonal_solver = radial_tridiagonal_solver_; @@ -218,6 +213,7 @@ void SmootherTake::buildTridiagonalSolverMatrices() ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_r, const int i_theta) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildTridiagonalSolverMatrices(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver, radial_tridiagonal_solver, arr, att, art, detDF, coeff_beta); }); @@ -231,6 +227,7 @@ void SmootherTake::buildTridiagonalSolverMatrices() ), // Kokkos lambda function to execute for each point in the index space KOKKOS_LAMBDA(const int i_theta, const int i_r) { + const double coeff_beta = level_cache.obtainBeta(i_r, i_theta, grid); nodeBuildTridiagonalSolverMatrices(i_r, i_theta, grid, DirBC_Interior, circle_tridiagonal_solver, radial_tridiagonal_solver, arr, att, art, detDF, coeff_beta); }); diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index 45271a301..ee3592b2e 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -31,8 +31,8 @@ cacheDomainGeometry=1 # caching is required for optimal performance, # so both density profile coefficients and domain geometry need to be cached. if [[ $stencilDistributionMethod -eq 0 ]]; then - if [[ $cacheDensityProfileCoefficients -ne 1 || $cacheDomainGeometry -ne 1 ]]; then - echo "Error: Caching must be enabled (set to 1) for both density profile coefficients and domain geometry in 'Take' approach (stencilDistributionMethod=0)." + if [[ $cacheDomainGeometry -ne 1 ]]; then + echo "Error: Caching of the domain geometry must be enabled for 'Take' approach." exit 1 fi fi diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eff7760db..d292366b7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -44,6 +44,10 @@ set(INTERPOLATION_SOURCES set(CONFIG_PARSER_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/config_parser.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/select_test_case.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/select_solver_circular.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/select_solver_shafranov.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/select_solver_czarny.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ConfigParser/select_solver_culham.cpp ) # Create the main library diff --git a/src/ConfigParser/select_solver_circular.cpp b/src/ConfigParser/select_solver_circular.cpp new file mode 100644 index 000000000..cd754c062 --- /dev/null +++ b/src/ConfigParser/select_solver_circular.cpp @@ -0,0 +1,243 @@ +#include "../include/ConfigParser/config_parser.h" +#include "../include/GMGPolar/gmgpolar.h" + +using namespace gmgpolar; + +template +void ConfigParser::solve(GMGPolar& solver) const +{ + if constexpr (std::is_same_v) { + + switch (problem_type_) { + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R2: { + CartesianR2_Boundary_CircularGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR2_Poisson_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Sonnendrucker_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_SonnendruckerGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Zoni_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_ZoniShifted_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniShiftedGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R6: { + CartesianR6_Boundary_CircularGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR6_Poisson_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Sonnendrucker_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_SonnendruckerGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Zoni_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_ZoniShifted_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniShiftedGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::POLAR_R6: { + PolarR6_Boundary_CircularGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + PolarR6_Poisson_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Sonnendrucker_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_SonnendruckerGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Zoni_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_ZoniShifted_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniShiftedGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::REFINED_RADIUS: { + Refined_Boundary_CircularGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ALPHA_INVERSE: { + Refined_ZoniShiftedGyro_CircularGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta for configuration.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha for configuration.\n"); + } + break; + } + + default: + throw std::runtime_error("Invalid problem.\n"); + } + } +} + +// Explicit instantiations — CircularGeometry for all coefficient types +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; \ No newline at end of file diff --git a/src/ConfigParser/select_solver_culham.cpp b/src/ConfigParser/select_solver_culham.cpp new file mode 100644 index 000000000..54633952c --- /dev/null +++ b/src/ConfigParser/select_solver_culham.cpp @@ -0,0 +1,78 @@ +#include "../include/ConfigParser/config_parser.h" +#include "../include/GMGPolar/gmgpolar.h" + +using namespace gmgpolar; + +template +void ConfigParser::solve(GMGPolar& solver) const +{ + if constexpr (std::is_same_v) { + + switch (problem_type_) { + + /* ------------------------------------------------------------------ */ + /* Culham only supports POLAR_R6 and REFINED_RADIUS. */ + /* CARTESIAN_R2 and CARTESIAN_R6 are not defined for this geometry. */ + /* ------------------------------------------------------------------ */ + + case ProblemType::POLAR_R6: { + PolarR6_Boundary_CulhamGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniShiftedGyro_CulhamGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta for configuration.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha for configuration.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::REFINED_RADIUS: { + Refined_Boundary_CulhamGeometry bc(Rmax_); + switch (alpha_type_) { + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ALPHA_INVERSE: { + Refined_ZoniShiftedGyro_CulhamGeometry src(grid_, Rmax_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta for configuration.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha for configuration.\n"); + } + break; + } + + default: + throw std::runtime_error("Invalid geometry for configuration.\n"); + } + } +} + +// Explicit instantiations — CulhamGeometry for all coefficient types +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; \ No newline at end of file diff --git a/src/ConfigParser/select_solver_czarny.cpp b/src/ConfigParser/select_solver_czarny.cpp new file mode 100644 index 000000000..3a347c5dc --- /dev/null +++ b/src/ConfigParser/select_solver_czarny.cpp @@ -0,0 +1,242 @@ +#include "../include/ConfigParser/config_parser.h" +#include "../include/GMGPolar/gmgpolar.h" + +using namespace gmgpolar; + +template +void ConfigParser::solve(GMGPolar& solver) const +{ + if constexpr (std::is_same_v) { + + switch (problem_type_) { + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R2: { + CartesianR2_Boundary_CzarnyGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR2_Poisson_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Sonnendrucker_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_SonnendruckerGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Zoni_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_ZoniShifted_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniShiftedGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R6: { + CartesianR6_Boundary_CzarnyGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR6_Poisson_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Sonnendrucker_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_SonnendruckerGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Zoni_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_ZoniShifted_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniShiftedGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::POLAR_R6: { + PolarR6_Boundary_CzarnyGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + PolarR6_Poisson_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Sonnendrucker_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_SonnendruckerGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Zoni_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_ZoniShifted_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniShiftedGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::REFINED_RADIUS: { + Refined_Boundary_CzarnyGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ALPHA_INVERSE: { + Refined_ZoniShiftedGyro_CzarnyGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta for configuration.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha for configuration.\n"); + } + break; + } + + default: + throw std::runtime_error("Invalid problem.\n"); + } + } +} + +// Explicit instantiations — CzarnyGeometry for all coefficient types +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; \ No newline at end of file diff --git a/src/ConfigParser/select_solver_shafranov.cpp b/src/ConfigParser/select_solver_shafranov.cpp new file mode 100644 index 000000000..d58b616ca --- /dev/null +++ b/src/ConfigParser/select_solver_shafranov.cpp @@ -0,0 +1,243 @@ +#include "../include/ConfigParser/config_parser.h" +#include "../include/GMGPolar/gmgpolar.h" + +using namespace gmgpolar; + +template +void ConfigParser::solve(GMGPolar& solver) const +{ + if constexpr (std::is_same_v) { + + switch (problem_type_) { + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R2: { + CartesianR2_Boundary_ShafranovGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR2_Poisson_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Sonnendrucker_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_SonnendruckerGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_Zoni_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR2_ZoniShifted_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR2_ZoniShiftedGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::CARTESIAN_R6: { + CartesianR6_Boundary_ShafranovGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + CartesianR6_Poisson_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Sonnendrucker_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_SonnendruckerGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_Zoni_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + CartesianR6_ZoniShifted_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + CartesianR6_ZoniShiftedGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::POLAR_R6: { + PolarR6_Boundary_ShafranovGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::POISSON: { + PolarR6_Poisson_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Sonnendrucker_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_SonnendruckerGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_Zoni_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ZERO: { + PolarR6_ZoniShifted_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + case BetaCoeff::ALPHA_INVERSE: { + PolarR6_ZoniShiftedGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + break; + } + + /* ------------------------------------------------------------------ */ + case ProblemType::REFINED_RADIUS: { + Refined_Boundary_ShafranovGeometry bc(Rmax_, kappa_eps_, delta_e_); + switch (alpha_type_) { + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type_) { + case BetaCoeff::ALPHA_INVERSE: { + Refined_ZoniShiftedGyro_ShafranovGeometry src(grid_, Rmax_, kappa_eps_, delta_e_); + solver.solve(bc, src); + break; + } + default: + throw std::runtime_error("Invalid beta for configuration.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha for configuration.\n"); + } + break; + } + + default: + throw std::runtime_error("Invalid problem.\n"); + } + } +} + +// Explicit instantiations — ShafranovGeometry all coefficient types +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void +ConfigParser::solve(GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; +template void ConfigParser::solve( + GMGPolar&) const; \ No newline at end of file diff --git a/src/ConfigParser/select_test_case.cpp b/src/ConfigParser/select_test_case.cpp index cf12c8972..a594fee31 100644 --- a/src/ConfigParser/select_test_case.cpp +++ b/src/ConfigParser/select_test_case.cpp @@ -194,778 +194,3 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble throw std::runtime_error("Invalid problem.\n"); } } - -template -void ConfigParser::solve(GMGPolar& solver) const -{ - assert(domain_geometry_ != nullptr); - /* ------------------------------------------ */ - /* Source Term (rhs_f) and BoundaryConditions */ - switch (problem_type_) { - case ProblemType::CARTESIAN_R2: - - switch (geometry_type_) { - case GeometryType::CIRCULAR: { - CartesianR2_Boundary_CircularGeometry boundary_conditions(Rmax_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR2_Poisson_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Sonnendrucker_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_SonnendruckerGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Zoni_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_ZoniShifted_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniShiftedGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::SHAFRANOV: { - CartesianR2_Boundary_ShafranovGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR2_Poisson_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Sonnendrucker_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_SonnendruckerGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Zoni_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_ZoniShifted_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniShiftedGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::CZARNY: { - CartesianR2_Boundary_CzarnyGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR2_Poisson_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Sonnendrucker_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_SonnendruckerGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_Zoni_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR2_ZoniShifted_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR2_ZoniShiftedGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::CARTESIAN_R6: - - switch (geometry_type_) { - case GeometryType::CIRCULAR: { - CartesianR6_Boundary_CircularGeometry boundary_conditions(Rmax_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR6_Poisson_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Sonnendrucker_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_SonnendruckerGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Zoni_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_ZoniShifted_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniShiftedGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::SHAFRANOV: { - CartesianR6_Boundary_ShafranovGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR6_Poisson_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Sonnendrucker_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_SonnendruckerGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Zoni_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_ZoniShifted_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniShiftedGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::CZARNY: { - CartesianR6_Boundary_CzarnyGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - CartesianR6_Poisson_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Sonnendrucker_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_SonnendruckerGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_Zoni_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - CartesianR6_ZoniShifted_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - CartesianR6_ZoniShiftedGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::POLAR_R6: - - switch (geometry_type_) { - case GeometryType::CIRCULAR: { - PolarR6_Boundary_CircularGeometry boundary_conditions(Rmax_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - PolarR6_Poisson_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Sonnendrucker_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_SonnendruckerGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Zoni_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_ZoniShifted_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniShiftedGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::SHAFRANOV: { - PolarR6_Boundary_ShafranovGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - PolarR6_Poisson_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Sonnendrucker_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_SonnendruckerGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Zoni_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_ZoniShifted_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniShiftedGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::CZARNY: { - PolarR6_Boundary_CzarnyGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::POISSON: { - PolarR6_Poisson_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Sonnendrucker_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_SonnendruckerGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_Zoni_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ZERO: { - PolarR6_ZoniShifted_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniShiftedGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - break; - } - - case GeometryType::CULHAM: { - PolarR6_Boundary_CulhamGeometry boundary_conditions(Rmax_); - switch (alpha_type_) { - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ALPHA_INVERSE: { - PolarR6_ZoniShiftedGyro_CulhamGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta for configuration.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha for configuration.\n"); - } - break; - } - - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::REFINED_RADIUS: - - switch (geometry_type_) { - case GeometryType::CIRCULAR: { - Refined_Boundary_CircularGeometry boundary_conditions(Rmax_); - - switch (alpha_type_) { - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ALPHA_INVERSE: { - Refined_ZoniShiftedGyro_CircularGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta for configuration.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha for configuration.\n"); - } - break; - } - - case GeometryType::SHAFRANOV: { - Refined_Boundary_ShafranovGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ALPHA_INVERSE: { - Refined_ZoniShiftedGyro_ShafranovGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta for configuration.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha for configuration.\n"); - } - break; - } - - case GeometryType::CZARNY: { - Refined_Boundary_CzarnyGeometry boundary_conditions(Rmax_, kappa_eps_, delta_e_); - - switch (alpha_type_) { - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ALPHA_INVERSE: { - Refined_ZoniShiftedGyro_CzarnyGeometry source_term(grid_, Rmax_, kappa_eps_, delta_e_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta for configuration.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha for configuration.\n"); - } - break; - } - - case GeometryType::CULHAM: { - Refined_Boundary_CulhamGeometry boundary_conditions(Rmax_); - - switch (alpha_type_) { - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type_) { - case BetaCoeff::ALPHA_INVERSE: { - Refined_ZoniShiftedGyro_CulhamGeometry source_term(grid_, Rmax_); - solver.solve(boundary_conditions, source_term); - break; - } - default: - throw std::runtime_error("Invalid beta for configuration.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha for configuration.\n"); - } - break; - } - - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - default: - throw std::runtime_error("Invalid problem.\n"); - } -} - -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; - -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; - -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; - -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void -ConfigParser::solve(GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; -template void ConfigParser::solve( - GMGPolar& solver) const; diff --git a/tests/GMGPolar/pcg_tests.cpp b/tests/GMGPolar/pcg_tests.cpp index 7fe96b1ee..6b1d82f92 100644 --- a/tests/GMGPolar/pcg_tests.cpp +++ b/tests/GMGPolar/pcg_tests.cpp @@ -290,7 +290,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // maxOpenMPThreads std::integral_constant, // DirBC_Interior std::integral_constant, // StencilDistributionMethod - std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDensityProfileCoefficient std::integral_constant, // cacheDomainGeometry std::integral_constant, // extrapolation std::integral_constant, // maxLevels diff --git a/tests/GMGPolar/solve_tests.cpp b/tests/GMGPolar/solve_tests.cpp index 9eab97833..9ae4eb763 100644 --- a/tests/GMGPolar/solve_tests.cpp +++ b/tests/GMGPolar/solve_tests.cpp @@ -254,7 +254,7 @@ using gmgpolar_test_suite = testing::Types< std::integral_constant, // maxOpenMPThreads std::integral_constant, // DirBC_Interior std::integral_constant, // StencilDistributionMethod - std::integral_constant, // cacheDensityProfileCoefficient + std::integral_constant, // cacheDensityProfileCoefficient std::integral_constant, // cacheDomainGeometry std::integral_constant, // extrapolation std::integral_constant, // maxLevels