Skip to content

Template GMGPolar class on DensityProfileCoefficients #137

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 27 commits into from
Jan 27, 2026
Merged

Template GMGPolar class on DensityProfileCoefficients #137

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
0defbc9
Add a concept describing a DomainGeometry
EmilyBourne Dec 15, 2025
45e599a
Template GMGPolar on DomainGeometry
EmilyBourne Dec 15, 2025
880415b
Add templating on class methods
EmilyBourne Dec 15, 2025
fd3152c
Modify concept so abstract class does not match
EmilyBourne Dec 17, 2025
7964abd
Clang format
EmilyBourne Dec 17, 2025
9b230da
Use inheritance for GMGPolar to allow main declaration with ConfigParser
EmilyBourne Dec 17, 2025
ae1aee5
Add a method to get a std::unique_ptr<IGMGPolar> from the ConfigParser
EmilyBourne Dec 17, 2025
bf8cc11
Clang format
EmilyBourne Dec 17, 2025
880d87b
Use more common naming strategy
EmilyBourne Dec 17, 2025
60a8222
Commit missing file
EmilyBourne Dec 17, 2025
c27dd92
Add DensityProfileCoefficients to template parameters of GMGPolar in …
EmilyBourne Dec 17, 2025
cda5211
Update class prototypes
EmilyBourne Dec 17, 2025
e1a5734
Ensure DensityProfileCoefficients type is available
EmilyBourne Dec 17, 2025
6d0685a
Clang formatting
EmilyBourne Dec 17, 2025
facae26
Clean up includes
EmilyBourne Dec 17, 2025
77dbbdd
Improve coverage
EmilyBourne Dec 18, 2025
166bc2d
Only check values that are set
EmilyBourne Dec 18, 2025
93ffb92
Merge remote-tracking branch 'GMGPolar/main' into ebourne_gmgpolar_te…
EmilyBourne Dec 20, 2025
a50ee04
Add explanatory comments
EmilyBourne Dec 20, 2025
1ca2ab5
Create variable and add comment
EmilyBourne Dec 20, 2025
e00f038
Clang format
EmilyBourne Dec 20, 2025
7f9f50e
Merge remote-tracking branch 'GMGPolar/main' into ebourne_gmgpolar_te…
EmilyBourne Jan 26, 2026
a2c5c31
Merge branch 'ebourne_gmgpolar_template' into ebourne_gmgpolar_templa…
EmilyBourne Jan 26, 2026
d1fcff7
Fix merge
EmilyBourne Jan 26, 2026
2377f33
Merge remote-tracking branch 'GMGPolar/main' into ebourne_gmgpolar_te…
EmilyBourne Jan 27, 2026
f30bd4a
Merge branch 'ebourne_gmgpolar_template' into ebourne_gmgpolar_templa…
EmilyBourne Jan 27, 2026
266dbb6
Merge remote-tracking branch 'GMGPolar/main' into ebourne_gmgpolar_te…
EmilyBourne Jan 27, 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/ConfigParser/config_parser.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ class ConfigParser

// Test Case
const DomainGeometryVariant& domainGeometry() const;
const DensityProfileCoefficients& densityProfileCoefficients() const;
const DensityProfileCoefficientsVariant& densityProfileCoefficients() const;
const BoundaryConditions& boundaryConditions() const;
const SourceTerm& sourceTerm() const;
const ExactSolution& exactSolution() const;
Expand Down Expand Up @@ -61,7 +61,7 @@ class ConfigParser
cmdline::parser parser_;
// Input Functions
std::unique_ptr<const DomainGeometryVariant> domain_geometry_;
std::unique_ptr<const DensityProfileCoefficients> density_profile_coefficients_;
std::unique_ptr<const DensityProfileCoefficientsVariant> density_profile_coefficients_;
std::unique_ptr<const BoundaryConditions> boundary_conditions_;
std::unique_ptr<const SourceTerm> source_term_;
std::unique_ptr<const ExactSolution> exact_solution_;
Expand Down
4 changes: 4 additions & 0 deletions include/ConfigParser/test_selection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,7 @@
#include "../../include/GMGPolar/test_cases.h"

using DomainGeometryVariant = std::variant<CircularGeometry, ShafranovGeometry, CzarnyGeometry, CulhamGeometry>;

using DensityProfileCoefficientsVariant =
std::variant<PoissonCoefficients, SonnendruckerCoefficients, SonnendruckerGyroCoefficients, ZoniCoefficients,
ZoniGyroCoefficients, ZoniShiftedCoefficients, ZoniShiftedGyroCoefficients>;
4 changes: 2 additions & 2 deletions include/GMGPolar/build_rhs_f.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#include "../../include/GMGPolar/gmgpolar.h"

