Skip to content

Commit d86c55a

Browse files
EmilyBournejulianlitz
EmilyBourne
and
julianlitz
authored
Template Level on DomainGeometry (#190)
Co-authored-by: julianlitz <julian.litz@icloud.com> Co-authored-by: Julian Litz <91479202+julianlitz@users.noreply.github.com>
1 parent b3b9c52 commit d86c55a

105 files changed

Lines changed: 4966 additions & 4236 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎include/Definitions/geometry_helper.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55

66
#include <cmath>
77

8+
template <concepts::DomainGeometry DomainGeometry>
89
inline void compute_jacobian_elements(const DomainGeometry& domain_geometry, double r, double theta, double coeff_alpha,
910
double& arr, double& att, double& art, double& detDF)
1011
{

‎include/DirectSolver/DirectSolver-COO-MUMPS-Give/applySymmetryShift.inl‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
#pragma once
2+
3+
#ifdef GMGPOLAR_USE_MUMPS
4+
5+
/* ----------------------- */
6+
/* Boundary Symmetry Shift */
7+
/* ----------------------- */
8+
9+
template <concepts::DomainGeometry DomainGeometry>
10+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::applySymmetryShiftInnerBoundary(Vector<double> x) const
11+
{
12+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
13+
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
14+
15+
assert(DirectSolver<DomainGeometry>::DirBC_Interior_);
16+
17+
int i_r;
18+
double r;
19+
int global_index;
20+
double h1, h2, k1, k2;
21+
double coeff1, coeff2;
22+
double coeff_beta, arr, att, art, detDF;
23+
24+
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
25+
const double theta = grid.theta(i_theta);
26+
/* -------------------------- */
27+
/* Node on the inner boundary */
28+
/* -------------------------- */
29+
i_r = 0;
30+
r = grid.radius(i_r);
31+
global_index = grid.index(i_r, i_theta);
32+
33+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
34+
35+
h2 = grid.radialSpacing(i_r);
36+
k1 = grid.angularSpacing(i_theta - 1);
37+
k2 = grid.angularSpacing(i_theta);
38+
39+
coeff2 = 0.5 * (k1 + k2) / h2;
40+
41+
/* Fill x(i+1,j) */
42+
x(grid.index(i_r + 1, i_theta)) -= -coeff2 * arr * x(grid.index(i_r, i_theta)) /* Left */
43+
+ 0.25 * art * x(grid.index(i_r, i_theta + 1)) /* Top Left */
44+
- 0.25 * art * x(grid.index(i_r, i_theta - 1)); /* Bottom Left */
45+
46+
/* --------------------------- */
47+
/* Node next to inner boundary */
48+
/* --------------------------- */
49+
i_r = 1;
50+
r = grid.radius(i_r);
51+
global_index = grid.index(i_r, i_theta);
52+
53+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
54+
55+
h1 = grid.radialSpacing(i_r - 1);
56+
k1 = grid.angularSpacing(i_theta - 1);
57+
k2 = grid.angularSpacing(i_theta);
58+
59+
coeff1 = 0.5 * (k1 + k2) / h1;
60+
61+
/* Fill x(i,j) */
62+
x(grid.index(i_r, i_theta)) -= -coeff1 * arr * x(grid.index(i_r - 1, i_theta)); /* Left */
63+
/* Fill x(i,j-1) */
64+
x(grid.index(i_r, i_theta - 1)) -= +0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
65+
/* Fill x(i,j+1) */
66+
x(grid.index(i_r, i_theta + 1)) -= -0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
67+
}
68+
}
69+
70+
template <concepts::DomainGeometry DomainGeometry>
71+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::applySymmetryShiftOuterBoundary(Vector<double> x) const
72+
{
73+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
74+
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
75+
76+
int i_r;
77+
double r;
78+
int global_index;
79+
double h1, h2, k1, k2;
80+
double coeff1, coeff2;
81+
double coeff_beta, arr, att, art, detDF;
82+
83+
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
84+
const double theta = grid.theta(i_theta);
85+
/* --------------------------- */
86+
/* Node next to outer boundary */
87+
/* --------------------------- */
88+
i_r = grid.nr() - 2;
89+
r = grid.radius(i_r);
90+
global_index = grid.index(i_r, i_theta);
91+
92+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
93+
94+
h2 = grid.radialSpacing(i_r);
95+
k1 = grid.angularSpacing(i_theta - 1);
96+
k2 = grid.angularSpacing(i_theta);
97+
98+
coeff2 = 0.5 * (k1 + k2) / h2;
99+
100+
/* Fill result(i,j) */
101+
x(grid.index(i_r, i_theta)) -= -coeff2 * arr * x(grid.index(i_r + 1, i_theta)); /* Right */
102+
/* Fill result(i,j-1) */
103+
x(grid.index(i_r, i_theta - 1)) -= -0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
104+
/* Fill result(i,j+1) */
105+
x(grid.index(i_r, i_theta + 1)) -= +0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
106+
107+
/* -------------------------- */
108+
/* Node on the outer boundary */
109+
/* -------------------------- */
110+
i_r = grid.nr() - 1;
111+
r = grid.radius(i_r);
112+
global_index = grid.index(i_r, i_theta);
113+
114+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
115+
116+
h1 = grid.radialSpacing(i_r - 1);
117+
k1 = grid.angularSpacing(i_theta - 1);
118+
k2 = grid.angularSpacing(i_theta);
119+
120+
coeff1 = 0.5 * (k1 + k2) / h1;
121+
122+
/* Fill result(i-1,j) */
123+
x(grid.index(i_r - 1, i_theta)) -= -coeff1 * arr * x(grid.index(i_r, i_theta)) /* Right */
124+
- 0.25 * art * x(grid.index(i_r, i_theta + 1)) /* Top Right */
125+
+ 0.25 * art * x(grid.index(i_r, i_theta - 1)); /* Bottom Right */
126+
}
127+
}
128+
129+
template <concepts::DomainGeometry DomainGeometry>
130+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::applySymmetryShift(Vector<double> x) const
131+
{
132+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
133+
const bool DirBC_Interior = DirectSolver<DomainGeometry>::DirBC_Interior_;
134+
const int num_omp_threads = DirectSolver<DomainGeometry>::num_omp_threads_;
135+
136+
assert(std::ssize(x) == grid.numberOfNodes());
137+
assert(grid.nr() >= 4);
138+
139+
#pragma omp parallel sections num_threads(num_omp_threads)
140+
{
141+
#pragma omp section
142+
{
143+
if (DirBC_Interior) {
144+
applySymmetryShiftInnerBoundary(x);
145+
}
146+
}
147+
148+
#pragma omp section
149+
{
150+
applySymmetryShiftOuterBoundary(x);
151+
}
152+
}
153+
}
154+
155+
#endif

