Skip to content

Commit d8f9ab0

Browse files
zingalezingale
authored
Merge branch 'development' into vode_const
2 parents c937733 + df3029b commit d8f9ab0

File tree

12 files changed

+215
-99
lines changed

12 files changed

+215
-99
lines changed

‎Docs/source/unit_tests.rst‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -117,7 +117,7 @@ by options in the input file.
117117
One-zone tests
118118
==============
119119

120-
.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, nse_net_cell, part_func_cell
120+
.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, nse_net_cell, part_func_cell, interp_cell
121121

122122
* ``burn_cell`` :
123123

@@ -138,6 +138,11 @@ One-zone tests
138138
given a $\rho$, $T$, and $X_k$, call the equation of state and print out
139139
the thermodynamic information. See :ref:`sec:eos_cell` for more information.
140140

141+
* ``interp_cell`` :
142+
143+
This tests the cubic interpolant used in pynucastro networks for interpolating
144+
a rate that is given as pairs of $(T, N_A \langle \sigma v \rangle)$.
145+
141146
* ``jac_cell`` :
142147

143148
for a single thermodynamic state, compute the analytic Jacobian

‎integration/VODE/vode_dvjust.H‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ void dvjust (int IORD, BurnT& state, DvodeT& vstate)
2424
const int NQM1 = vstate.NQ - 1;
2525
const int NQM2 = vstate.NQ - 2;
2626

27-
amrex::Real HSUM{}, XI{}, ALPH0{}, ALPH1{}, PROD{}, XIOLD{}, T1{};
27+
amrex::Real XI{};
2828

2929
// Check to see if the order is being increased or decreased.
3030

@@ -35,7 +35,7 @@ void dvjust (int IORD, BurnT& state, DvodeT& vstate)
3535
vstate.el(j) = 0.0_rt;
3636
}
3737
vstate.el(3) = 1.0_rt;
38-
HSUM = 0.0_rt;
38+
amrex::Real HSUM = 0.0_rt;
3939
for (int j = 1; j <= NQM2; ++j) {
4040
// Construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)).
4141
HSUM += vstate.tau(j);
@@ -62,11 +62,11 @@ void dvjust (int IORD, BurnT& state, DvodeT& vstate)
6262
}
6363

6464
vstate.el(3) = 1.0_rt;
65-
ALPH0 = -1.0_rt;
66-
ALPH1 = 1.0_rt;
67-
PROD = 1.0_rt;
68-
XIOLD = 1.0_rt;
69-
HSUM = vstate.HSCAL;
65+
amrex::Real ALPH0 = -1.0_rt;
66+
amrex::Real ALPH1 = 1.0_rt;
67+
amrex::Real PROD = 1.0_rt;
68+
amrex::Real XIOLD = 1.0_rt;
69+
amrex::Real HSUM = vstate.HSCAL;
7070

7171
if (vstate.NQ != 1) {
7272
for (int j = 1; j <= NQM1; ++j) {
@@ -84,7 +84,7 @@ void dvjust (int IORD, BurnT& state, DvodeT& vstate)
8484
}
8585
}
8686