template <concepts::DomainGeometry DomainGeometry>
void GMGPolar<DomainGeometry>::discretize_rhs_f(const Level& level, Vector<double> rhs_f)
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::discretize_rhs_f(const Level& level, Vector<double> rhs_f)
{
const PolarGrid& grid = level.grid();
assert(std::ssize(rhs_f) == grid.numberOfNodes());
Expand Down
23 changes: 5 additions & 18 deletions include/GMGPolar/gmgpolar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,31 +1,16 @@
#pragma once

#include <chrono>
#include <filesystem>
#include <iostream>
#include <memory>
#include <omp.h>
#include <optional>
#include <utility>

class Level;
class LevelCache;

#include "../InputFunctions/boundaryConditions.h"
#include "../InputFunctions/densityProfileCoefficients.h"
#include "../InputFunctions/domainGeometry.h"
#include "../InputFunctions/exactSolution.h"
#include "../InputFunctions/sourceTerm.h"
#include "../Interpolation/interpolation.h"
#include "../Level/level.h"
#include "../LinearAlgebra/vector.h"
#include "../LinearAlgebra/vector_operations.h"
#include "../PolarGrid/polargrid.h"
#include "../common/global_definitions.h"
#include "test_cases.h"

#include "igmgpolar.h"

template <concepts::DomainGeometry DomainGeometry>
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
class GMGPolar : public IGMGPolar
{
public:
Expand All @@ -41,8 +26,9 @@ class GMGPolar : public IGMGPolar
// - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE.
GMGPolar(const PolarGrid& grid, const DomainGeometry& domain_geometry,
const DensityProfileCoefficients& density_profile_coefficients)
: IGMGPolar(grid, density_profile_coefficients)
: IGMGPolar(grid)
, domain_geometry_(domain_geometry)
, density_profile_coefficients_(density_profile_coefficients)
{
}

Expand All @@ -57,6 +43,7 @@ class GMGPolar : public IGMGPolar
/* Grid Configuration & Input Functions */
/* ------------------------------------ */
const DomainGeometry& domain_geometry_;
const DensityProfileCoefficients& density_profile_coefficients_;

/* --------------- */
/* Setup Functions */
Expand Down
4 changes: 1 addition & 3 deletions include/GMGPolar/igmgpolar.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,7 @@ class IGMGPolar
// with Dirichlet boundary condition u = u_D on \partial \Omega.
// Parameters:
// - grid: Cartesian mesh discretizing the computational domain.
// - density_profile_coefficients: Coefficients \alpha and \beta defining the PDE.
IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients);
IGMGPolar(const PolarGrid& grid);

/* ---------------------------------------------------------------------- */
/* General output & visualization */
Expand Down Expand Up @@ -190,7 +189,6 @@ class IGMGPolar
/* Grid Configuration & Input Functions */
/* ------------------------------------ */
const PolarGrid& grid_;
const DensityProfileCoefficients& density_profile_coefficients_;
const ExactSolution* exact_solution_; // Optional exact solution for validation

/* ------------------ */
Expand Down
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
@@ -1,6 +1,6 @@

template <concepts::DomainGeometry DomainGeometry>
void GMGPolar<DomainGeometry>::setup()
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::setup()
{
LIKWID_START("Setup");
auto start_setup = std::chrono::high_resolution_clock::now();
Expand Down
12 changes: 7 additions & 5 deletions include/GMGPolar/writeToVTK.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
#include "../../include/GMGPolar/gmgpolar.h"

template <concepts::DomainGeometry DomainGeometry>
void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path, const PolarGrid& grid)
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std::filesystem::path& file_path,
const PolarGrid& grid)
{
const auto filename = file_path.stem().string() + ".vtu";

Expand Down Expand Up @@ -59,9 +60,10 @@ void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path
<< "</VTKFile>\n";
}

template <concepts::DomainGeometry DomainGeometry>
void GMGPolar<DomainGeometry>::writeToVTK(const std::filesystem::path& file_path, const Level& level,
ConstVector<double> grid_function)
template <concepts::DomainGeometry DomainGeometry, concepts::DensityProfileCoefficients DensityProfileCoefficients>
void GMGPolar<DomainGeometry, DensityProfileCoefficients>::writeToVTK(const std::filesystem::path& file_path,
const Level& level,
ConstVector<double> grid_function)
{
const PolarGrid& grid = level.grid();
const LevelCache& level_cache = level.levelCache();
Expand Down
15 changes: 14 additions & 1 deletion include/InputFunctions/densityProfileCoefficients.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,17 @@ class DensityProfileCoefficients

// Only used in custom mesh generation -> refinement_radius
virtual double getAlphaJump() const = 0;
};
};

namespace concepts
{

template <typename T>
concept DensityProfileCoefficients =
!std::same_as<T, 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>;
{ coeffs.getAlphaJump() } -> std::convertible_to<double>;
};

} // namespace concepts
2 changes: 1 addition & 1 deletion src/ConfigParser/config_parser.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,7 @@ const DomainGeometryVariant& ConfigParser::domainGeometry() const
return *domain_geometry_.get();
}

