Skip to content

Commit f40e59f

Browse files
haykhhaykh
authored
Merge pull request entity-toolkit#123 from entity-toolkit/dev/vector-potential
MPI support for A_phi calculation in GR
2 parents 52cec30 + ec47ff8 commit f40e59f

6 files changed

Lines changed: 267 additions & 106 deletions

File tree

‎dev/nix/kokkos.nix‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ let
1010
pversion = "4.6.01";
1111
compilerPkgs = {
1212
"HIP" = with pkgs.rocmPackages; [
13+
llvm.rocm-merged-llvm
1314
rocm-core
1415
clr
1516
rocthrust
@@ -18,11 +19,13 @@ let
1819
rocm-smi
1920
];
2021
"CUDA" = with pkgs.cudaPackages; [
22+
llvmPackages_18.clang-tools
2123
cudatoolkit
2224
cuda_cudart
2325
pkgs.gcc13
2426
];
2527
"NONE" = [
28+
pkgs.llvmPackages_18.clang-tools
2629
pkgs.gcc13
2730
];
2831
};
@@ -36,9 +39,7 @@ let
3639
"HIP" = [
3740
"-D Kokkos_ENABLE_HIP=ON"
3841
"-D Kokkos_ARCH_${getArch { }}=ON"
39-
# remove leading AMD_
4042
"-D AMDGPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}"
41-
"-D CMAKE_C_COMPILER=hipcc"
4243
"-D CMAKE_CXX_COMPILER=hipcc"
4344
];
4445
"CUDA" = [

‎dev/nix/shell.nix‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,6 @@ pkgs.mkShell {
4242
zlib
4343
cmake
4444

45-
llvmPackages_18.clang-tools
46-
libgcc
47-
4845
adios2Pkg
4946
kokkosPkg
5047

‎pgens/wald/pgen.hpp‎

Lines changed: 133 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#include "arch/kokkos_aliases.h"
88
#include "arch/traits.h"
99
#include "utils/comparators.h"
10+
#include "utils/error.h"
1011
#include "utils/formatting.h"
1112
#include "utils/log.h"
1213
#include "utils/numeric.h"
@@ -17,14 +18,31 @@
1718
#include "framework/domain/domain.h"
1819
#include "framework/domain/metadomain.h"
1920

21+
#include <string>
2022
#include <vector>
2123

24+
enum InitFieldGeometry {
25+
Wald,
26+
Vertical,
27+
};
28+
2229
namespace user {
2330
using namespace ntt;
2431

2532
template <class M, Dimension D>
2633
struct InitFields {
27-
InitFields(M metric_) : metric { metric_ } {}
34+
InitFields(M metric_, const std::string& init_field_geometry)
35+
: metric { metric_ } {
36+
if (init_field_geometry == "wald") {
37+
field_geometry = InitFieldGeometry::Wald;
38+
} else if (init_field_geometry == "vertical") {
39+
field_geometry = InitFieldGeometry::Vertical;
40+
} else {
41+
raise::Error(fmt::format("Unrecognized field geometry: %s",
42+
init_field_geometry.c_str()),
43+
HERE);
44+
}
45+
}
2846

2947
Inline auto A_3(const coord_t<D>& x_Cd) const -> real_t {
3048
return HALF * (metric.template h_<3, 3>(x_Cd) +
@@ -83,102 +101,132 @@ namespace user {
83101

84102
Inline auto bx3(
85103
const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j + HALF )
86-
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
87-
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
88-
89-
x0m[0] = xi[0];
90-
x0m[1] = xi[1] - HALF;
91-
x0p[0] = xi[0];
92-
x0p[1] = xi[1] + HALF;
93-
94-
real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
95-
return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP;
104+
if (field_geometry == InitFieldGeometry::Wald) {
105+
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
106+
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
107+
108+
x0m[0] = xi[0];
109+
x0m[1] = xi[1] - HALF;
110+
x0p[0] = xi[0];
111+
x0p[1] = xi[1] + HALF;
112+
113+
real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
114+
return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP;
115+
} else if (field_geometry == InitFieldGeometry::Vertical) {
116+
return ZERO;
117+
} else {
118+
raise::KernelError(HERE, "Unrecognized field geometry");
119+
return ZERO;
120+
}
96121
}
97122

98123
Inline auto dx1(const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j )
99-
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
100-
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
101-
102-
real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) };
103-
real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
104-
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
105-
real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) };
106-
real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) };
107-
108-
// D1 at ( i + HALF , j )
109-
x0m[0] = xi[0] - HALF;
110-
x0m[1] = xi[1];
111-
x0p[0] = xi[0] + HALF;
112-
x0p[1] = xi[1];
113-
real_t E1d { (A_0(x0p) - A_0(x0m)) };
114-
real_t D1d { E1d / alpha_iPj };
115-
116-
// D3 at ( i , j )
117-
x0m[0] = xi[0] - HALF - HALF;
118-
x0m[1] = xi[1];
119-
x0p[0] = xi[0] - HALF + HALF;
120-
x0p[1] = xi[1];
121-
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
122-
123-
real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d +
124-
metric.template h<1, 3>({ xi[0], xi[1] }) * D3d };
125-
126-
return D1u;
124+
if (field_geometry == InitFieldGeometry::Wald) {
125+
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
126+
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
127+
128+
real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) };
129+
real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
130+
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
131+
real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) };
132+
real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) };
133+
134+
// D1 at ( i + HALF , j )
135+
x0m[0] = xi[0] - HALF;
136+
x0m[1] = xi[1];
137+
x0p[0] = xi[0] + HALF;
138+
x0p[1] = xi[1];
139+
real_t E1d { (A_0(x0p) - A_0(x0m)) };
140+
real_t D1d { E1d / alpha_iPj };
141+
142+
// D3 at ( i , j )
143+
x0m[0] = xi[0] - HALF - HALF;
144+
x0m[1] = xi[1];
145+
x0p[0] = xi[0] - HALF + HALF;
146+
x0p[1] = xi[1];
147+
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
148+
149+
real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d +
150+
metric.template h<1, 3>({ xi[0], xi[1] }) * D3d };
151+
152+
return D1u;
153+
} else if (field_geometry == InitFieldGeometry::Vertical) {
154+
return ZERO;
155+
} else {
156+
raise::KernelError(HERE, "Unrecognized field geometry");
157+
return ZERO;
158+
}
127159
}
128160