src/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.cpp renamed to include/DirectSolver/DirectSolver-COO-MUMPS-Give/buildSolverMatrix.inl

Lines changed: 57 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,10 @@
1-
#include "../../../include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h"
2-
3-
#include "../../../include/Definitions/geometry_helper.h"
1+
#pragma once
42

53
#ifdef GMGPOLAR_USE_MUMPS
64

5+
namespace direct_solver_coo_mumps_give
6+
{
7+
78
static inline void updateMatrixElement(SparseMatrixCOO<double>& matrix, int ptr, int offset, int row, int col,
89
double val)
910
{
@@ -12,11 +13,16 @@ static inline void updateMatrixElement(SparseMatrixCOO<double>& matrix, int ptr,
1213
matrix.value(ptr + offset) += val;
1314
}
1415

15-
void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid,
16-
bool DirBC_Interior, SparseMatrixCOO<double>& solver_matrix,
17-
double arr, double att, double art, double detDF,
18-
double coeff_beta)
16+
} // namespace direct_solver_coo_mumps_give
17+
18+
template <concepts::DomainGeometry DomainGeometry>
19+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::nodeBuildSolverMatrixGive(int i_r, int i_theta, const PolarGrid& grid,
20+
bool DirBC_Interior,
21+
SparseMatrixCOO<double>& solver_matrix,
22+
double arr, double att, double art,
23+
double detDF, double coeff_beta)
1924
{
25+
using direct_solver_coo_mumps_give::updateMatrixElement;
2026
int ptr, offset;
2127
int row, col;
2228
double val;
@@ -251,7 +257,7 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int i_r, int i_theta
251257
const int i_theta_M1 = grid.wrapThetaIndex(i_theta - 1);
252258
const int i_theta_P1 = grid.wrapThetaIndex(i_theta + 1);
253259

254-
assert(grid_.ntheta() % 2 == 0);
260+
assert(grid.ntheta() % 2 == 0);
255261
const int i_theta_AcrossOrigin = grid.wrapThetaIndex(i_theta + grid.ntheta() / 2);
256262

257263
double h1 = 2.0 * grid.radius(0);
@@ -773,87 +779,100 @@ void DirectSolver_COO_MUMPS_Give::nodeBuildSolverMatrixGive(int i_r, int i_theta
773779
}
774780
}
775781