const DensityProfileCoefficients& ConfigParser::densityProfileCoefficients() const
const DensityProfileCoefficientsVariant& ConfigParser::densityProfileCoefficients() const
{
return *density_profile_coefficients_.get();
}
Expand Down
37 changes: 23 additions & 14 deletions src/ConfigParser/select_test_case.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,21 +5,23 @@ std::unique_ptr<IGMGPolar> ConfigParser::solver() const
{
// Create local aliases so the class doesn't need to be captured by the lamda
// These are references, not copies.
const PolarGrid& grid = grid_;
const DensityProfileCoefficients& density_profile_coefficients = *density_profile_coefficients_;
const PolarGrid& grid = grid_;

// Create a solver specialized to the active domain geometry.
return std::visit(
[&grid, &density_profile_coefficients](auto const& domain_geometry) {
// Deduce the concrete geometry type
using DomainGeomType = std::decay_t<decltype(domain_geometry)>;
[&grid](auto const& domain_geometry, auto const& density_profile_coefficients) {
using DomainGeomType = std::decay_t<decltype(domain_geometry)>;
using DensityProfileCoefficientsType = std::decay_t<decltype(density_profile_coefficients)>;

// Construct the solver specialized for this geometry type.
std::unique_ptr<IGMGPolar> solver =
std::make_unique<GMGPolar<DomainGeomType>>(grid, domain_geometry, density_profile_coefficients);
std::make_unique<GMGPolar<DomainGeomType, DensityProfileCoefficientsType>>(
grid, domain_geometry, density_profile_coefficients);

// The lambdas must return objects of identical type
return solver;
},
*domain_geometry_);
*domain_geometry_, *density_profile_coefficients_);
}

void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType problem_type, AlphaCoeff alpha_type,
Expand Down Expand Up @@ -52,16 +54,19 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
/* Density Profile Coefficients */
switch (alpha_type) {
case AlphaCoeff::POISSON:
density_profile_coefficients_ = std::make_unique<PoissonCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(PoissonCoefficients(Rmax, alpha_jump));
break;

case AlphaCoeff::SONNENDRUCKER:
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ = std::make_unique<SonnendruckerCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerCoefficients(Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ = std::make_unique<SonnendruckerGyroCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(SonnendruckerGyroCoefficients(Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand All @@ -71,10 +76,12 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
case AlphaCoeff::ZONI:
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ = std::make_unique<ZoniCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniCoefficients(Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ = std::make_unique<ZoniGyroCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniGyroCoefficients(Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand All @@ -84,10 +91,12 @@ void ConfigParser::selectTestCase(GeometryType geometry_type, ProblemType proble
case AlphaCoeff::ZONI_SHIFTED:
switch (beta_type) {
case BetaCoeff::ZERO:
density_profile_coefficients_ = std::make_unique<ZoniShiftedCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedCoefficients(Rmax, alpha_jump));
break;
case BetaCoeff::ALPHA_INVERSE:
density_profile_coefficients_ = std::make_unique<ZoniShiftedGyroCoefficients>(Rmax, alpha_jump);
density_profile_coefficients_ =
std::make_unique<DensityProfileCoefficientsVariant>(ZoniShiftedGyroCoefficients(Rmax, alpha_jump));
break;
default:
throw std::runtime_error("Invalid beta.\n");
Expand Down
3 changes: 1 addition & 2 deletions src/GMGPolar/gmgpolar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,8 @@
/* ---------------------------------------------------------------------- */
/* Constructor & Initialization */
/* ---------------------------------------------------------------------- */
IGMGPolar::IGMGPolar(const PolarGrid& grid, const DensityProfileCoefficients& density_profile_coefficients)
IGMGPolar::IGMGPolar(const PolarGrid& grid)
: grid_(grid)
, density_profile_coefficients_(density_profile_coefficients)
, exact_solution_(nullptr)
// General solver output and visualization settings
, verbose_(0)
Expand Down
4 changes: 4 additions & 0 deletions tests/ConfigParser/config_parser.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,10 @@ TEST_P(ConfigParserTest, ParseAllGeometryAndProblemCombinations)
EXPECT_DOUBLE_EQ(parser.absoluteTolerance().value(), absoluteTolerance);
ASSERT_TRUE(parser.relativeTolerance().has_value());
EXPECT_DOUBLE_EQ(parser.relativeTolerance().value(), relativeTolerance);

// Solver
std::unique_ptr<IGMGPolar> solver = parser.solver();
EXPECT_EQ(&solver->grid(), &parser.grid());
}

// Define test cases covering all combinations
Expand Down
10 changes: 5 additions & 5 deletions tests/GMGPolar/convergence_order.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,11 +152,11 @@ std::vector<double> refine(std::vector<double> const& original_points)
return refined;
}

std::tuple<double, double> get_gmgpolar_error(PolarGrid const& grid, CzarnyGeometry const& domain_geometry,
DensityProfileCoefficients const& coefficients,
BoundaryConditions const& boundary_conditions,
SourceTerm const& source_term, ExactSolution const& solution,
ExtrapolationType extrapolation)
template <class DensityProfileCoefficients>
std::tuple<double, double>
get_gmgpolar_error(PolarGrid const& grid, CzarnyGeometry const& domain_geometry,
DensityProfileCoefficients const& coefficients, BoundaryConditions const& boundary_conditions,
SourceTerm const& source_term, ExactSolution const& solution, ExtrapolationType extrapolation)
{
GMGPolar gmgpolar(grid, domain_geometry, coefficients);
gmgpolar.setSolution(&solution);
Expand Down