129161
Inline auto dx2(const coord_t<D>& x_Ph) const -> real_t { // at ( i , j + HALF )
130-
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
131-
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
132-
x0m[0] = xi[0];
133-
x0m[1] = xi[1] - HALF;
134-
x0p[0] = xi[0];
135-
x0p[1] = xi[1] + HALF;
136-
real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
137-
real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) };
138-
real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) };
139-
real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) };
140-
141-
real_t E2d { (A_0(x0p) - A_0(x0m)) };
142-
real_t D2d { E2d / alpha_ijP - (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP };
143-
real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d };
144-
145-
return D2u;
162+
if (field_geometry == InitFieldGeometry::Wald) {
163+
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
164+
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
165+
x0m[0] = xi[0];
166+
x0m[1] = xi[1] - HALF;
167+
x0p[0] = xi[0];
168+
x0p[1] = xi[1] + HALF;
169+
real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
170+
real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) };
171+
real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) };
172+
real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) };
173+
174+
real_t E2d { (A_0(x0p) - A_0(x0m)) };
175+
real_t D2d { E2d / alpha_ijP -
176+
(A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP };
177+
real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d };
178+
179+
return D2u;
180+
} else if (field_geometry == InitFieldGeometry::Vertical) {
181+
return ZERO;
182+
} else {
183+
raise::KernelError(HERE, "Unrecognized field geometry");
184+
return ZERO;
185+
}
146186
}
147187

