Skip to content

Commit dfda389

Browse files
dradicedradice
authored
Merge branch 'main' into cargrid-output
2 parents 2a0853b + 49909af commit dfda389

16 files changed

+2227
-80
lines changed

‎inputs/z4c/awa/z4c_gauge_wave.athinput

+76
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
# AthenaXXX input file for Z4c linear wave tests
2+
3+
<comment>
4+
problem = z4c linear waves
5+
reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
6+
7+
<job>
8+
basename = gauge # problem ID: basename of output filenames
9+
10+
<mesh>
11+
nghost = 2 # Number of ghost cells
12+
nx1 = 50 # Number of zones in X1-direction
13+
x1min = 0.0 # minimum value of X1
14+
x1max = 1.0 # maximum value of X1
15+
ix1_bc = periodic # inner-X1 boundary flag
16+
ox1_bc = periodic # outer-X1 boundary flag
17+
18+
nx2 = 4 # Number of zones in X2-direction
19+
x2min = 0.0 # minimum value of X2
20+
x2max = 1.0 # maximum value of X2
21+
ix2_bc = periodic # inner-X2 boundary flag
22+
ox2_bc = periodic # outer-X2 boundary flag
23+
24+
nx3 = 4 # Number of zones in X3-direction
25+
x3min = 0.0 # minimum value of X3
26+
x3max = 1.0 # maximum value of X3
27+
ix3_bc = periodic # inner-X3 boundary flag
28+
ox3_bc = periodic # outer-X3 boundary flag
29+
30+
<meshblock>
31+
nx1 = 50 # Number of cells in each MeshBlock, X1-dir
32+
nx2 = 4 # Number of cells in each MeshBlock, X2-dir
33+
nx3 = 4 # Number of cells in each MeshBlock, X3-dir
34+
35+
<time>
36+
evolution = dynamic # dynamic/kinematic/static
37+
integrator = rk4 # time integration algorithm
38+
cfl_number = 0.5 # The Courant, Friedrichs, & Lewy (CFL) Number
39+
nlim = -1 # cycle limit (no limit if <0)
40+
tlim = 1000.0 # time limit
41+
ndiag = 1 # cycles between diagostic output
42+
43+
<z4c>
44+
lapse_harmonic = 0.0 # Harmonic lapse parameter mu_L
45+
lapse_oplog = 2.0 # 1+log lapse parameter
46+
shift_eta = 2.0 # Shift damping term
47+
diss = 0.02 # Kreiss-Oliger dissipation parameter
48+
chi_div_floor = 0.00001 # Floor on the conformal factor
49+
damp_kappa1 = 0.02 # Constraint damping factor 1
50+
damp_kappa2 = 0.0 # Constraint damping factor 2
51+
52+
<problem>
53+
pgen_name = z4c_gauge_wave # problem generator name
54+
amp = 0.01
55+
56+
<output1>
57+
file_type = tab # legacy VTK output
58+
variable = z4c # variables to be output
59+
dt = 1 # time increment between outputs
60+
slice_x2 = 0.5
61+
slice_x3 = 0.5
62+
ghost_zones = true # switch to output ghost cells
63+
64+
<output2>
65+
file_type = hst # history data dump
66+
user_hist = false
67+
data_format = %12.5e # Optional data format string
68+
dt = 1 # time increment between outputs
69+
70+
#<output3>
71+
#file_type = tab # legacy VTK output
72+
#variable = z4c # variables to be output
73+
#dt = 0.025 # time increment between outputs
74+
#slice_x2 = 0.5
75+
#slice_x3 = 0.5
76+
#ghost_zones = true # switch to output ghost cells

‎inputs/z4c/awa/z4c_stability.athinput

