Skip to content

Provide density profile coefficients from grid index #313

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

Merged
EmilyBourne merged 7 commits into from
Aug 11, 2026
Merged

Provide density profile coefficients from grid index #313

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ef83903
Remove default constructor
EmilyBourne Aug 10, 2026
3eaf86a
Add PolarGrid element to DensityProfileCoefficients
EmilyBourne Aug 10, 2026
2fb91c9
Update constructor calls
EmilyBourne Aug 10, 2026
96ce97a
Update alpha/beta interface to use indices
EmilyBourne Aug 10, 2026
df48d66
Clang format
EmilyBourne Aug 10, 2026
3f0fdac
Must use global indexing
EmilyBourne Aug 10, 2026
3c27bcd
Pass level depth to LevelCache
EmilyBourne Aug 10, 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: 2 additions & 2 deletions include/GMGPolar/setup.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,14 +42,14 @@ void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()

auto finest_levelCache = std::make_unique<LevelCache<DomainGeometry, DensityProfileCoefficients>>(
*finest_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_,
cache_domain_geometry_);
cache_domain_geometry_, 0);
levels_.emplace_back(0, std::move(finest_grid), std::move(finest_levelCache), extrapolation_, FMG_, PCG_FMG_);

for (int level_depth = 1; level_depth < number_of_levels_; level_depth++) {
auto current_grid = std::make_unique<PolarGrid>(coarseningGrid(levels_[level_depth - 1].grid()));
auto current_levelCache = std::make_unique<LevelCache<DomainGeometry, DensityProfileCoefficients>>(
*current_grid, density_profile_coefficients_, domain_geometry_, cache_density_profile_coefficients_,
cache_domain_geometry_);
cache_domain_geometry_, level_depth);
levels_.emplace_back(level_depth, std::move(current_grid), std::move(current_levelCache), extrapolation_, FMG_,
PCG_FMG_);
}
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/poissonCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class PoissonCoefficients
{
public:
PoissonCoefficients() = default;
explicit PoissonCoefficients(double Rmax, double alpha);
explicit PoissonCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
Comment thread
julianlitz marked this conversation as resolved.
const double Rmax = 1.3;
const double alpha_jump = 0.5 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/sonnendruckerCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class SonnendruckerCoefficients
{
public:
SonnendruckerCoefficients() = default;
explicit SonnendruckerCoefficients(double Rmax, double alpha);
explicit SonnendruckerCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
Comment thread
julianlitz marked this conversation as resolved.
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.66 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/sonnendruckerGyroCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class SonnendruckerGyroCoefficients
{
public:
SonnendruckerGyroCoefficients() = default;
explicit SonnendruckerGyroCoefficients(double Rmax, double alpha);
explicit SonnendruckerGyroCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.66 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/zoniCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class ZoniCoefficients
{
public:
ZoniCoefficients() = default;
explicit ZoniCoefficients(double Rmax, double alpha);
explicit ZoniCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.4837 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/zoniGyroCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class ZoniGyroCoefficients
{
public:
ZoniGyroCoefficients() = default;
explicit ZoniGyroCoefficients(double Rmax, double alpha);
explicit ZoniGyroCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.4837 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/zoniShiftedCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class ZoniShiftedCoefficients
{
public:
ZoniShiftedCoefficients() = default;
explicit ZoniShiftedCoefficients(double Rmax, double alpha);
explicit ZoniShiftedCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.678 * 1.3;
};
Expand Down
9 changes: 5 additions & 4 deletions include/InputFunctions/DensityProfileCoefficients/zoniShiftedGyroCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#pragma once
#include "../../PolarGrid/polargrid.h"

#include <cmath>
#include <Kokkos_Core.hpp>
Expand All @@ -11,15 +12,15 @@ namespace gmgpolar
class ZoniShiftedGyroCoefficients
{
public:
ZoniShiftedGyroCoefficients() = default;
explicit ZoniShiftedGyroCoefficients(double Rmax, double alpha);
explicit ZoniShiftedGyroCoefficients(const PolarGrid& grid, double Rmax, double alpha);

KOKKOS_FUNCTION double alpha(double r, double theta) const;
KOKKOS_FUNCTION double beta(double r, double theta) const;
KOKKOS_FUNCTION double alpha(int i_r, int i_theta) const;
KOKKOS_FUNCTION double beta(int i_r, int i_theta) const;

KOKKOS_FUNCTION double getAlphaJump() const;

private:
PolarGrid grid_;
const double Rmax = 1.3;
const double alpha_jump = 0.678 * 1.3;
};
Expand Down
6 changes: 3 additions & 3 deletions include/InputFunctions/densityProfileCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,9 @@ namespace concepts
{

template <typename T>
concept DensityProfileCoefficients = requires(const T coeffs, double r, double theta) {
{ coeffs.alpha(r, theta) } -> std::convertible_to<double>;
{ coeffs.beta(r, theta) } -> std::convertible_to<double>;
concept DensityProfileCoefficients = requires(const T coeffs, int i_r, int i_theta) {
{ coeffs.alpha(i_r, i_theta) } -> std::convertible_to<double>;
{ coeffs.beta(i_r, i_theta) } -> std::convertible_to<double>;
{ coeffs.getAlphaJump() } -> std::convertible_to<double>;
};

Expand Down
14 changes: 9 additions & 5 deletions include/Level/level.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ class LevelCache
public:
explicit LevelCache(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients,
const DomainGeometry& domain_geometry, const bool cache_density_profile_coefficients,
const bool cache_domain_geometry);
const bool cache_domain_geometry, const int level_depth);

const DomainGeometry& domainGeometry() const;
const DensityProfileCoefficients& densityProfileCoefficients() const;
Expand All @@ -155,8 +155,10 @@ class LevelCache
double theta, double& coeff_beta, double& arr, double& att, double& art,
double& detDF) const
{
coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index]
: density_profile_coefficients_.beta(r, theta);
const int i_r_glob = i_r << level_depth_;
const int i_theta_glob = i_theta << level_depth_;
coeff_beta = cache_density_profile_coefficients_ ? coeff_beta_[global_index]
: density_profile_coefficients_.beta(i_r_glob, i_theta_glob);

if (cache_domain_geometry_) {
arr = arr_[global_index];
Expand All @@ -165,14 +167,16 @@ class LevelCache
detDF = detDF_[global_index];
}
else {
double coeff_alpha = cache_density_profile_coefficients_ ? coeff_alpha_[global_index]
: density_profile_coefficients_.alpha(r, theta);
double coeff_alpha = cache_density_profile_coefficients_
? coeff_alpha_[global_index]
: density_profile_coefficients_.alpha(i_r_glob, i_theta_glob);

compute_jacobian_elements(domain_geometry_, r, theta, coeff_alpha, arr, att, art, detDF);
}
}

private:
const int level_depth_;
const DomainGeometry domain_geometry_;
const DensityProfileCoefficients density_profile_coefficients_;

Expand Down
31 changes: 18 additions & 13 deletions include/Level/levelCache.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ template <concepts::DensityProfileCoefficients DensityProfileCoefficients>
static void cache_density_profile_coefficients(const PolarGrid& grid,
const DensityProfileCoefficients& density_profile_coefficients,
const Vector<double>& coeff_alpha, const Vector<double>& coeff_beta,
const bool cache_domain_geometry)
const bool cache_domain_geometry, int level_depth)
{
Kokkos::parallel_for(
"Cache density profile coefficients",
Expand All @@ -16,21 +16,21 @@ static void cache_density_profile_coefficients(const PolarGrid& grid,
),
// Kokkos lambda function to execute for each point in the index space
KOKKOS_LAMBDA(const int i_r, const int i_theta) {
const double r = grid.radius(i_r);
const double theta = grid.theta(i_theta);
const int index = grid.index(i_r, i_theta);
const int index = grid.index(i_r, i_theta);
const int i_r_glob = i_r << level_depth;
const int i_theta_glob = i_theta << level_depth;
if (!cache_domain_geometry) {
coeff_alpha(index) = density_profile_coefficients.alpha(r, theta);
coeff_alpha(index) = density_profile_coefficients.alpha(i_r_glob, i_theta_glob);
}
coeff_beta(index) = density_profile_coefficients.beta(r, theta);
coeff_beta(index) = density_profile_coefficients.beta(i_r_glob, i_theta_glob);
});
}

template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
static void cache_domain_geometry(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients,
const DomainGeometry& domain_geometry, const Vector<double>& vec_arr,
const Vector<double>& vec_att, const Vector<double>& vec_art,
const Vector<double>& vec_detDF)
const Vector<double>& vec_detDF, int level_depth)
{
// We split the loops into two regions to better respect the
// access patterns of the smoother and improve cache locality
Expand All @@ -47,7 +47,9 @@ static void cache_domain_geometry(const PolarGrid& grid, const DensityProfileCoe
const double r = grid.radius(i_r);
const double theta = grid.theta(i_theta);
const int index = grid.index(i_r, i_theta);
const double coeff_alpha = density_profile_coefficients.alpha(r, theta);
const int i_r_glob = i_r << level_depth;
const int i_theta_glob = i_theta << level_depth;
const double coeff_alpha = density_profile_coefficients.alpha(i_r_glob, i_theta_glob);

double arr, att, art, detDF;
compute_jacobian_elements(domain_geometry, r, theta, coeff_alpha, arr, att, art, detDF);
Expand All @@ -68,7 +70,9 @@ static void cache_domain_geometry(const PolarGrid& grid, const DensityProfileCoe
const double theta = grid.theta(i_theta);
const double r = grid.radius(i_r);
const int index = grid.index(i_r, i_theta);
const double coeff_alpha = density_profile_coefficients.alpha(r, theta);
const int i_r_glob = i_r << level_depth;
const int i_theta_glob = i_theta << level_depth;
const double coeff_alpha = density_profile_coefficients.alpha(i_r_glob, i_theta_glob);

double arr, att, art, detDF;
compute_jacobian_elements(domain_geometry, r, theta, coeff_alpha, arr, att, art, detDF);
Expand All @@ -85,8 +89,9 @@ template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoeff
LevelCache<DomainGeometry, DensityProfileCoefficients>::LevelCache(
const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients,
const DomainGeometry& domain_geometry, const bool cache_density_profile_coefficients,
const bool cache_domain_geometry)
: domain_geometry_(domain_geometry)
const bool cache_domain_geometry, const int level_depth)
: level_depth_(level_depth)
, domain_geometry_(domain_geometry)
, density_profile_coefficients_(density_profile_coefficients)
, cache_density_profile_coefficients_(cache_density_profile_coefficients)
// If the domain geometry is cached, we don't need to cache the alpha coefficient
Expand All @@ -103,14 +108,14 @@ LevelCache<DomainGeometry, DensityProfileCoefficients>::LevelCache(
// repeated expensive evaluations during runtime computations
if (cache_density_profile_coefficients_) {
level_cache_helpers::cache_density_profile_coefficients(grid, density_profile_coefficients, coeff_alpha_,
coeff_beta_, cache_domain_geometry);
coeff_beta_, cache_domain_geometry, level_depth);
}

// Pre-compute and store Jacobian matrix elements (arr, att, art, detDF) at all grid nodes
// to avoid repeated coordinate transformation calculations during domain operations
if (cache_domain_geometry_) {
level_cache_helpers::cache_domain_geometry(grid, density_profile_coefficients, domain_geometry, arr_, att_,
art_, detDF_);
art_, detDF_, level_depth);
}
Kokkos::fence();
}
Expand Down
18 changes: 9 additions & 9 deletions src/ConfigParser/select_test_case.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,18 +65,18 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
switch (alpha_type) {
case AlphaCoeff::POISSON:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(PoissonCoefficients(Rmax, alpha_jump));
std::make_unique<DensityProfileCoefficientsVariant>(PoissonCoefficients(grid_, Rmax, alpha_jump));
break;

case AlphaCoeff::SONNENDRUCKER:
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerCoefficients(Rmax, alpha_jump));
std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerCoefficients(grid_, Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerGyroCoefficients(Rmax, alpha_jump));
density_profile_coefficients_ = std::make_unique<DensityProfileCoefficientsVariant>(
SonnendruckerGyroCoefficients(grid_, Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand All @@ -87,11 +87,11 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniCoefficients(Rmax, alpha_jump));
std::make_unique<DensityProfileCoefficientsVariant>(ZoniCoefficients(grid_, Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniGyroCoefficients(Rmax, alpha_jump));
std::make_unique<DensityProfileCoefficientsVariant>(ZoniGyroCoefficients(grid_, Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand All @@ -102,11 +102,11 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedCoefficients(Rmax, alpha_jump));
std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedCoefficients(grid_, Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedGyroCoefficients(Rmax, alpha_jump));
density_profile_coefficients_ = std::make_unique<DensityProfileCoefficientsVariant>(
ZoniShiftedGyroCoefficients(grid_, Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand Down
Loading
Loading