From c3152b193e3a6eec2e3359591d4544f6ba20f961 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 18 May 2026 22:25:21 +0200 Subject: [PATCH 1/4] Test --- include/Definitions/geometry_helper.h | 12 ++-- .../DirectSolverTake/applySymmetryShift.inl | 4 -- .../DirectSolverTake/buildSolverMatrix.inl | 68 +++++++++---------- .../applyAscOrtho.inl | 63 ++++++----------- .../buildInnerBoundaryAsc.inl | 18 +++-- .../buildTridiagonalAsc.inl | 67 +++++++++--------- include/GMGPolar/setup.h | 8 +-- include/Level/level.h | 14 ++++ include/Residual/ResidualTake/applyATake.inl | 25 ++++--- .../Smoother/SmootherTake/applyAscOrtho.inl | 59 ++++++---------- .../SmootherTake/buildInnerBoundaryAsc.inl | 29 ++++---- .../SmootherTake/buildTridiagonalAsc.inl | 37 +++++----- scripts/tutorial/run.sh | 6 +- 13 files changed, 183 insertions(+), 227 deletions(-) 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..a6a8b4fe9 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -21,7 +21,7 @@ maxOpenMPThreads=1 # Stencil distribution method: # 0 - CPU "Take": Each node independently applies the stencil # 1 - CPU "Give": The stencil operation is distributed across adjacent neighboring nodes -stencilDistributionMethod=0 +stencilDistributionMethod=1 # Caching behavior: # 0 - Recompute values on each iteration: Uses less memory but results in slower execution. # 1 - Reuse cached values: Consumes more memory but significantly improves performance. @@ -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 From 64b786e8e0a5ca44b3ca629afb881c008c8620fb Mon Sep 17 00:00:00 2001 From: julianlitz Date: Mon, 18 May 2026 22:29:53 +0200 Subject: [PATCH 2/4] Test --- scripts/tutorial/run.sh | 2 +- tests/GMGPolar/pcg_tests.cpp | 2 +- tests/GMGPolar/solve_tests.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index a6a8b4fe9..5a0dc98e3 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -25,7 +25,7 @@ stencilDistributionMethod=1 # Caching behavior: # 0 - Recompute values on each iteration: Uses less memory but results in slower execution. # 1 - Reuse cached values: Consumes more memory but significantly improves performance. -cacheDensityProfileCoefficients=1 +cacheDensityProfileCoefficients=0 cacheDomainGeometry=1 # Note: In the "Take" approach (stencilDistributionMethod=0), # caching is required for optimal performance, 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 From 6f61588c439fa2e8ea80fd0c9af2297c86348e63 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 19 May 2026 01:07:25 +0200 Subject: [PATCH 3/4] Compilation --- include/ConfigParser/config_parser.h | 4 +- src/CMakeLists.txt | 4 + src/ConfigParser/config_parser.cpp | 209 ++++- src/ConfigParser/select_solver_circular.cpp | 243 ++++++ src/ConfigParser/select_solver_culham.cpp | 78 ++ src/ConfigParser/select_solver_czarny.cpp | 242 ++++++ src/ConfigParser/select_solver_shafranov.cpp | 243 ++++++ src/ConfigParser/select_test_case.cpp | 775 ------------------- 8 files changed, 1008 insertions(+), 790 deletions(-) create mode 100644 src/ConfigParser/select_solver_circular.cpp create mode 100644 src/ConfigParser/select_solver_culham.cpp create mode 100644 src/ConfigParser/select_solver_czarny.cpp create mode 100644 src/ConfigParser/select_solver_shafranov.cpp 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/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/config_parser.cpp b/src/ConfigParser/config_parser.cpp index d9739af09..7ed16ee92 100644 --- a/src/ConfigParser/config_parser.cpp +++ b/src/ConfigParser/config_parser.cpp @@ -1,4 +1,6 @@ #include "../include/ConfigParser/config_parser.h" +#include "../include/GMGPolar/gmgpolar.h" + using namespace gmgpolar; ConfigParser::ConfigParser() @@ -67,10 +69,8 @@ ConfigParser::ConfigParser() parse(0, nullptr); } -// Parses command-line arguments and initializes all parameters bool ConfigParser::parse(int argc, char* argv[]) { - if (argc != 0) { try { parser_.parse_check(argc, argv); @@ -82,7 +82,7 @@ bool ConfigParser::parse(int argc, char* argv[]) } } - // Parse general parameters from command-line arguments + // Parse general parameters verbose_ = parser_.get("verbose"); paraview_ = parser_.get("paraview") != 0; max_omp_threads_ = parser_.get("maxOpenMPThreads"); @@ -98,7 +98,7 @@ bool ConfigParser::parse(int argc, char* argv[]) cache_density_profile_coefficients_ = parser_.get("cacheDensityProfileCoefficients") != 0; cache_domain_geometry_ = parser_.get("cacheDomainGeometry") != 0; - // Parse grid parameters from command-line arguments + // Parse grid parameters double R0 = parser_.get("R0"); double Rmax = parser_.get("Rmax"); int nr_exp = parser_.get("nr_exp"); @@ -106,7 +106,7 @@ bool ConfigParser::parse(int argc, char* argv[]) int anisotropic_factor = parser_.get("anisotropic_factor"); int divideBy2 = parser_.get("divideBy2"); - // Parse multigrid parameters from command-line arguments + // Parse multigrid parameters FMG_ = parser_.get("FMG") != 0; FMG_iterations_ = parser_.get("FMG_iterations"); const int FMG_cycleValue = parser_.get("FMG_cycle"); @@ -315,7 +315,195 @@ bool ConfigParser::parse(int argc, char* argv[]) return true; } -// Control Parameters +/* -------------------------------------------------------------------------- */ +/* solver() */ +/* -------------------------------------------------------------------------- */ +std::unique_ptr ConfigParser::solver() const +{ + const PolarGrid& grid = grid_; + + return std::visit( + [&grid](auto const& domain_geometry, auto const& density_profile_coefficients) { + using DomainGeomType = std::decay_t; + using DensityProfileCoefficientsType = std::decay_t; + + std::unique_ptr solver = + std::make_unique>( + grid, domain_geometry, density_profile_coefficients); + + return solver; + }, + *domain_geometry_, *density_profile_coefficients_); +} + +/* -------------------------------------------------------------------------- */ +/* selectTestCase() */ +/* -------------------------------------------------------------------------- */ +void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType problem_type, AlphaCoeff alpha_type, + BetaCoeff beta_type, double Rmax, double kappa_eps, double delta_e, double alpha_jump) +{ + geometry_type_ = geometry_type; + problem_type_ = problem_type; + alpha_type_ = alpha_type; + beta_type_ = beta_type; + Rmax_ = Rmax; + kappa_eps_ = kappa_eps; + delta_e_ = delta_e; + alpha_jump_ = alpha_jump; + + /* --------------- */ + /* Domain Geometry */ + switch (geometry_type) { + case GeometryType::CIRCULAR: + domain_geometry_ = std::make_unique(CircularGeometry(Rmax)); + break; + case GeometryType::SHAFRANOV: + domain_geometry_ = std::make_unique(ShafranovGeometry(Rmax, kappa_eps, delta_e)); + break; + case GeometryType::CZARNY: + domain_geometry_ = std::make_unique(CzarnyGeometry(Rmax, kappa_eps, delta_e)); + break; + case GeometryType::CULHAM: + domain_geometry_ = std::make_unique(CulhamGeometry(Rmax)); + break; + default: + throw std::runtime_error("Invalid geometry.\n"); + } + + /* ---------------------------- */ + /* Density Profile Coefficients */ + switch (alpha_type) { + case AlphaCoeff::POISSON: + density_profile_coefficients_ = + std::make_unique(PoissonCoefficients(Rmax, alpha_jump)); + break; + case AlphaCoeff::SONNENDRUCKER: + switch (beta_type) { + case BetaCoeff::ZERO: + density_profile_coefficients_ = + std::make_unique(SonnendruckerCoefficients(Rmax, alpha_jump)); + break; + case BetaCoeff::ALPHA_INVERSE: + density_profile_coefficients_ = + std::make_unique(SonnendruckerGyroCoefficients(Rmax, alpha_jump)); + break; + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI: + switch (beta_type) { + case BetaCoeff::ZERO: + density_profile_coefficients_ = + std::make_unique(ZoniCoefficients(Rmax, alpha_jump)); + break; + case BetaCoeff::ALPHA_INVERSE: + density_profile_coefficients_ = + std::make_unique(ZoniGyroCoefficients(Rmax, alpha_jump)); + break; + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + case AlphaCoeff::ZONI_SHIFTED: + switch (beta_type) { + case BetaCoeff::ZERO: + density_profile_coefficients_ = + std::make_unique(ZoniShiftedCoefficients(Rmax, alpha_jump)); + break; + case BetaCoeff::ALPHA_INVERSE: + density_profile_coefficients_ = + std::make_unique(ZoniShiftedGyroCoefficients(Rmax, alpha_jump)); + break; + default: + throw std::runtime_error("Invalid beta.\n"); + } + break; + default: + throw std::runtime_error("Invalid alpha.\n"); + } + + /* -------------- */ + /* Exact Solution */ + switch (problem_type) { + case ProblemType::CARTESIAN_R2: + switch (geometry_type) { + case GeometryType::CIRCULAR: + exact_solution_ = std::make_unique(Rmax); + break; + case GeometryType::SHAFRANOV: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CZARNY: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + default: + throw std::runtime_error("Invalid geometry for configuration.\n"); + } + break; + + case ProblemType::CARTESIAN_R6: + switch (geometry_type) { + case GeometryType::CIRCULAR: + exact_solution_ = std::make_unique(Rmax); + break; + case GeometryType::SHAFRANOV: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CZARNY: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + default: + throw std::runtime_error("Invalid geometry for configuration.\n"); + } + break; + + case ProblemType::POLAR_R6: + switch (geometry_type) { + case GeometryType::CIRCULAR: + exact_solution_ = std::make_unique(Rmax); + break; + case GeometryType::SHAFRANOV: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CZARNY: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CULHAM: + exact_solution_ = std::make_unique(Rmax); + break; + default: + throw std::runtime_error("Invalid geometry for configuration.\n"); + } + break; + + case ProblemType::REFINED_RADIUS: + switch (geometry_type) { + case GeometryType::CIRCULAR: + exact_solution_ = std::make_unique(Rmax); + break; + case GeometryType::SHAFRANOV: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CZARNY: + exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); + break; + case GeometryType::CULHAM: + exact_solution_ = std::make_unique(Rmax); + break; + default: + throw std::runtime_error("Invalid geometry for configuration.\n"); + } + break; + + default: + throw std::runtime_error("Invalid problem.\n"); + } +} + +/* -------------------------------------------------------------------------- */ +/* Getters */ +/* -------------------------------------------------------------------------- */ int ConfigParser::verbose() const { return verbose_; @@ -324,17 +512,14 @@ bool ConfigParser::paraview() const { return paraview_; } - int ConfigParser::maxOpenMPThreads() const { return max_omp_threads_; } - bool ConfigParser::DirBC_Interior() const { return DirBC_Interior_; } - StencilDistributionMethod ConfigParser::stencilDistributionMethod() const { return stencil_distribution_method_; @@ -347,13 +532,10 @@ bool ConfigParser::cacheDomainGeometry() const { return cache_domain_geometry_; } - const PolarGrid& ConfigParser::grid() const { return grid_; } - -// Multigrid Parameters bool ConfigParser::FMG() const { return FMG_; @@ -366,7 +548,6 @@ MultigridCycleType ConfigParser::FMG_cycle() const { return FMG_cycle_; } - bool ConfigParser::PCG() const { return PCG_; @@ -441,4 +622,4 @@ const DensityProfileCoefficientsVariant& ConfigParser::densityProfileCoefficient const ExactSolution& ConfigParser::exactSolution() const { return *exact_solution_.get(); -} +} \ No newline at end of file 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; From 94f358221e0f231beebf6620238d61406d946582 Mon Sep 17 00:00:00 2001 From: julianlitz Date: Tue, 19 May 2026 01:11:48 +0200 Subject: [PATCH 4/4] fix --- scripts/tutorial/run.sh | 4 +- src/ConfigParser/config_parser.cpp | 209 ++--------------------------- 2 files changed, 16 insertions(+), 197 deletions(-) diff --git a/scripts/tutorial/run.sh b/scripts/tutorial/run.sh index 5a0dc98e3..ee3592b2e 100755 --- a/scripts/tutorial/run.sh +++ b/scripts/tutorial/run.sh @@ -21,11 +21,11 @@ maxOpenMPThreads=1 # Stencil distribution method: # 0 - CPU "Take": Each node independently applies the stencil # 1 - CPU "Give": The stencil operation is distributed across adjacent neighboring nodes -stencilDistributionMethod=1 +stencilDistributionMethod=0 # Caching behavior: # 0 - Recompute values on each iteration: Uses less memory but results in slower execution. # 1 - Reuse cached values: Consumes more memory but significantly improves performance. -cacheDensityProfileCoefficients=0 +cacheDensityProfileCoefficients=1 cacheDomainGeometry=1 # Note: In the "Take" approach (stencilDistributionMethod=0), # caching is required for optimal performance, diff --git a/src/ConfigParser/config_parser.cpp b/src/ConfigParser/config_parser.cpp index 7ed16ee92..d9739af09 100644 --- a/src/ConfigParser/config_parser.cpp +++ b/src/ConfigParser/config_parser.cpp @@ -1,6 +1,4 @@ #include "../include/ConfigParser/config_parser.h" -#include "../include/GMGPolar/gmgpolar.h" - using namespace gmgpolar; ConfigParser::ConfigParser() @@ -69,8 +67,10 @@ ConfigParser::ConfigParser() parse(0, nullptr); } +// Parses command-line arguments and initializes all parameters bool ConfigParser::parse(int argc, char* argv[]) { + if (argc != 0) { try { parser_.parse_check(argc, argv); @@ -82,7 +82,7 @@ bool ConfigParser::parse(int argc, char* argv[]) } } - // Parse general parameters + // Parse general parameters from command-line arguments verbose_ = parser_.get("verbose"); paraview_ = parser_.get("paraview") != 0; max_omp_threads_ = parser_.get("maxOpenMPThreads"); @@ -98,7 +98,7 @@ bool ConfigParser::parse(int argc, char* argv[]) cache_density_profile_coefficients_ = parser_.get("cacheDensityProfileCoefficients") != 0; cache_domain_geometry_ = parser_.get("cacheDomainGeometry") != 0; - // Parse grid parameters + // Parse grid parameters from command-line arguments double R0 = parser_.get("R0"); double Rmax = parser_.get("Rmax"); int nr_exp = parser_.get("nr_exp"); @@ -106,7 +106,7 @@ bool ConfigParser::parse(int argc, char* argv[]) int anisotropic_factor = parser_.get("anisotropic_factor"); int divideBy2 = parser_.get("divideBy2"); - // Parse multigrid parameters + // Parse multigrid parameters from command-line arguments FMG_ = parser_.get("FMG") != 0; FMG_iterations_ = parser_.get("FMG_iterations"); const int FMG_cycleValue = parser_.get("FMG_cycle"); @@ -315,195 +315,7 @@ bool ConfigParser::parse(int argc, char* argv[]) return true; } -/* -------------------------------------------------------------------------- */ -/* solver() */ -/* -------------------------------------------------------------------------- */ -std::unique_ptr ConfigParser::solver() const -{ - const PolarGrid& grid = grid_; - - return std::visit( - [&grid](auto const& domain_geometry, auto const& density_profile_coefficients) { - using DomainGeomType = std::decay_t; - using DensityProfileCoefficientsType = std::decay_t; - - std::unique_ptr solver = - std::make_unique>( - grid, domain_geometry, density_profile_coefficients); - - return solver; - }, - *domain_geometry_, *density_profile_coefficients_); -} - -/* -------------------------------------------------------------------------- */ -/* selectTestCase() */ -/* -------------------------------------------------------------------------- */ -void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType problem_type, AlphaCoeff alpha_type, - BetaCoeff beta_type, double Rmax, double kappa_eps, double delta_e, double alpha_jump) -{ - geometry_type_ = geometry_type; - problem_type_ = problem_type; - alpha_type_ = alpha_type; - beta_type_ = beta_type; - Rmax_ = Rmax; - kappa_eps_ = kappa_eps; - delta_e_ = delta_e; - alpha_jump_ = alpha_jump; - - /* --------------- */ - /* Domain Geometry */ - switch (geometry_type) { - case GeometryType::CIRCULAR: - domain_geometry_ = std::make_unique(CircularGeometry(Rmax)); - break; - case GeometryType::SHAFRANOV: - domain_geometry_ = std::make_unique(ShafranovGeometry(Rmax, kappa_eps, delta_e)); - break; - case GeometryType::CZARNY: - domain_geometry_ = std::make_unique(CzarnyGeometry(Rmax, kappa_eps, delta_e)); - break; - case GeometryType::CULHAM: - domain_geometry_ = std::make_unique(CulhamGeometry(Rmax)); - break; - default: - throw std::runtime_error("Invalid geometry.\n"); - } - - /* ---------------------------- */ - /* Density Profile Coefficients */ - switch (alpha_type) { - case AlphaCoeff::POISSON: - density_profile_coefficients_ = - std::make_unique(PoissonCoefficients(Rmax, alpha_jump)); - break; - case AlphaCoeff::SONNENDRUCKER: - switch (beta_type) { - case BetaCoeff::ZERO: - density_profile_coefficients_ = - std::make_unique(SonnendruckerCoefficients(Rmax, alpha_jump)); - break; - case BetaCoeff::ALPHA_INVERSE: - density_profile_coefficients_ = - std::make_unique(SonnendruckerGyroCoefficients(Rmax, alpha_jump)); - break; - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI: - switch (beta_type) { - case BetaCoeff::ZERO: - density_profile_coefficients_ = - std::make_unique(ZoniCoefficients(Rmax, alpha_jump)); - break; - case BetaCoeff::ALPHA_INVERSE: - density_profile_coefficients_ = - std::make_unique(ZoniGyroCoefficients(Rmax, alpha_jump)); - break; - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - case AlphaCoeff::ZONI_SHIFTED: - switch (beta_type) { - case BetaCoeff::ZERO: - density_profile_coefficients_ = - std::make_unique(ZoniShiftedCoefficients(Rmax, alpha_jump)); - break; - case BetaCoeff::ALPHA_INVERSE: - density_profile_coefficients_ = - std::make_unique(ZoniShiftedGyroCoefficients(Rmax, alpha_jump)); - break; - default: - throw std::runtime_error("Invalid beta.\n"); - } - break; - default: - throw std::runtime_error("Invalid alpha.\n"); - } - - /* -------------- */ - /* Exact Solution */ - switch (problem_type) { - case ProblemType::CARTESIAN_R2: - switch (geometry_type) { - case GeometryType::CIRCULAR: - exact_solution_ = std::make_unique(Rmax); - break; - case GeometryType::SHAFRANOV: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CZARNY: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::CARTESIAN_R6: - switch (geometry_type) { - case GeometryType::CIRCULAR: - exact_solution_ = std::make_unique(Rmax); - break; - case GeometryType::SHAFRANOV: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CZARNY: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::POLAR_R6: - switch (geometry_type) { - case GeometryType::CIRCULAR: - exact_solution_ = std::make_unique(Rmax); - break; - case GeometryType::SHAFRANOV: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CZARNY: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CULHAM: - exact_solution_ = std::make_unique(Rmax); - break; - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - case ProblemType::REFINED_RADIUS: - switch (geometry_type) { - case GeometryType::CIRCULAR: - exact_solution_ = std::make_unique(Rmax); - break; - case GeometryType::SHAFRANOV: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CZARNY: - exact_solution_ = std::make_unique(Rmax, kappa_eps, delta_e); - break; - case GeometryType::CULHAM: - exact_solution_ = std::make_unique(Rmax); - break; - default: - throw std::runtime_error("Invalid geometry for configuration.\n"); - } - break; - - default: - throw std::runtime_error("Invalid problem.\n"); - } -} - -/* -------------------------------------------------------------------------- */ -/* Getters */ -/* -------------------------------------------------------------------------- */ +// Control Parameters int ConfigParser::verbose() const { return verbose_; @@ -512,14 +324,17 @@ bool ConfigParser::paraview() const { return paraview_; } + int ConfigParser::maxOpenMPThreads() const { return max_omp_threads_; } + bool ConfigParser::DirBC_Interior() const { return DirBC_Interior_; } + StencilDistributionMethod ConfigParser::stencilDistributionMethod() const { return stencil_distribution_method_; @@ -532,10 +347,13 @@ bool ConfigParser::cacheDomainGeometry() const { return cache_domain_geometry_; } + const PolarGrid& ConfigParser::grid() const { return grid_; } + +// Multigrid Parameters bool ConfigParser::FMG() const { return FMG_; @@ -548,6 +366,7 @@ MultigridCycleType ConfigParser::FMG_cycle() const { return FMG_cycle_; } + bool ConfigParser::PCG() const { return PCG_; @@ -622,4 +441,4 @@ const DensityProfileCoefficientsVariant& ConfigParser::densityProfileCoefficient const ExactSolution& ConfigParser::exactSolution() const { return *exact_solution_.get(); -} \ No newline at end of file +}