+76
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
# AthenaXXX input file for Z4c linear wave tests
2+
3+
<comment>
4+
problem = z4c linear waves
5+
reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
6+
7+
<job>
8+
basename = stability # problem ID: basename of output filenames
9+
10+
<mesh>
11+
nghost = 2 # Number of ghost cells
12+
nx1 = 4 # Number of zones in X1-direction
13+
x1min = 0.0 # minimum value of X1
14+
x1max = 1.0 # maximum value of X1
15+
ix1_bc = periodic # inner-X1 boundary flag
16+
ox1_bc = periodic # outer-X1 boundary flag
17+
18+
nx2 = 4 # Number of zones in X2-direction
19+
x2min = 0.0 # minimum value of X2
20+
x2max = 1.0 # maximum value of X2
21+
ix2_bc = periodic # inner-X2 boundary flag
22+
ox2_bc = periodic # outer-X2 boundary flag
23+
24+
nx3 = 50 # Number of zones in X3-direction
25+
x3min = 0.0 # minimum value of X3
26+
x3max = 1.0 # maximum value of X3
27+
ix3_bc = periodic # inner-X3 boundary flag
28+
ox3_bc = periodic # outer-X3 boundary flag
29+
30+
<meshblock>
31+
nx1 = 4 # Number of cells in each MeshBlock, X1-dir
32+
nx2 = 4 # Number of cells in each MeshBlock, X2-dir
33+
nx3 = 50 # Number of cells in each MeshBlock, X3-dir
34+
35+
<time>
36+
evolution = dynamic # dynamic/kinematic/static
37+
integrator = rk4 # time integration algorithm
38+
cfl_number = 0.5 # The Courant, Friedrichs, & Lewy (CFL) Number
39+
nlim = -1 # cycle limit (no limit if <0)
40+
tlim = 1000.0 # time limit
41+
ndiag = 1 # cycles between diagostic output
42+
43+
<z4c>
44+
diss = 0.02
45+
lapse_harmonic = 0.0 # Harmonic lapse parameter mu_L
46+
lapse_oplog = 1.0 # 1+log lapse parameter
47+
shift_eta = 2.0 # Shift damping term
48+
chi_div_floor = 0.00001
49+
damp_kappa1 = 0.00 # Constraint damping factor 1
50+
damp_kappa2 = 0.02
51+
52+
53+
<problem>
54+
pgen_name = z4c_stability # problem generator name
55+
user_hist = true # enroll user-defined history function
56+
rho = 1
57+
58+
<output1>
59+
file_type = tab # legacy VTK output
60+
variable = z4c # variables to be output
61+
dt = 100 # time increment between outputs
62+
slice_x1 = 0.5
63+
slice_x2 = 0.5
64+
ghost_zones = true # switch to output ghost cells
65+
66+
<output2>
67+
file_type = hst # history data dump
68+
user_hist = true
69+
data_format = %12.5e # Optional data format string
70+
dt = 1 # time increment between outputs
71+
72+
<output2>
73+
file_type = hst # history data dump
74+
user_hist = false
75+
data_format = %12.5e # Optional data format string
76+
dt = 1 # time increment between outputs

‎inputs/z4c/spectre_bbh/z4c_spectre_bbh.athinput

+3-3
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ basename = bbh # problem ID: basename of output filenames
66

77
<problem>
88
pgen_name = z4c_spectre_bbh
9-
id_filename_glob = BbhVolume*.h5
9+
id_filename_glob = /home/nilsvu/2024-09-bbh-athenak/ID/BbhVolume*.h5
1010
id_subfile_name = VolumeData
1111
id_observation_step = -1
1212

@@ -41,7 +41,7 @@ integrator = rk4 # time integration algorithm
4141
cfl_number = 0.1 # The Courant, Friedrichs, & Lewy (CFL) Number
4242
nlim = 50 # cycle limit
4343
tlim = 50 # time limit
44-
ndiag = 1 # cycles between diagostic output
44+
ndiag = 1 # cycles between diagnostic output
4545