87-
T1 = (-ALPH0 - ALPH1) / PROD;
87+
amrex::Real T1 = (-ALPH0 - ALPH1) / PROD;
8888
// Load column L + 1 in YH array.
8989
for (int i = 1; i <= int_neqs; ++i) {
9090
vstate.yh(i,vstate.L+1) = T1 * vstate.yh(i,VODE_LMAX);

‎integration/VODE/vode_dvstep.H‎

Lines changed: 14 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -82,14 +82,13 @@ int dvstep (BurnT& state, DvodeT& vstate)
8282
constexpr amrex::Real ONEPSM = 1.00001e0_rt;
8383
constexpr amrex::Real THRESH = 1.5e0_rt;
8484

85-
amrex::Real CNQUOT{}, DDN{}, DSM{}, DUP{}, TOLD{};
86-
amrex::Real FLOTL{}, R{};
85+
amrex::Real DSM{};
8786
int NCF{}, NFLAG{};
8887

8988
constexpr int int_neqs = integrator_neqs<BurnT>();
9089

9190
int kflag = 0;
92-
TOLD = vstate.tn;
91+
amrex::Real TOLD = vstate.tn;
9392
NCF = 0;
9493
vstate.JCUR = 0;
9594
NFLAG = 0;
@@ -117,7 +116,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
117116

118117
// Rescale the history array for a change in H by a factor of ETA.
119118

120-
R = 1.0_rt;
119+
amrex::Real R = 1.0_rt;
121120

122121
for (int j = 2; j <= vstate.L; ++j) {
123122
R *= vstate.ETA;
@@ -201,7 +200,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
201200
vstate.ETA = amrex::max(vstate.ETA, HMIN / std::abs(vstate.H));
202201

203202
// Rescale the history array for a change in H by a factor of ETA.
204-
R = 1.0_rt;
203+
amrex::Real R = 1.0_rt;
205204

206205
for (int j = 2; j <= vstate.L; ++j) {
207206
R *= vstate.ETA;
@@ -357,7 +356,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
357356
vstate.ETAMAX = ETAMX2;
358357
}
359358

360-
R = 1.0_rt / vstate.tq(2);
359+
amrex::Real R = 1.0_rt / vstate.tq(2);
361360
for (int i = 1; i <= int_neqs; ++i) {
362361
vstate.acor(i) *= R;
363362
}
@@ -389,15 +388,15 @@ int dvstep (BurnT& state, DvodeT& vstate)
389388
if (kflag > KFC) {
390389

391390
// Compute ratio of new H to current H at the current order.
392-
FLOTL = vstate.L;
391+
amrex::Real FLOTL = vstate.L;
393392
vstate.ETA = 1.0_rt / (std::pow(BIAS2 * DSM, 1.0_rt / FLOTL) + ADDON);
394393
vstate.ETA = amrex::max(vstate.ETA, HMIN / std::abs(vstate.H), ETAMIN);
395394
if ((kflag <= -2) && (vstate.ETA > ETAMXF)) {
396395
vstate.ETA = ETAMXF;
397396
}
398397

399398
// Rescale the history array for a change in H by a factor of ETA.
400-
R = 1.0_rt;
399+
amrex::Real R = 1.0_rt;
401400

402401
for (int j = 2; j <= vstate.L; ++j) {
403402
R *= vstate.ETA;
@@ -434,7 +433,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
434433
vstate.NQWAIT = vstate.L;
435434

436435
// Rescale the history array for a change in H by a factor of ETA.
437-
R = 1.0_rt;
436+
amrex::Real R = 1.0_rt;
438437

439438
for (int j = 2; j <= vstate.L; ++j) {
440439
R *= vstate.ETA;
@@ -477,7 +476,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
477476
bool already_set_eta = false;
478477

479478
// Compute ratio of new H to current H at the current order.
480-
FLOTL = vstate.L;
479+
amrex::Real FLOTL = vstate.L;
481480
const amrex::Real ETAQ = 1.0_rt / (std::pow(BIAS2 * DSM, 1.0_rt / FLOTL) + ADDON);
482481
if (vstate.NQWAIT == 0) {
483482
vstate.NQWAIT = 2;
@@ -486,7 +485,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
486485

487486
if (vstate.NQ != 1) {
488487
// Compute ratio of new H to current H at the current order less one.
489-
DDN = 0.0_rt;
488+
amrex::Real DDN = 0.0_rt;
490489
for (int i = 1; i <= int_neqs; ++i) {
491490
DDN += (vstate.yh(i,vstate.L) * vstate.ewt(i)) * (vstate.yh(i,vstate.L) * vstate.ewt(i));
492491
}
@@ -498,11 +497,11 @@ int dvstep (BurnT& state, DvodeT& vstate)
498497

499498
if (vstate.L != VODE_LMAX) {
500499
// Compute ratio of new H to current H at current order plus one.
501-
CNQUOT = (vstate.tq(5) / vstate.CONP) * std::pow(vstate.H / vstate.tau(2), vstate.L);
500+
amrex::Real CNQUOT = (vstate.tq(5) / vstate.CONP) * std::pow(vstate.H / vstate.tau(2), vstate.L);
502501
for (int i = 1; i <= int_neqs; ++i) {
503502
vstate.savf(i) = vstate.acor(i) - CNQUOT * vstate.yh(i,VODE_LMAX);
504503
}
505-
DUP = 0.0_rt;
504+
amrex::Real DUP = 0.0_rt;
506505
for (int i = 1; i <= int_neqs; ++i) {
507506
DUP += (vstate.savf(i) * vstate.ewt(i)) * (vstate.savf(i) * vstate.ewt(i));
508507
}
@@ -547,7 +546,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
547546
if (vstate.n_step <= 10) {
548547
vstate.ETAMAX = ETAMX2;
549548
}
550-
R = 1.0_rt / vstate.tq(2);
549+
amrex::Real R = 1.0_rt / vstate.tq(2);
551550
for (int i = 1; i <= int_neqs; ++i) {
552551
vstate.acor(i) *= R;
553552
}
@@ -561,7 +560,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
561560
if (vstate.n_step <= 10) {
562561
vstate.ETAMAX = ETAMX2;
563562
}
564-
R = 1.0_rt / vstate.tq(2);
563+
amrex::Real R = 1.0_rt / vstate.tq(2);
565564
for (int i = 1; i <= int_neqs; ++i) {
566565
vstate.acor(i) *= R;
567566
}

‎screening/screen.H‎

Lines changed: 25 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -185,7 +185,7 @@ number_t debye_huckel (const plasma_state_t<number_t>& state,
185185
template <typename number_t>
186186
AMREX_GPU_HOST_DEVICE AMREX_INLINE
187187
number_t actual_screen5 (const plasma_state_t<number_t>& state,
188-
const scrn::screen_factors_t& scn_fac)
188+
const scrn::screen_factors_t& scn_fac)
189189
{
190190
// this subroutine calculates screening factors and their derivatives
191191
// for nuclear reaction rates in the weak, intermediate and strong regimes.
@@ -200,17 +200,18 @@ number_t actual_screen5 (const plasma_state_t<number_t>& state,
200200

201201

202202
// fact = 2^(1/3)
203-
const amrex::Real fact = 1.25992104989487e0_rt;
204-
const amrex::Real gamefx = 0.3e0_rt; // lower gamma limit for intermediate screening
205-
const amrex::Real gamefs = 0.8e0_rt; // upper gamma limit for intermediate screening
206-
const amrex::Real h12_max = 300.e0_rt;
203+
constexpr amrex::Real fact = 1.25992104989487e0_rt;
204+
constexpr amrex::Real gamefx = 0.3e0_rt; // lower gamma limit for intermediate screening
205+
constexpr amrex::Real gamefs = 0.8e0_rt; // upper gamma limit for intermediate screening
206+
constexpr amrex::Real inv_dg = 1.0_rt / (gamefs - gamefx);
207+
constexpr amrex::Real h12_max = 300.e0_rt;
207208

208209
// Get the ion data based on the input index
209-
amrex::Real z1 = scn_fac.z1;
210-
amrex::Real z2 = scn_fac.z2;
210+
const amrex::Real z1 = scn_fac.z1;
211+
const amrex::Real z2 = scn_fac.z2;
211212

212213
// calculate individual screening factors
213-
amrex::Real bb = z1 * z2;
214+
const amrex::Real bb = z1 * z2;
214215
number_t gamp = state.aa;
215216

216217
// In Eq.4 in Itoh:1979, this term is 2*Z_1*Z_2/(Z_1^(1/3) + Z_2^(1/3))
@@ -252,77 +253,64 @@ number_t actual_screen5 (const plasma_state_t<number_t>& state,
252253
// Full version of Eq. 19 in Graboske:1973 by considering weak regime
253254
// and Wallace:1982 Eq. A14. Here the degeneracy factor is assumed to be 1.
254255

255-
number_t h12w = bb * state.qlam0z;
256+
const number_t h12w = bb * state.qlam0z;
256257

257258
number_t h12 = h12w;
258259

259260
// intermediate and strong sceening regime
260261

261262
if (gamef > gamefx) {
262263

263-
// gamma_ij^(1/4)
264-
265-
number_t gamp14 = admath::pow(gamp, 0.25_rt);
264+
// gamma_ij^(1/4)
265+
const number_t gamp14 = admath::sqrt(admath::sqrt(gamp));
266266

267267
// Here we follow Eq. A9 in Wallace:1982
268268
// See Eq. 25 Alastuey:1978, Eq. 16 and 17 in Jancovici:1977 for reference
269-
number_t cc = 0.896434e0_rt * gamp * scn_fac.zhat
269+
const number_t cc = 0.896434e0_rt * gamp * scn_fac.zhat
270270
- 3.44740e0_rt * gamp14 * scn_fac.zhat2
271271
- 0.5551e0_rt * (admath::log(gamp) + scn_fac.lzav)
272272
- 2.996e0_rt;
273273

274274
// (3gamma_ij/tau_ij)^3
275-
number_t a3 = alph12 * alph12 * alph12;
275+
const number_t a3 = alph12 * alph12 * alph12;
276276

277277
// Part of Eq. 28 in Alastuey:1978
278-
number_t rr = (5.0_rt/32.0_rt) - alph12*(0.014e0_rt + 0.0128e0_rt*alph12);
278+
const number_t rr = (5.0_rt/32.0_rt) - alph12*(0.014e0_rt + 0.0128e0_rt*alph12);
279279

280280
// Part of Eq. 28 in Alastuey:1978
281-
number_t ss = tau12*rr;
281+
const number_t ss = tau12*rr;
282282

283283
// Part of Eq. 31 in Alastuey:1978
284-
number_t tt = -0.0098e0_rt + 0.0048e0_rt*alph12;
284+
const number_t tt = -0.0098e0_rt + 0.0048e0_rt*alph12;
285285

286286
// Part of Eq. 31 in Alastuey:1978
287-
number_t uu = 0.0055e0_rt + alph12*tt;
287+
const number_t uu = 0.0055e0_rt + alph12*tt;
288288

289289
// Part of Eq. 31 in Alastuey:1978
290-
number_t vv = gamef * alph12 * uu;
290+
const number_t vv = gamef * alph12 * uu;
291291

292292
// Exponent of Eq. 32 in Alastuey:1978, which uses Eq.28 and Eq.31
293293
// Strong screening factor
294294
h12 = cc - a3 * (ss + vv);
295295

296296
// See conclusion and Eq. 34 in Alastuey:1978
297297
// This is an extra factor to account for quantum effects
298-
rr = 1.0_rt - 0.0562e0_rt*a3;
299-
300-
number_t xlgfac;
298+
// In extreme case, limit to 0.77, see conclusion in Alastuey:1978
301299

302-
// In extreme case, rr is 0.77, see conclusion in Alastuey:1978
303-
if (rr >= 0.77e0_rt) {
304-
xlgfac = rr;
305-
} else {
306-
xlgfac = 0.77e0_rt;
307-
}
300+
const number_t xlgfac = admath::max(1.0_rt - 0.0562e0_rt*a3, 0.77_rt);
308301

309302
// Include the extra factor that accounts for quantum effects
310-
h12 = admath::log(xlgfac) + h12;
303+
h12 += admath::log(xlgfac);
311304

312305
// If gamma_ij < upper limit of intermediate regime
313306
// then it is in the intermediate regime, else strong screening.
314307
if (gamef <= gamefs) {
315-
amrex::Real dgamma = 1.0e0_rt/(gamefs - gamefx);
316-
317-
rr = dgamma*(gamefs - gamef);
318-
319-
ss = dgamma*(gamef - gamefx);
320-
321-
vv = h12;
308+
const number_t weight1 = inv_dg * (gamefs - gamef);
309+
const number_t weight2 = inv_dg * (gamef - gamefx);
322310

323311
// Then the screening factor is a combination
324312
// of the strong and weak screening factor.
325-
h12 = h12w*rr + vv*ss;
313+
h12 = h12w * weight1 + h12 * weight2;
326314
}
327315

328316
// end of intermediate and strong screening

‎unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
Initializing AMReX (24.10-178-g5a2390e00837)...
2-
AMReX (24.10-178-g5a2390e00837) initialized
1+
Initializing AMReX (25.12-25-g10b67a795031)...
2+
AMReX (25.12-25-g10b67a795031) initialized
33
starting the single zone burn...
44
Maximum Time (s): 0.01585
55
State Density (g/cm^3): 1000000000
@@ -13,21 +13,21 @@ RHS at t = 0
1313
ash 0.01230280576
1414
------------------------------------
1515
successful? 1
16-
- Hnuc = 5.277066077e+17
17-
- added e = 8.364149732e+15
18-
- final T = 1433686831
16+
- Hnuc = 5.277014299e+17
17+
- added e = 8.364067664e+15
18+
- final T = 1433682844
1919
------------------------------------
2020
e initial = 1.253426044e+18
21-
e final = 1.261790194e+18
21+
e final = 1.261790112e+18
2222
------------------------------------
2323
new mass fractions:
24-
C12 0.9657916864
24+
C12 0.965792022
2525
O16 1e-30
26-
ash 0.03420831361
26+
ash 0.03420797796
2727
------------------------------------
2828
species creation rates:
29-
omegadot(C12): -2.158253225
30-
omegadot(O16): 1.989224949e-43
31-
omegadot(ash): 2.158253225
29+
omegadot(C12): -2.158232048
30+
omegadot(O16): 7.735874803e-44
31+
omegadot(ash): 2.158232048
3232
number of steps taken: 381
33-
AMReX (24.10-178-g5a2390e00837) finalized
33+
AMReX (25.12-25-g10b67a795031) finalized

0 commit comments

Comments
 (0)