Skip to content

Template GMGPolar class on DomainGeometry #136

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

Open
EmilyBourne wants to merge 14 commits into
base: main
Choose a base branch
from
Open

Template GMGPolar class on DomainGeometry #136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 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
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
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
7 changes: 5 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 @@ -9,6 +9,8 @@
#include "../../include/common/global_definitions.h"
#include "../../include/PolarGrid/polargrid.h"
#include "../../include/GMGPolar/test_cases.h"
#include "../../include/GMGPolar/igmgpolar.h"
#include "test_selection.h"

class ConfigParser
{
Expand All @@ -22,11 +24,12 @@ class ConfigParser
bool parse(int argc, char* argv[]);

// Test Case
const DomainGeometry& domainGeometry() const;
const DomainGeometryVariant& domainGeometry() const;
const DensityProfileCoefficients& densityProfileCoefficients() const;
const BoundaryConditions& boundaryConditions() const;
const SourceTerm& sourceTerm() const;
const ExactSolution& exactSolution() const;
std::unique_ptr<IGMGPolar> solver() const;

// Control Parameters
int verbose() const;
Expand Down Expand Up @@ -58,7 +61,7 @@ class ConfigParser
// Parse command-line arguments to extract problem configuration
cmdline::parser parser_;
// Input Functions
std::unique_ptr<const DomainGeometry> domain_geometry_;
std::unique_ptr<const DomainGeometryVariant> domain_geometry_;
std::unique_ptr<const DensityProfileCoefficients> density_profile_coefficients_;
std::unique_ptr<const BoundaryConditions> boundary_conditions_;
std::unique_ptr<const SourceTerm> source_term_;
Expand Down
5 changes: 5 additions & 0 deletions include/ConfigParser/test_selection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
#pragma once
#include <variant>
#include "../../include/GMGPolar/test_cases.h"

using DomainGeometryVariant = std::variant<CircularGeometry, ShafranovGeometry, CzarnyGeometry, CulhamGeometry>;
141 changes: 141 additions & 0 deletions include/GMGPolar/build_rhs_f.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,141 @@
#include "../../include/GMGPolar/gmgpolar.h"

template <concepts::DomainGeometry DomainGeometry>
void GMGPolar<DomainGeometry>::discretize_rhs_f(const Level& level, Vector<double> rhs_f)
{
const PolarGrid& grid = level.grid();
assert(rhs_f.size() == static_cast<uint>(grid.numberOfNodes()));

if (level.levelCache().cacheDomainGeometry()) {
/* DomainGeometry is cached */
const auto& detDF_cache = level.levelCache().detDF();
#pragma omp parallel
{
// ---------------------------------------------- //
// Discretize rhs values (circular index section) //
// ---------------------------------------------- //
#pragma omp for nowait
for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) {
double r = grid.radius(i_r);
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
double theta = grid.theta(i_theta);
if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) {
double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
const double detDF = detDF_cache[grid.index(i_r, i_theta)];
rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF);
}
else if (i_r == 0 && DirBC_Interior_) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
else if (i_r == grid.nr() - 1) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
}
}

// -------------------------------------------- //
// Discretize rhs values (radial index section) //
// -------------------------------------------- //
#pragma omp for nowait
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
double theta = grid.theta(i_theta);
for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {
double r = grid.radius(i_r);
if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) {
double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
const double detDF = detDF_cache[grid.index(i_r, i_theta)];
rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF);
}
else if (i_r == 0 && DirBC_Interior_) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
else if (i_r == grid.nr() - 1) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
}
}
}
}
else {
/* DomainGeometry is not cached */

#pragma omp parallel
{
// ---------------------------------------------- //
// Discretize rhs values (circular index section) //
// ---------------------------------------------- //
#pragma omp for nowait
for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) {
double r = grid.radius(i_r);
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
double theta = grid.theta(i_theta);

if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) {
double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
/* Calculate the elements of the Jacobian matrix for the transformation mapping */
/* The Jacobian matrix is: */
/* [Jrr, Jrt] */
/* [Jtr, Jtt] */
double Jrr = domain_geometry_.dFx_dr(r, theta);
double Jtr = domain_geometry_.dFy_dr(r, theta);
double Jrt = domain_geometry_.dFx_dt(r, theta);
double Jtt = domain_geometry_.dFy_dt(r, theta);
/* Compute the determinant of the Jacobian matrix */
double detDF = Jrr * Jtt - Jrt * Jtr;
rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF);
}
else if (i_r == 0 && DirBC_Interior_) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
else if (i_r == grid.nr() - 1) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
}
}

// -------------------------------------------- //
// Discretize rhs values (radial index section) //
// -------------------------------------------- //
#pragma omp for nowait
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
double theta = grid.theta(i_theta);

for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {
double r = grid.radius(i_r);
if ((0 < i_r && i_r < grid.nr() - 1) || (i_r == 0 && !DirBC_Interior_)) {
double h1 = (i_r == 0) ? 2.0 * grid.radius(0) : grid.radialSpacing(i_r - 1);
double h2 = grid.radialSpacing(i_r);
double k1 = grid.angularSpacing(i_theta - 1);
double k2 = grid.angularSpacing(i_theta);
/* Calculate the elements of the Jacobian matrix for the transformation mapping */
/* The Jacobian matrix is: */
/* [Jrr, Jrt] */
/* [Jtr, Jtt] */
double Jrr = domain_geometry_.dFx_dr(r, theta);
double Jtr = domain_geometry_.dFy_dr(r, theta);
double Jrt = domain_geometry_.dFx_dt(r, theta);
double Jtt = domain_geometry_.dFy_dt(r, theta);
/* Compute the determinant of the Jacobian matrix */
double detDF = Jrr * Jtt - Jrt * Jtr;
rhs_f[grid.index(i_r, i_theta)] *= 0.25 * (h1 + h2) * (k1 + k2) * fabs(detDF);
}
else if (i_r == 0 && DirBC_Interior_) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
else if (i_r == grid.nr() - 1) {
rhs_f[grid.index(i_r, i_theta)] *= 1.0;
}
}
}
}
}
}
Loading