Skip to content

Commit

Permalink
use std::clamp to make some code clearer (#2905)
Browse files Browse the repository at this point in the history 
also swtich to std::min/max where appropriate (instead of amrex::)
  • Loading branch information
@zingale
zingale authored Jul 10, 2024
1 parent 0b1a3d5 commit 8b289c1
Show file tree
Hide file tree
Showing 6 changed files with 58 additions and 62 deletions.
4 changes: 2 additions & 2 deletions Source/hydro/Castro_ctu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,8 +360,8 @@ Castro::add_sdc_source_to_states(const Box& bx, const int idir, const Real dt,

if (n >= QFS && n <= QFS-1+NumSpec) {
// mass fractions should be in [0, 1]
qleft(i,j,k,n) = amrex::max(0.0_rt, amrex::min(1.0_rt, qleft(i,j,k,n)));
qright(i,j,k,n) = amrex::max(0.0_rt, amrex::min(1.0_rt, qright(i,j,k,n)));
qleft(i,j,k,n) = std::clamp(qleft(i,j,k,n), 0.0_rt, 1.0_rt);
qright(i,j,k,n) = std::clamp(qright(i,j,k,n), 0.0_rt, 1.0_rt);
}
}

Expand Down
41 changes: 20 additions & 21 deletions Source/hydro/HLL_solvers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,15 +231,15 @@ namespace HLL {


// signal speeds (E91, eq. 4.5)
Real bl = amrex::min(a1, ql[ivel] - cl);
Real br = amrex::max(a4, qr[ivel] + cr);
Real bl = std::min(a1, ql[ivel] - cl);
Real br = std::max(a4, qr[ivel] + cr);

Real bm = amrex::min(0.0_rt, bl);
Real bp = amrex::max(0.0_rt, br);
Real bm = std::min(0.0_rt, bl);
Real bp = std::max(0.0_rt, br);

Real bd = bp - bm;

if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) {
if (std::abs(bd) < small_hll * std::max(std::abs(bm), std::abs(bp))) {
return;
}

Expand Down Expand Up @@ -391,24 +391,23 @@ namespace HLL {
}


Real rl = amrex::max(ql(i,j,k,QRHO), small_dens);
Real rl = std::max(ql(i,j,k,QRHO), small_dens);

// pick left velocities based on direction
Real ul = ql(i,j,k,iu);

Real pl = amrex::max(ql(i,j,k,QPRES), small_pres);
Real pl = std::max(ql(i,j,k,QPRES), small_pres);

Real rr = amrex::max(qr(i,j,k,QRHO), small_dens);
Real rr = std::max(qr(i,j,k,QRHO), small_dens);

// pick right velocities based on direction
Real ur = qr(i,j,k,iu);

Real pr = amrex::max(qr(i,j,k,QPRES), small_pres);
Real pr = std::max(qr(i,j,k,QPRES), small_pres);

// now we essentially do the CGF solver to get p and u on the
// interface, but we won't use these in any flux construction.
Real csmall = amrex::max(small, amrex::max(small * qaux_arr(i,j,k,QC),
small * qaux_arr(i-sx,j-sy,k-sz,QC)));
Real csmall = amrex::max(small, small * qaux_arr(i,j,k,QC), small * qaux_arr(i-sx,j-sy,k-sz,QC));
Real cavg = 0.5_rt*(qaux_arr(i,j,k,QC) + qaux_arr(i-sx,j-sy,k-sz,QC));

Real gamcl = qaux_arr(i-sx,j-sy,k-sz,QGAMC);
Expand All @@ -422,14 +421,14 @@ namespace HLL {
#endif

Real wsmall = small_dens*csmall;
Real wl = amrex::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl)));
Real wr = amrex::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr)));
Real wl = std::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl)));
Real wr = std::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr)));

Real wwinv = 1.0_rt/(wl + wr);
Real pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))*wwinv;
Real ustar = ((wl*ul + wr*ur) + (pl - pr))*wwinv;

