Skip to content

Allow Optional Runtime Computation of the Density Profile in the Take Approach #275

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
julianlitz wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Allow Optional Runtime Computation of the Density Profile in the Take Approach #275

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c3152b1
Test
julianlitz May 18, 2026
64b786e
Test
julianlitz May 18, 2026
6f61588
Compilation
julianlitz May 18, 2026
94f3582
fix
julianlitz May 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion include/ConfigParser/config_parser.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
class GMGPolar;

class ConfigParser
{
public:
Expand Down
12 changes: 6 additions & 6 deletions include/Definitions/geometry_helper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@
#include <Kokkos_Core.hpp>

template <concepts::DomainGeometry DomainGeometry>
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: */
Expand All @@ -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] */
Expand Down
4 changes: 0 additions & 4 deletions include/DirectSolver/DirectSolverTake/applySymmetryShift.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,7 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(HostVecto

assert(DirectSolver<LevelCacheType>::DirBC_Interior_);

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
Expand Down Expand Up @@ -55,9 +53,7 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(HostVecto
const PolarGrid& grid = DirectSolver<LevelCacheType>::grid_;
const LevelCacheType& level_cache = DirectSolver<LevelCacheType>::level_cache_;

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
Expand Down
68 changes: 31 additions & 37 deletions include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,7 @@ template <typename SystemMatrix>
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<double>& arr, HostConstVector<double>& att,
HostConstVector<double>& art, HostConstVector<double>& detDF,
HostConstVector<double>& coeff_beta)
HostConstVector<double>& art, HostConstVector<double>& detDF, const double coeff_beta)
{
int ptr, offset;
int row, column;
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -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 */
Expand Down Expand Up @@ -503,14 +498,11 @@ typename DirectSolverTake<LevelCacheType>::SystemMatrix DirectSolverTake<LevelCa
SparseMatrixCSR<double, Kokkos::HostSpace> solver_matrix(n, n, nnz_per_row);
#endif

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();
HostConstVector<double> coeff_beta = level_cache.coeff_beta();
HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();

/* We split the loops into two regions to better respect the */
/* access patterns of the smoother and improve cache locality. */
Expand All @@ -524,6 +516,7 @@ typename DirectSolverTake<LevelCacheType>::SystemMatrix DirectSolverTake<LevelCa
),
// 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);
nodeBuildSolverMatrixTake(i_r, i_theta, grid, DirBC_Interior, solver_matrix, arr, att, art, detDF,
coeff_beta);
});
Expand All @@ -537,6 +530,7 @@ typename DirectSolverTake<LevelCacheType>::SystemMatrix DirectSolverTake<LevelCa
),
// 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);
nodeBuildSolverMatrixTake(i_r, i_theta, grid, DirBC_Interior, solver_matrix, arr, att, art, detDF,
coeff_beta);
});
Expand Down
63 changes: 22 additions & 41 deletions include/ExtrapolatedSmoother/ExtrapolatedSmootherTake/applyAscOrtho.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ static KOKKOS_INLINE_FUNCTION void
nodeApplyAscOrthoCircleTake(const int i_r, const int i_theta, const PolarGrid& grid, const bool DirBC_Interior,
HostConstVector<double>& x, HostConstVector<double>& rhs, HostVector<double>& result,
HostConstVector<double>& arr, HostConstVector<double>& att, HostConstVector<double>& art,
HostConstVector<double>& detDF, HostConstVector<double>& coeff_beta)
HostConstVector<double>& detDF)
{
KOKKOS_ASSERT(i_r >= 0 && i_r <= grid.numberSmootherCircles());

Expand Down Expand Up @@ -189,8 +189,7 @@ static KOKKOS_INLINE_FUNCTION void
nodeApplyAscOrthoRadialTake(int i_r, int i_theta, const PolarGrid& grid, bool DirBC_Interior,
HostConstVector<double>& x, HostConstVector<double>& rhs, HostVector<double>& result,
HostConstVector<double>& arr, const HostConstVector<double>& att,
HostConstVector<double>& art, const HostConstVector<double>& detDF,
HostConstVector<double>& coeff_beta)
HostConstVector<double>& art, const HostConstVector<double>& detDF)
{
assert(i_r >= grid.numberSmootherCircles() - 1 && i_r < grid.nr());

Expand Down Expand Up @@ -472,14 +471,11 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackCircleSection(H
const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
const bool DirBC_Interior = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();
HostConstVector<double> coeff_beta = level_cache.coeff_beta();
HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> 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;
Expand All @@ -494,8 +490,7 @@ void ExtrapolatedSmootherTake<LevelCacheType>::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();
Expand All @@ -511,16 +506,12 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteCircleSection(H
const PolarGrid& grid = ExtrapolatedSmoother<LevelCacheType>::grid_;
const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
const bool DirBC_Interior = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;
const int num_omp_threads = ExtrapolatedSmoother<LevelCacheType>::num_omp_threads_;

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();
HostConstVector<double> coeff_beta = level_cache.coeff_beta();
HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> 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;
Expand All @@ -535,8 +526,7 @@ void ExtrapolatedSmootherTake<LevelCacheType>::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();
Expand All @@ -552,16 +542,12 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoBlackRadialSection(H
const PolarGrid& grid = ExtrapolatedSmoother<LevelCacheType>::grid_;
const LevelCacheType& level_cache = ExtrapolatedSmoother<LevelCacheType>::level_cache_;
const bool DirBC_Interior = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;
const int num_omp_threads = ExtrapolatedSmoother<LevelCacheType>::num_omp_threads_;

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();
HostConstVector<double> coeff_beta = level_cache.coeff_beta();
HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();

assert(grid.ntheta() % 2 == 0);
const int start_black_radials = 0;
Expand All @@ -576,8 +562,7 @@ void ExtrapolatedSmootherTake<LevelCacheType>::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();
Expand All @@ -595,14 +580,11 @@ void ExtrapolatedSmootherTake<LevelCacheType>::applyAscOrthoWhiteRadialSection(H
const bool DirBC_Interior = ExtrapolatedSmoother<LevelCacheType>::DirBC_Interior_;
const int num_omp_threads = ExtrapolatedSmoother<LevelCacheType>::num_omp_threads_;

assert(level_cache.cacheDensityProfileCoefficients());
assert(level_cache.cacheDomainGeometry());

HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();
HostConstVector<double> coeff_beta = level_cache.coeff_beta();
HostConstVector<double> arr = level_cache.arr();
HostConstVector<double> att = level_cache.att();
HostConstVector<double> art = level_cache.art();
HostConstVector<double> detDF = level_cache.detDF();

assert(grid.ntheta() % 2 == 0);
const int start_white_radials = 1;
Expand All @@ -617,8 +599,7 @@ void ExtrapolatedSmootherTake<LevelCacheType>::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();
Expand Down
Loading
Loading