4646
<mesh_refinement>
4747
refinement = static #adaptive # type of refinement
@@ -52,7 +52,7 @@ refinement_interval = 1
5252

5353
# amr chi method
5454
<z4c_amr>
55-
method = chi_min
55+
method = chi
5656
chi_min = 0.2
5757

5858
<refinement1>

‎src/driver/driver.cpp

+44-29
Original file line numberDiff line numberDiff line change
@@ -165,8 +165,8 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
165165
nimp_stages = 3;
166166
nexp_stages = 2;
167167
cfl_limit = 1.0;
168-
gam0[0] = 0.0;
169-
gam1[0] = 1.0;
168+
gam0[0] = 1.0;
169+
gam1[0] = 0.0;
170170
beta[0] = 1.0;
171171

172172
gam0[1] = 0.5;
@@ -185,6 +185,48 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
185185
a_twid[2][1] = 0.25;
186186
a_twid[2][2] = 0.25;
187187
a_impl = 0.5;
188+
} else if (integrator == "imex2+") {
189+
// IMEX(4,3,2): Krapp et al. (2024, arXiv:2310.04435), Eq.30.
190+
// three-stage explicit, four-stage implicit, second-order ImEx
191+
// two implicit stages added, adapting Athenak's overall architecture
192+
// Note explicit steps may not reduce to RK2 based on the parameters chosen
193+
nimp_stages = 4;
194+
nexp_stages = 3;
195+
cfl_limit = 1.0;
196+
gamma = 1.707106781186547; //1+1/sqrt(2)
197+
gam0[0] = 1.0;
198+
gam1[0] = 0.0;
199+
beta[0] = gamma;
200+
201+
gam0[1] = (2.0*gamma-1.0)/(2.0*gamma*gamma);
202+
gam1[1] = (1.0-(2.0*gamma-1.0)/(2.0*gamma*gamma));
203+
beta[1] = 1.0/(2.0*gamma);
204+
205+
gam0[2] = 1.0;
206+
gam1[2] = 0.0;
207+
beta[2] = 0.0;
208+
209+
a_twid[0][0] = 0.0;
210+
a_twid[0][1] = 0.0;
211+
a_twid[0][2] = 0.0;
212+
a_twid[0][3] = 0.0;
213+
214+
a_twid[1][0] = 0.0;
215+
a_twid[1][1] = 0.0;
216+
a_twid[1][2] = 0.0;
217+
a_twid[1][3] = 0.0;
218+
219+
a_twid[2][0] = 0.0;
220+
a_twid[2][1] = 0.0;
221+
a_twid[2][2] = (1.0-2.0*gamma*gamma)/2.0/gamma;
222+
a_twid[2][3] = 0.0;
223+
224+
a_twid[3][0] = 0.0;
225+
a_twid[3][1] = 0.0;
226+
a_twid[3][2] = 0.0;
227+
a_twid[3][3] = 0.0;
228+
229+
a_impl = gamma;
188230
} else if (integrator == "imex3") {
189231
// IMEX-SSP3(4,3,3): Pareschi & Russo (2005) Table VI.
190232
// three-stage explicit, four-stage implicit, third-order ImEx
@@ -227,33 +269,6 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
227269
a_twid[3][2] = (4.0*(b + e + a) - 1.0)/6.0;
228270
a_twid[3][3] = 2.0*(1.0 - a)/3.0;
229271
a_impl = a;
230-
} else if (integrator == "imex+") {
231-
// IMEX(2,3,2): Krapp et al. (2024, arXiv:2310.04435), Eq.30.
232-
// three-stage explicit, two-stage implicit, second-order ImEx
233-
// Note explicit steps may not reduce to RK2 based on the parameters chosen
234-
nimp_stages = 2;
235-
nexp_stages = 3;
236-
cfl_limit = 1.0;
237-
Real gamma = 1.7071067811865475; // 1 + 1/sqrt(2)
238-
gam0[0] = 0.0;
239-
gam1[0] = 1.0;
240-
beta[0] = gamma;
241-
242-
gam0[1] = (2.0*gamma-1.0)/(2.0*gamma*gamma);
243-
gam1[1] = 1.0-(2.0*gamma-1.0)/(2.0*gamma*gamma);
244-
beta[1] = 1.0/(2.0*gamma);
245-
246-
gam0[2] = 1.0;
247-
gam1[2] = 0.0;
248-
beta[2] = 0.0;
249-
250-
a_twid[0][0] = (1.0-2*gamma*gamma)/(2.0*gamma);
251-
a_twid[0][1] = 0.0;
252-
253-
a_twid[1][0] = 0.0;
254-
a_twid[1][1] = 0.0;
255-
256-
a_impl = gamma;
257272
// Error, unrecognized integrator name.
258273
} else {
259274
std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__

‎src/driver/driver.hpp

+1
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@ class Driver {
4444
Real delta[4]; // weights for updating the intermediate stage (u1)
4545
Real a_twid[4][4], a_impl; // matrix elements for implicit stages in ImEx
4646
Real cfl_limit; // maximum CFL number for integrator
47+
Real gamma; // gamma value for the IMEX_new integrator
4748
Kokkos::Timer* pwall_clock_; // timer for tracking the wall clock
4849
Real wall_time;
4950

‎src/ion-neutral/ion-neutral.hpp

+2
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,8 @@ struct IonNeutralTaskIDs {
5757
TaskID n_newdt;
5858
TaskID i_clear;
5959
TaskID n_clear;
60+
TaskID i_srctrms;
61+
TaskID n_srctrms;
6062
};
6163

6264
namespace ion_neutral {

‎src/ion-neutral/ion-neutral_tasks.cpp

+22-5
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,15 @@ void IonNeutral::AssembleIonNeutralTasks(
4747
id.i_sendf = tl["stagen"]->AddTask(&MHD::SendFlux, pmhd, id.i_flux);
4848
id.i_recvf = tl["stagen"]->AddTask(&MHD::RecvFlux, pmhd, id.i_sendf);
4949
id.i_rkupdt = tl["stagen"]->AddTask(&MHD::RKUpdate, pmhd, id.i_recvf);
50+
id.i_srctrms = tl["stagen"]->AddTask(&MHD::MHDSrcTerms, pmhd, id.i_rkupdt);
5051

51-
id.n_flux = tl["stagen"]->AddTask(&Hydro::Fluxes, phyd, id.i_rkupdt);
52+
id.n_flux = tl["stagen"]->AddTask(&Hydro::Fluxes, phyd, id.i_srctrms);
5253
id.n_sendf = tl["stagen"]->AddTask(&Hydro::SendFlux, phyd, id.n_flux);
5354
id.n_recvf = tl["stagen"]->AddTask(&Hydro::RecvFlux, phyd, id.n_sendf);
5455
id.n_rkupdt = tl["stagen"]->AddTask(&Hydro::RKUpdate, phyd, id.n_recvf);
56+
id.n_srctrms = tl["stagen"]->AddTask(&Hydro::HydroSrcTerms, phyd, id.n_rkupdt);
5557

56-
id.impl = tl["stagen"]->AddTask(&IonNeutral::ImpRKUpdate, this, id.n_rkupdt);
58+
id.impl = tl["stagen"]->AddTask(&IonNeutral::ImpRKUpdate, this, id.n_srctrms);
5759
id.i_restu = tl["stagen"]->AddTask(&MHD::RestrictU, pmhd, id.impl);
5860
id.n_restu = tl["stagen"]->AddTask(&Hydro::RestrictU, phyd, id.i_restu);
5961

@@ -185,9 +187,24 @@ TaskStatus IonNeutral::ImpRKUpdate(Driver *pdriver, int estage) {
185187
// equations for ion-neutral drag.
186188
// Only required for istage = (1,2,3,[4])
187189
if (estage < pdriver->nexp_stages) {
188-
Real gamma_adt = drag_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
189-
Real xi_adt = ionization_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
190-
Real alpha_adt = recombination_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
190+
Real gamma_adt;
191+
Real xi_adt;
192+
Real alpha_adt;
193+
194+
// Condition to set gamma_adt, xi_adt, and alpha_adt to zero
195+
if (istage < 3 && pdriver->integrator == "imex2+") {
196+
gamma_adt = 0.0;
197+
xi_adt = 0.0;
198+
alpha_adt = 0.0;
199+
} else {
200+
gamma_adt = drag_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
201+
xi_adt = ionization_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
202+
alpha_adt = recombination_coeff * (pdriver->a_impl) * (pmy_pack->pmesh->dt);
203+
}
204+
205+
//Real gamma_adt = drag_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
206+
//Real xi_adt = ionization_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt);
207+
//Real alpha_adt = recombination_coeff*(pdriver->a_impl)*(pmy_pack->pmesh->dt)
191208
auto ui = pmhd->u0;
192209
auto un = phyd->u0;
193210
par_for("imex_imp",DevExeSpace(),0,nmb1,0,(n3-1),0,(n2-1),0,(n1-1),

‎src/outputs/binary.cpp

+3-2
Original file line numberDiff line numberDiff line change
@@ -247,15 +247,16 @@ void MeshBinaryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) {
247247
std::size_t myoffset=header_offset+data_size*ns_mbs+data_size*m;
248248
// every rank has a MB to write, so write collectively
249249
if (m < noutmbs_min) {
250-
if (binfile.Write_any_type_at_all(pdata,(data_size),myoffset,"byte") != 1) {
250+
if (binfile.Write_any_type_at_all(pdata,(data_size),myoffset,"byte")
251+
!= data_size) {
251252
std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__
252253
<< std::endl << "binary data not written correctly to binary file, "
253254
<< "binary file is broken." << std::endl;
254255
exit(EXIT_FAILURE);
255256
}
256257
// some ranks are finished writing, so use non-collective write
257258
} else if (m < pm->nmb_thisrank) {
258-
if (binfile.Write_any_type_at(pdata,(data_size),myoffset,"byte") != 1) {
259+
if (binfile.Write_any_type_at(pdata,(data_size),myoffset,"byte") != data_size) {
259260
std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__
260261
<< std::endl << "binary data not written correctly to binary file, "
261262
<< "binary file is broken." << std::endl;

‎src/outputs/coarsened_binary.cpp

+4-2
Original file line numberDiff line numberDiff line change
@@ -507,15 +507,17 @@ void CoarsenedBinaryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) {
507507
std::size_t myoffset=header_offset+data_size*ns_mbs+data_size*m;
508508
// every rank has a MB to write, so write collectively
509509
if (m < noutmbs_min) {
510-
if (cbinfile.Write_any_type_at_all(pdata,(data_size),myoffset,"byte") != 1) {
510+
if (cbinfile.Write_any_type_at_all(pdata,(data_size),myoffset,"byte")
511+
!= data_size) {
511512
std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__
512513
<< std::endl << "binary data not written correctly to binary file, "
513514
<< "binary file is broken." << std::endl;
514515
exit(EXIT_FAILURE);
515516
}
516517
// some ranks are finished writing, so use non-collective write
517518
} else if (m < pm->nmb_thisrank) {
518-
if (cbinfile.Write_any_type_at(pdata,(data_size),myoffset,"byte") != 1) {
519+
if (cbinfile.Write_any_type_at(pdata,(data_size),myoffset,"byte")
520+
!= data_size) {
519521
std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__
520522
<< std::endl << "binary data not written correctly to binary file, "
521523
<< "binary file is broken." << std::endl;

0 commit comments

Comments
 (0)