pstar = amrex::max(pstar, small_pres);
pstar = std::max(pstar, small_pres);
// for symmetry preservation, if ustar is really small, then we
// set it to zero
if (std::abs(ustar) < riemann_constants::smallu*0.5_rt*(std::abs(ul) + std::abs(ur))){
Expand Down Expand Up @@ -460,18 +459,18 @@ namespace HLL {
gamco = 0.5_rt*(gamcl + gamcr);
}

ro = amrex::max(small_dens, ro);
ro = std::max(small_dens, ro);

Real roinv = 1.0_rt/ro;
Real co = std::sqrt(std::abs(gamco*po*roinv));
co = amrex::max(csmall, co);
co = std::max(csmall, co);
Real co2inv = 1.0_rt/(co*co);

Real rstar = ro + (pstar - po)*co2inv;
rstar = amrex::max(small_dens, rstar);
rstar = std::max(small_dens, rstar);

Real cstar = std::sqrt(std::abs(gamco*pstar/rstar));
cstar = max(cstar, csmall);
cstar = std::max(cstar, csmall);

Real sgnm = std::copysign(1.0_rt, ustar);
Real spout = co - sgnm*uo;
Expand All @@ -489,7 +488,7 @@ namespace HLL {
}

Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt;
frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac));
frac = std::clamp(frac, 0.0_rt, 1.0_rt);

Real qint[NQ] = {0.0};

Expand All @@ -501,8 +500,8 @@ namespace HLL {
// now we do the HLLC construction

// use the simplest estimates of the wave speeds
Real S_l = amrex::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr));
Real S_r = amrex::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr));
Real S_l = std::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr));
Real S_r = std::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr));

// estimate of the contact speed -- this is Toro Eq. 10.8
Real S_c = (pr - pl + rl*ul*(S_l - ul) - rr*ur*(S_r - ur))/
Expand Down
42 changes: 21 additions & 21 deletions Source/hydro/flatten.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,16 +22,16 @@ Real flatten(int i, int j, int k,

int ishft = dp > 0.0_rt ? 1 : -1;

Real denom = amrex::max(small_pres, std::abs(q_arr(i+2,j,k,pres_comp) - q_arr(i-2,j,k,pres_comp)));
Real denom = std::max(small_pres, std::abs(q_arr(i+2,j,k,pres_comp) - q_arr(i-2,j,k,pres_comp)));
Real zeta = std::abs(dp) / denom;
Real z = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
Real z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

Real tst = 0.0_rt;
if (q_arr(i-1,j,k,QU) - q_arr(i+1,j,k,QU) >= 0.0_rt) {
tst = 1.0_rt;
}

Real tmp = amrex::min(q_arr(i+1,j,k,pres_comp), q_arr(i-1,j,k,pres_comp));
Real tmp = std::min(q_arr(i+1,j,k,pres_comp), q_arr(i-1,j,k,pres_comp));

Real chi = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
Expand All @@ -41,23 +41,23 @@ Real flatten(int i, int j, int k,

dp = q_arr(i+1-ishft,j,k,pres_comp) - q_arr(i-1-ishft,j,k,pres_comp);

denom = amrex::max(small_pres, std::abs(q_arr(i+2-ishft,j,k,pres_comp)-q_arr(i-2-ishft,j,k,pres_comp)));
denom = std::max(small_pres, std::abs(q_arr(i+2-ishft,j,k,pres_comp)-q_arr(i-2-ishft,j,k,pres_comp)));
zeta = std::abs(dp) / denom;
Real z2 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
Real z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

tst = 0.0_rt;
if (q_arr(i-1-ishft,j,k,QU) - q_arr(i+1-ishft,j,k,QU) >= 0.0_rt) {
tst = 1.0_rt;
}

tmp = amrex::min(q_arr(i+1-ishft,j,k,pres_comp), q_arr(i-1-ishft,j,k,pres_comp));
tmp = std::min(q_arr(i+1-ishft,j,k,pres_comp), q_arr(i-1-ishft,j,k,pres_comp));

Real chi2 = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
chi2 = tst;
}

Real flatn = 1.0_rt - amrex::max(chi2 * z2, chi * z);
Real flatn = 1.0_rt - std::max(chi2 * z2, chi * z);


#if AMREX_SPACEDIM >= 2
Expand All @@ -67,16 +67,16 @@ Real flatten(int i, int j, int k,

ishft = dp > 0.0_rt ? 1 : -1;

denom = amrex::max(small_pres, std::abs(q_arr(i,j+2,k,pres_comp) - q_arr(i,j-2,k,pres_comp)));
denom = std::max(small_pres, std::abs(q_arr(i,j+2,k,pres_comp) - q_arr(i,j-2,k,pres_comp)));
zeta = std::abs(dp) / denom;
z = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

tst = 0.0_rt;
if (q_arr(i,j-1,k,QV) - q_arr(i,j+1,k,QV) >= 0.0_rt) {
tst = 1.0_rt;
}

tmp = amrex::min(q_arr(i,j+1,k,pres_comp), q_arr(i,j-1,k,pres_comp));
tmp = std::min(q_arr(i,j+1,k,pres_comp), q_arr(i,j-1,k,pres_comp));

chi = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
Expand All @@ -86,23 +86,23 @@ Real flatten(int i, int j, int k,

dp = q_arr(i,j+1-ishft,k,pres_comp) - q_arr(i,j-1-ishft,k,pres_comp);

denom = amrex::max(small_pres, std::abs(q_arr(i,j+2-ishft,k,pres_comp) - q_arr(i,j-2-ishft,k,pres_comp)));
denom = std::max(small_pres, std::abs(q_arr(i,j+2-ishft,k,pres_comp) - q_arr(i,j-2-ishft,k,pres_comp)));
zeta = std::abs(dp) / denom;
z2 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

tst = 0.0_rt;
if (q_arr(i,j-1-ishft,k,QV) - q_arr(i,j+1-ishft,k,QV) >= 0.0_rt) {
tst = 1.0_rt;
}

tmp = amrex::min(q_arr(i,j+1-ishft,k,pres_comp), q_arr(i,j-1-ishft,k,pres_comp));
tmp = std::min(q_arr(i,j+1-ishft,k,pres_comp), q_arr(i,j-1-ishft,k,pres_comp));

chi2 = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
chi2 = tst;
}

flatn = amrex::min(flatn, 1.0_rt - amrex::max(chi2 * z2, chi * z));
flatn = std::min(flatn, 1.0_rt - std::max(chi2 * z2, chi * z));
#endif


Expand All @@ -113,16 +113,16 @@ Real flatten(int i, int j, int k,

ishft = dp > 0.0_rt ? 1: -1;

denom = amrex::max(small_pres, std::abs(q_arr(i,j,k+2,pres_comp) - q_arr(i,j,k-2,pres_comp)));
denom = std::max(small_pres, std::abs(q_arr(i,j,k+2,pres_comp) - q_arr(i,j,k-2,pres_comp)));
zeta = std::abs(dp) / denom;
z = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

tst = 0.0_rt;
if (q_arr(i,j,k-1,QW) - q_arr(i,j,k+1,QW) >= 0.0_rt) {
tst = 1.0_rt;
}

tmp = amrex::min(q_arr(i,j,k+1,pres_comp), q_arr(i,j,k-1,pres_comp));
tmp = std::min(q_arr(i,j,k+1,pres_comp), q_arr(i,j,k-1,pres_comp));

chi = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
Expand All @@ -132,23 +132,23 @@ Real flatten(int i, int j, int k,

dp = q_arr(i,j,k+1-ishft,pres_comp) - q_arr(i,j,k-1-ishft,pres_comp);

denom = amrex::max(small_pres, std::abs(q_arr(i,j,k+2-ishft,pres_comp) - q_arr(i,j,k-2-ishft,pres_comp)));
denom = std::max(small_pres, std::abs(q_arr(i,j,k+2-ishft,pres_comp) - q_arr(i,j,k-2-ishft,pres_comp)));
zeta = std::abs(dp) / denom;
z2 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1)));
z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt);