776-
void DirectSolver_COO_MUMPS_Give::buildSolverMatrixCircleSection(const int i_r, SparseMatrixCOO<double>& solver_matrix)
782+
template <concepts::DomainGeometry DomainGeometry>
783+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixCircleSection(const int i_r,
784+
SparseMatrixCOO<double>& solver_matrix)
777785
{
778-
const double r = grid_.radius(i_r);
779-
for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) {
780-
const int global_index = grid_.index(i_r, i_theta);
781-
const double theta = grid_.theta(i_theta);
786+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
787+
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
788+
789+
const double r = grid.radius(i_r);
790+
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
791+
const int global_index = grid.index(i_r, i_theta);
792+
const double theta = grid.theta(i_theta);
782793

783794
double coeff_beta, arr, att, art, detDF;
784-
level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
795+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
785796

786797
// Build solver matrix at the current node
787-
nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF,
788-
coeff_beta);
798+
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirectSolver<DomainGeometry>::DirBC_Interior_, solver_matrix, arr,
799+
att, art, detDF, coeff_beta);
789800
}
790801
}
791802

792-
void DirectSolver_COO_MUMPS_Give::buildSolverMatrixRadialSection(const int i_theta,
793-
SparseMatrixCOO<double>& solver_matrix)
803+
template <concepts::DomainGeometry DomainGeometry>
804+
void DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrixRadialSection(const int i_theta,
805+
SparseMatrixCOO<double>& solver_matrix)
794806
{
795-
const double theta = grid_.theta(i_theta);
796-
for (int i_r = grid_.numberSmootherCircles(); i_r < grid_.nr(); i_r++) {
797-
const int global_index = grid_.index(i_r, i_theta);
798-
const double r = grid_.radius(i_r);
807+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
808+
const LevelCache<DomainGeometry>& level_cache = DirectSolver<DomainGeometry>::level_cache_;
809+
810+
const double theta = grid.theta(i_theta);
811+
for (int i_r = grid.numberSmootherCircles(); i_r < grid.nr(); i_r++) {
812+
const int global_index = grid.index(i_r, i_theta);
813+
const double r = grid.radius(i_r);
799814

800815
double coeff_beta, arr, att, art, detDF;
801-
level_cache_.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
816+
level_cache.obtainValues(i_r, i_theta, global_index, r, theta, coeff_beta, arr, att, art, detDF);
802817

803818
// Build solver matrix at the current node
804-
nodeBuildSolverMatrixGive(i_r, i_theta, grid_, DirBC_Interior_, solver_matrix, arr, att, art, detDF,
805-
coeff_beta);
819+
nodeBuildSolverMatrixGive(i_r, i_theta, grid, DirectSolver<DomainGeometry>::DirBC_Interior_, solver_matrix, arr,
820+
att, art, detDF, coeff_beta);
806821
}
807822
}
808823

809824
// clang-format off
810825

811826
/* ------------------------------------------------------------------------ */
812827
/* If the indexing is not smoother-based, please adjust the access patterns */
813-
SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give::buildSolverMatrix()
828+
template <concepts::DomainGeometry DomainGeometry>
829+
SparseMatrixCOO<double> DirectSolver_COO_MUMPS_Give<DomainGeometry>::buildSolverMatrix()
814830
{
815-
const int n = grid_.numberOfNodes();
831+
const PolarGrid& grid = DirectSolver<DomainGeometry>::grid_;
832+
const int num_omp_threads = DirectSolver<DomainGeometry>::num_omp_threads_;
833+
834+
const int n = grid.numberOfNodes();
816835
const int nnz = getNonZeroCountSolverMatrix();
817836

818837
// Although the matrix is symmetric, we need to store all its entries, so we disable the symmetry.
819838
SparseMatrixCOO<double> solver_matrix(n, n, nnz);
820839
solver_matrix.is_symmetric(false);
821840

822-
#pragma omp parallel for num_threads(num_omp_threads_) if (nnz > 10'000)
841+
#pragma omp parallel for num_threads(num_omp_threads) if (nnz > 10'000)
823842
for (int i = 0; i < nnz; i++) {
824843
solver_matrix.value(i) = 0.0;
825844
}
826845

827-
if (num_omp_threads_ == 1) {
846+
if (num_omp_threads == 1) {
828847
/* Single-threaded execution */
829-
for (int i_r = 0; i_r < grid_.numberSmootherCircles(); i_r++) {
848+
for (int i_r = 0; i_r < grid.numberSmootherCircles(); i_r++) {
830849
buildSolverMatrixCircleSection(i_r, solver_matrix);
831850
}
832-
for (int i_theta = 0; i_theta < grid_.ntheta(); i_theta++) {
851+
for (int i_theta = 0; i_theta < grid.ntheta(); i_theta++) {
833852
buildSolverMatrixRadialSection(i_theta, solver_matrix);
834853
}
835854
}
836855
else {
837856
/* Multi-threaded execution: For Loops */
838-
const int num_circle_tasks = grid_.numberSmootherCircles();
839-
const int additional_radial_tasks = grid_.ntheta() % 3;
840-
const int num_radial_tasks = grid_.ntheta() - additional_radial_tasks;
857+
const int num_circle_tasks = grid.numberSmootherCircles();
858+
const int additional_radial_tasks = grid.ntheta() % 3;
859+
const int num_radial_tasks = grid.ntheta() - additional_radial_tasks;
841860

842-
#pragma omp parallel num_threads(num_omp_threads_)
861+
#pragma omp parallel num_threads(num_omp_threads)
843862
{
844863
#pragma omp for
845864
for (int circle_task = 0; circle_task < num_circle_tasks; circle_task += 3) {
846-
int i_r = grid_.numberSmootherCircles() - circle_task - 1;
865+
int i_r = grid.numberSmootherCircles() - circle_task - 1;
847866
buildSolverMatrixCircleSection(i_r, solver_matrix);
848867
}
849868
#pragma omp for
850869
for (int circle_task = 1; circle_task < num_circle_tasks; circle_task += 3) {
851-
int i_r = grid_.numberSmootherCircles() - circle_task - 1;
870+
int i_r = grid.numberSmootherCircles() - circle_task - 1;
852871
buildSolverMatrixCircleSection(i_r, solver_matrix);
853872
}
854873
#pragma omp for nowait
855874
for (int circle_task = 2; circle_task < num_circle_tasks; circle_task += 3) {
856-
int i_r = grid_.numberSmootherCircles() - circle_task - 1;
875+
int i_r = grid.numberSmootherCircles() - circle_task - 1;
857876
buildSolverMatrixCircleSection(i_r, solver_matrix);
858877
}
859878

‎include/DirectSolver/DirectSolver-COO-MUMPS-Give/directSolverGive.h‎

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,11 @@
44

55
#ifdef GMGPOLAR_USE_MUMPS
66

7-
class DirectSolver_COO_MUMPS_Give : public DirectSolver
7+
template <concepts::DomainGeometry DomainGeometry>
8+
class DirectSolver_COO_MUMPS_Give : public DirectSolver<DomainGeometry>
89
{
910
public:
10-
explicit DirectSolver_COO_MUMPS_Give(const PolarGrid& grid, const LevelCache& level_cache,
11+
explicit DirectSolver_COO_MUMPS_Give(const PolarGrid& grid, const LevelCache<DomainGeometry>& level_cache,
1112
const DomainGeometry& domain_geometry,
1213
const DensityProfileCoefficients& density_profile_coefficients,
1314
bool DirBC_Interior, int num_omp_threads);
@@ -89,4 +90,10 @@ class DirectSolver_COO_MUMPS_Give : public DirectSolver
8990
double detDF, double coeff_beta);
9091
};
9192

92-
#endif
93+
#include "applySymmetryShift.inl"
94+
#include "buildSolverMatrix.inl"
95+
#include "directSolverGive.inl"
96+
#include "initializeMumps.inl"
97+
#include "matrixStencil.inl"
98+
99+
#endif

0 commit comments

Comments
 (0)