148188
Inline auto dx3(const coord_t<D>& x_Ph) const -> real_t { // at ( i , j )
149-
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
150-
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
151-
real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
152-
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0], xi[1] }) };
153-
real_t beta_ij { metric.beta1({ xi[0], xi[1] }) };
154-
real_t alpha_ij { metric.alpha({ xi[0], xi[1] }) };
155-
real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) };
156-
157-
// D3 at ( i , j )
158-
x0m[0] = xi[0] - HALF;
159-
x0m[1] = xi[1];
160-
x0p[0] = xi[0] + HALF;
161-
x0p[1] = xi[1];
162-
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
163-
164-
// D1 at ( i + HALF , j )
165-
x0m[0] = xi[0] + HALF - HALF;
166-
x0m[1] = xi[1];
167-
x0p[0] = xi[0] + HALF + HALF;
168-
x0p[1] = xi[1];
169-
real_t E1d { (A_0(x0p) - A_0(x0m)) };
170-
real_t D1d { E1d / alpha_iPj };
171-
172-
if (cmp::AlmostZero(x_Ph[1])) {
173-
return metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
189+
if (field_geometry == InitFieldGeometry::Wald) {
190+
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
191+
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
192+
real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
193+
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0], xi[1] }) };
194+
real_t beta_ij { metric.beta1({ xi[0], xi[1] }) };
195+
real_t alpha_ij { metric.alpha({ xi[0], xi[1] }) };
196+
real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) };
197+
198+
// D3 at ( i , j )
199+
x0m[0] = xi[0] - HALF;
200+
x0m[1] = xi[1];
201+
x0p[0] = xi[0] + HALF;
202+
x0p[1] = xi[1];
203+
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
204+
205+
// D1 at ( i + HALF , j )
206+
x0m[0] = xi[0] + HALF - HALF;
207+
x0m[1] = xi[1];
208+
x0p[0] = xi[0] + HALF + HALF;
209+
x0p[1] = xi[1];
210+
real_t E1d { (A_0(x0p) - A_0(x0m)) };
211+
real_t D1d { E1d / alpha_iPj };
212+
213+
if (cmp::AlmostZero(x_Ph[1])) {
214+
return metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
215+
} else {
216+
return metric.template h<3, 3>({ xi[0], xi[1] }) * D3d +
217+
metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
218+
}
219+
} else if (field_geometry == InitFieldGeometry::Vertical) {
220+
return ZERO;
174221
} else {
175-
return metric.template h<3, 3>({ xi[0], xi[1] }) * D3d +
176-
metric.template h<1, 3>({ xi[0], xi[1] }) * D1d;
222+
raise::KernelError(HERE, "Unrecognized field geometry");
223+
return ZERO;
177224
}
178225
}
179226

180227
private:
181-
const M metric;
228+
const M metric;
229+
InitFieldGeometry field_geometry;
182230
};
183231

184232
template <SimEngine::type S, class M>
@@ -201,7 +249,8 @@ namespace user {
201249
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
202250
: arch::ProblemGenerator<S, M> { p }
203251
, global_domain { m }
204-
, init_flds { m.mesh().metric } {}
252+
, init_flds { m.mesh().metric,
253+
p.template get<std::string>("setup.init_field", "wald") } {}
205254
};
206255

207256
} // namespace user

‎pgens/wald/wald.toml‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,11 +4,11 @@
44
runtime = 100.0
55

66
[grid]
7-
resolution = [128, 128]
7+
resolution = [512, 512]
88
extent = [[1.0, 10.0]]
99

1010
[grid.metric]
11-
metric = "kerr_schild"
11+
metric = "qkerr_schild"
1212
qsph_r0 = 0.0
1313
qsph_h = 0.0
1414
ks_a = 0.95
@@ -38,10 +38,10 @@
3838
ppc0 = 2.0
3939

4040
[setup]
41+
init_field = "wald" # or "vertical"
4142

4243
[output]
43-
format = "hdf5"
44-
separate_files = false
44+
format = "hdf5"
4545

4646
[output.fields]
4747
interval_time = 1.0

‎src/framework/domain/metadomain.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -261,6 +261,9 @@ namespace ntt {
261261
#if defined(OUTPUT_ENABLED)
262262
out::Writer g_writer;
263263
checkpoint::Writer g_checkpoint_writer;
264+
#if defined(MPI_ENABLED)
265+
void CommunicateVectorPotential(unsigned short);
266+
#endif
264267
#endif
265268

266269
#if defined(MPI_ENABLED)

0 commit comments

Comments
 (0)