tst = 0.0_rt;
if (q_arr(i,j,k-1-ishft,QW) - q_arr(i,j,k+1-ishft,QW) >= 0.0_rt) {
tst = 1.0_rt;
}

tmp = amrex::min(q_arr(i,j,k+1-ishft,pres_comp), q_arr(i,j,k-1-ishft,pres_comp));
tmp = std::min(q_arr(i,j,k+1-ishft,pres_comp), q_arr(i,j,k-1-ishft,pres_comp));

chi2 = 0.0_rt;
if (std::abs(dp) > shktst*tmp) {
chi2 = tst;
}

flatn = amrex::min(flatn, 1.0_rt - amrex::max(chi2 * z2, chi * z));
flatn = std::min(flatn, 1.0_rt - std::max(chi2 * z2, chi * z));
#endif

return flatn;
Expand Down
6 changes: 3 additions & 3 deletions Source/hydro/fourth_order.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,7 @@ Castro::states(const Box& bx,
Real d3a_min = amrex::min(d3am1, d3a0, d3ap1, d3ap2);
Real d3a_max = amrex::max(d3am1, d3a0, d3ap1, d3ap2);

if (C3 * amrex::max(std::abs(d3a_min), std::abs(d3a_max)) <=
if (C3 * std::max(std::abs(d3a_min), std::abs(d3a_max)) <=
(d3a_max - d3a_min)) {
// limit
if (dafm*dafp < 0.0_rt) {
Expand Down Expand Up @@ -426,7 +426,7 @@ Castro::states(const Box& bx,
Real d3a_min = amrex::min(d3am1, d3a0, d3ap1, d3ap2);
Real d3a_max = amrex::max(d3am1, d3a0, d3ap1, d3ap2);

if (C3 * amrex::max(std::abs(d3a_min), std::abs(d3a_max)) <=
if (C3 * std::max(std::abs(d3a_min), std::abs(d3a_max)) <=
(d3a_max - d3a_min)) {
// limit
if (dafm*dafp < 0.0_rt) {
Expand Down Expand Up @@ -590,7 +590,7 @@ Castro::states(const Box& bx,
Real d3a_min = amrex::min(d3am1, d3a0, d3ap1, d3ap2);
Real d3a_max = amrex::max(d3am1, d3a0, d3ap1, d3ap2);

if (C3 * amrex::max(std::abs(d3a_min), std::abs(d3a_max)) <=
if (C3 * std::max(std::abs(d3a_min), std::abs(d3a_max)) <=
(d3a_max - d3a_min)) {
// limit
if (dafm*dafp < 0.0_rt) {
Expand Down
15 changes: 6 additions & 9 deletions Source/hydro/ppm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ ppm_reconstruct(const Real* s,
Real dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[i0] - s[im2]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
}

dsl = 2.0_rt * (s[i0] - s[im1]);
Expand All @@ -76,7 +76,7 @@ ppm_reconstruct(const Real* s,
Real dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl),std::abs(dsr));
}

// Interpolate s to edges
Expand All @@ -85,9 +85,7 @@ ppm_reconstruct(const Real* s,

// Make sure sedge lies in between adjacent cell-centered values

sm = amrex::max(sm, amrex::min(s[i0], s[im1]));
sm = amrex::min(sm, amrex::max(s[i0], s[im1]));

sm = std::clamp(sm, std::min(s[i0], s[im1]), std::max(s[i0], s[im1]));

// Compute van Leer slopes

Expand All @@ -97,7 +95,7 @@ ppm_reconstruct(const Real* s,
dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
}

dsl = 2.0_rt * (s[ip1] - s[i0]);
Expand All @@ -106,7 +104,7 @@ ppm_reconstruct(const Real* s,
dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip2] - s[i0]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
}

// Interpolate s to edges
Expand All @@ -115,8 +113,7 @@ ppm_reconstruct(const Real* s,

// Make sure sedge lies in between adjacent cell-centered values

sp = amrex::max(sp, amrex::min(s[ip1], s[i0]));
sp = amrex::min(sp, amrex::max(s[ip1], s[i0]));
sp = std::clamp(sp, std::min(s[ip1], s[i0]), std::max(s[ip1], s[i0]));


// Flatten the parabola
Expand Down
12 changes: 6 additions & 6 deletions Source/hydro/riemann_2shock_solvers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ namespace TwoShock {

// CG Eq. 31
gstar = (pstar-p)*gdot/(pstar+p) + gam;
gstar = amrex::max(gmin, amrex::min(gmax, gstar));
gstar = std::clamp(gstar, gmin, gmax);

// Now use that predicted value of game with the R-H jump conditions
// to compute the wave speed.
Expand Down Expand Up @@ -190,8 +190,8 @@ namespace TwoShock {

// these should consider a wider average of the cell-centered
// gammas
amrex::Real gmin = amrex::min(amrex::min(gamel, gamer), 1.0_rt);
amrex::Real gmax = amrex::max(amrex::max(gamel, gamer), 2.0_rt);
amrex::Real gmin = amrex::min(gamel, gamer, 1.0_rt);
amrex::Real gmax = amrex::max(gamel, gamer, 2.0_rt);

amrex::Real game_bar = 0.5_rt*(gamel + gamer);
amrex::Real gamc_bar = 0.5_rt*(ql.gamc + qr.gamc);
Expand Down Expand Up @@ -330,8 +330,8 @@ namespace TwoShock {
amrex::Real pstarl = 1.e200;
amrex::Real pstaru = -1.e200;
for (int n = riemann_shock_maxiter-6; n < riemann_shock_maxiter; n++) {
pstarl = amrex::min(pstarl, pstar_hist[n]);
pstaru = amrex::max(pstaru, pstar_hist[n]);
pstarl = std::min(pstarl, pstar_hist[n]);
pstaru = std::max(pstaru, pstar_hist[n]);
}

pstarl = amrex::max(pstarl, small_pres);
Expand Down Expand Up @@ -666,7 +666,7 @@ namespace TwoShock {

// interpolate for the case that we are in a rarefaction
amrex::Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt;
frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac));
frac = std::clamp(frac, 0.0_rt, 1.0_rt);

qint.rho = frac*rstar + (1.0_rt - frac)*ro;
qint.un = frac*ustar + (1.0_rt - frac)*uo;
Expand Down

0 comments on commit 8b289c1

Please sign in to comment.