Skip to content

Commit c8d2f6a

Browse files
sbryngelsonsbryngelson
authored
Merge branch 'master' into flame-mixlyr-3D
2 parents e7da4f7 + 808df61 commit c8d2f6a

67 files changed

Lines changed: 4762 additions & 282 deletions

Some content is hidden

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

‎.claude/rules/common-pitfalls.md‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,19 @@ covered in `docs/documentation/contributing.md`.
6060
QBMM/viscous and MHD HLLD, while both Lagrange bubble cases *complete* with out-of-tolerance
6161
answers. Measured 2026-08-29 on MI210. A compile-only check returns green, so any future attempt to
6262
drop these must run the tests, not just build.
63+
- **CCE OpenACC (19.0.0 through 21.0.2, `-O2`; `-O0`/`-O1` correct; OpenMP offload unaffected): a device
64+
routine that contains any `GPU_LOOP` (itself or in anything it calls) must be called with scalars,
65+
never with an array element as an actual argument.** Every `routine` level is affected, including the
66+
conforming `loop vector` inside `routine vector`. With both ingredients present the
67+
element is misaddressed: an `intent(in)` element reads as garbage, an `intent(out)` element is
68+
never written. Either ingredient alone is fine, which is why master's `s_compute_pressure(q%sf(j,k,l),...)`
69+
works (no loop) and `s_compute_mixture_coefficients` works (scalar actuals). PR #1811 added the
70+
Newton and RK4 loops to the EOS helpers and every call that passed `%sf(j,k,l)` or `blkmod1(k,l,q)`
71+
ended in `NaN(s) in timestep output` on the Frontier CCE OpenACC lanes only, bit-identical on every
72+
other backend. Fix: copy elements to locals before the call, receive into a local. 37-line
73+
reproducer and the bisection: sbryngelson/compiler-bugs `cce/acc-routine-element-by-reference`,
74+
MFC #1815. Do not "fix" it by deleting the `seq` directives instead: they are the idiom master
75+
uses in every device routine.
6376
- The same "call it from the loop body" rule covers `m_thermochem`: calling `get_species_*` from
6477
inside a `GPU_ROUTINE` rather than from the kernel gave CCE OpenMP a runtime
6578
`Memory access fault by GPU node-N ... Reason: Unknown` on the first step (exit 134), while every

‎.github/workflows/common/build.sh‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,9 +41,16 @@ source .github/scripts/retry-build.sh
4141
# aborts in MPI_Init ("OPAL ERROR: Unreachable in file ext3x_client.c"). mpirun
4242
# is unaffected, and it is how MFC launches every binary anyway. Output is left
4343
# on stdout so a future failure is diagnosable from the CI log.
44+
#
45+
# --bind-to none is required, not cosmetic: on a busy GPU node SLURM hands the job
46+
# an offset/partial cpuset (e.g. cores 32-47), and Open MPI's default --bind-to
47+
# core then dies with "hwloc_set_cpubind Error for bitmap" before the binary even
48+
# starts. The run templates, coverage_build.py and preflight.sh all already pass
49+
# it; this smoke-test was the lone launch that did not, and it red-crossed jobs
50+
# that backfilled onto shared GPU nodes.
4451
validate_cmd=""
4552
if [ "$job_cluster" = "phoenix" ]; then
46-
validate_cmd='syscheck_bin=$(find build/install -name syscheck -type f 2>/dev/null | head -1); [ -z "$syscheck_bin" ] || mpirun -np 1 "$syscheck_bin"'
53+
validate_cmd='syscheck_bin=$(find build/install -name syscheck -type f 2>/dev/null | head -1); [ -z "$syscheck_bin" ] || mpirun --bind-to none -np 1 "$syscheck_bin"'
4754
fi
4855

4956
# --- Variant selection ---

‎docs/documentation/case.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -739,6 +739,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
739739
| `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database |
740740
| `pres_inf_wrt` | Logical | Add the liquid stiffness to the formatted database |
741741
| `c_wrt` | Logical | Add the sound speed to the database |
742+
| `T_wrt` | Logical | Add each fluid's temperature to the database (needs `cv` > 0) |
742743
| `omega_wrt(i)` | Logical | Add the $i$-direction vorticity to the database |
743744
| `schlieren_wrt` | Logical | Add the numerical schlieren to the database|
744745
| `qm_wrt` | Logical | Add the Q-criterion to the database|

‎docs/documentation/contributing.md‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -464,8 +464,8 @@ means supplying these, not grepping for `gammas`:
464464
| `f_pressure` / `s_compute_energy` | \f$p(e)\f$ and \f$E(p)\f$ |
465465
| `f_bulk_modulus` | \f$K(p)\f$ - every sound speed in MFC is \f$K/\rho\f$, differing only in how phases are mixed |
466466
| `s_compute_speed_of_sound` / `_avg` | that mixing: Wood's law, 6-equation, bubble-diluted |
467-
| `f_phase_internal_energy` | per-phase internal energy (6-equation model) |
468-
| `f_isentrope_exponent` / `f_isentrope_pressure` / `f_pressure_on_isentrope` | the isentrope \f$p + B = \textrm{const}\,\rho^n\f$ |
467+
| `s_phase_internal_energy` | per-phase internal energy (6-equation model) |
468+
| `f_isentrope_exponent` / `f_isentrope_pressure` | the isentrope \f$p + B = \textrm{const}\,\rho^n\f$ |
469469
| `f_sg_thermal` | the thermal law \f$p + B = (n-1)c_v\rho T\f$ |
470470

471471
The first six are *mechanical* - they need only \f$p, \rho, e, c\f$. The last two are *caloric* and

‎examples/1D_isentropic_release/case.py‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
"""
2+
Isentropic release of a fluid whose reference curve is an isentrope (JWL or Vinet): a Riemann problem
3+
between two states of the same fluid. Everything left of the contact keeps the left state's entropy,
4+
so the fan and the left star state must lie on the closed-form isentrope through (rho0, p0).
5+
"""
6+
7+
import argparse
8+
import json
9+
import math
10+
11+
parser = argparse.ArgumentParser(description="1D JWL isentropic release")
12+
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
13+
parser.add_argument("-N", type=int, default=400)
14+
parser.add_argument("--cfl", type=float, default=0.4)
15+
parser.add_argument("--eos", choices=["jwl", "vinet"], default="jwl", help="reference curve: both are isentropes")
16+
args = parser.parse_args()
17+
18+
rho0, p0, rho_r, p_r = 1.0, 1.0, 0.3, 0.1
19+
if args.eos == "jwl":
20+
fluid = {f"fluid_pp(1)%jwl_{k}": v for k, v in {"a": 6.0, "b": 0.15, "r1": 4.0, "r2": 1.0, "omega": 0.3, "rho0": rho0}.items()}
21+
else:
22+
fluid = {f"fluid_pp(1)%vinet_{k}": v for k, v in {"k0": 2.0, "k0p": 4.0, "gruneisen": 0.3, "rho0": rho0}.items()}
23+
N, L, T_end = args.N, 1.0, 0.15
24+
c_max = math.sqrt((1.3 * p0 + 6.15) / rho0) # generous bound on c + |u| for either fit
25+
dt = args.cfl * (L / N) / c_max
26+
Nt = math.ceil(T_end / dt)
27+
dt = T_end / Nt
28+
29+
case = {
30+
"run_time_info": "F",
31+
"x_domain%beg": 0.0,
32+
"x_domain%end": L,
33+
"m": N - 1,
34+
"n": 0,
35+
"p": 0,
36+
"dt": dt,
37+
"t_step_start": 0,
38+
"t_step_stop": Nt,
39+
"t_step_save": Nt,
40+
"num_patches": 2,
41+
"model_eqns": 2,
42+
"num_fluids": 1,
43+
"time_stepper": 3,
44+
"recon_type": "weno",
45+
"weno_order": 5,
46+
"weno_eps": 1.0e-16,
47+
"mapped_weno": "T",
48+
"riemann_solver": 2,
49+
"wave_speeds": 1,
50+
"avg_state": 2,
51+
"bc_x%beg": -3,
52+
"bc_x%end": -3,
53+
"format": 1,
54+
"precision": 2,
55+
"prim_vars_wrt": "T",
56+
"parallel_io": "F",
57+
"fluid_pp(1)%eos": args.eos,
58+
**fluid,
59+
}
60+
for pid, (x_c, rho, pres) in enumerate([(0.25, rho0, p0), (0.75, rho_r, p_r)], start=1):
61+
case.update(
62+
{
63+
f"patch_icpp({pid})%geometry": 1,
64+
f"patch_icpp({pid})%x_centroid": x_c,
65+
f"patch_icpp({pid})%length_x": 0.5 * L,
66+
f"patch_icpp({pid})%alpha_rho(1)": rho,
67+
f"patch_icpp({pid})%alpha(1)": 1.0,
68+
f"patch_icpp({pid})%vel(1)": 0.0,
69+
f"patch_icpp({pid})%pres": pres,
70+
}
71+
)
72+
print(json.dumps(case))

‎examples/1D_mg_acoustic/case.py‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
"""
2+
Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state.
3+
A rectangle patch carries a simple-wave perturbation (drho, dp = c^2 drho, du = c drho/rho0), so
4+
only the right-going characteristic is excited; the harness tracks the centroid of drho against
5+
the general analytic c. Patches rather than an analytic IC: an IC expression is compiled in, and
6+
every distinct one costs the test suite a full rebuild.
7+
"""
8+
9+
import argparse
10+
import json
11+
import math
12+
13+
parser = argparse.ArgumentParser(description="1D Mie-Gruneisen acoustic pulse")
14+
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
15+
parser.add_argument("-N", type=int, default=200)
16+
parser.add_argument("--cfl", type=float, default=0.4)
17+
parser.add_argument("--a", type=float, default=0.0, help="Gruneisen slope: Gamma_G = Gamma_0 + a mu")
18+
args = parser.parse_args()
19+
20+
rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4
21+
c = math.sqrt(c0**2 + (1.0 + gruneisen) * p0 / rho0 + args.a * p0 / (rho0 * gruneisen)) # the frozen speed at the reference state
22+
amp, x0, width = 1.0e-4, 0.25, 0.1
23+
N, L, T_end = args.N, 1.0, 0.3
24+
dt = args.cfl * (L / N) / c
25+
Nt = math.ceil(T_end / dt)
26+
dt = T_end / Nt
27+
28+
print(
29+
json.dumps(
30+
{
31+
"run_time_info": "F",
32+
"x_domain%beg": 0.0,
33+
"x_domain%end": L,
34+
"m": N - 1,
35+
"n": 0,
36+
"p": 0,
37+
"dt": dt,
38+
"t_step_start": 0,
39+
"t_step_stop": Nt,
40+
"t_step_save": Nt,
41+
"model_eqns": 2,
42+
"num_fluids": 1,
43+
"time_stepper": 3,
44+
"recon_type": "weno",
45+
"weno_order": 5,
46+
"weno_eps": 1.0e-16,
47+
"mapped_weno": "T",
48+
"riemann_solver": 2,
49+
"wave_speeds": 1,
50+
"avg_state": 2,
51+
"bc_x%beg": -3,
52+
"bc_x%end": -3,
53+
"format": 1,
54+
"precision": 2,
55+
"prim_vars_wrt": "T",
56+
"parallel_io": "F",
57+
"num_patches": 2,
58+
"patch_icpp(1)%geometry": 1,
59+
"patch_icpp(1)%x_centroid": 0.5,
60+
"patch_icpp(1)%length_x": L,
61+
"patch_icpp(1)%alpha_rho(1)": rho0,
62+
"patch_icpp(1)%alpha(1)": 1.0,
63+
"patch_icpp(1)%vel(1)": 0.0,
64+
"patch_icpp(1)%pres": p0,
65+
"patch_icpp(2)%geometry": 1,
66+
"patch_icpp(2)%alter_patch(1)": "T",
67+
"patch_icpp(2)%x_centroid": x0,
68+
"patch_icpp(2)%length_x": width,
69+
"patch_icpp(2)%alpha_rho(1)": rho0 + amp,
70+
"patch_icpp(2)%alpha(1)": 1.0,
71+
"patch_icpp(2)%vel(1)": c / rho0 * amp,
72+
"patch_icpp(2)%pres": p0 + c**2 * amp,
73+
"fluid_pp(1)%eos": "mie_gruneisen",
74+
"fluid_pp(1)%mg_rho0": rho0,
75+
"fluid_pp(1)%mg_c0": c0,
76+
"fluid_pp(1)%mg_s": s,
77+
"fluid_pp(1)%mg_gruneisen": gruneisen,
78+
**({"fluid_pp(1)%mg_gruneisen_a": args.a} if args.a else {}),
79+
}
80+
)
81+
)

‎examples/1D_mg_impact/case.py‎

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
"""
2+
Symmetric impact of two Mie-Gruneisen slabs approaching at relative speed U.
3+
Each slab is brought to rest by a shock with particle-velocity jump U/2, so the shock state
4+
lies exactly on the Hugoniot u_s = c0 + s u_p; the harness checks the shock speed and the
5+
plateau density against that relation.
6+
"""
7+
8+
import argparse
9+
import json
10+
import math
11+
12+
parser = argparse.ArgumentParser(description="1D Mie-Gruneisen symmetric impact")
13+
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
14+
parser.add_argument("-N", type=int, default=800)
15+
parser.add_argument("--U", type=float, default=1.0, help="closing speed of the two slabs")
16+
parser.add_argument("--cfl", type=float, default=0.4)
17+
parser.add_argument("--s2", type=float, default=0.0, help="quadratic Hugoniot coefficient u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3")
18+
parser.add_argument("--s3", type=float, default=0.0, help="cubic Hugoniot coefficient")
19+
args = parser.parse_args()
20+
21+
rho0, p0, c0, s, gruneisen = 1.0, 1.0e-3, 1.0, 1.5, 0.4
22+
N, L, T_end = args.N, 1.0, 0.2
23+
dt = args.cfl * (L / N) / (c0 + (s + 1.0 + (args.s2 + args.s3 * args.U) * args.U) * args.U)
24+
Nt = math.ceil(T_end / dt)
25+
dt = T_end / Nt
26+
27+
case = {
28+
"run_time_info": "F",
29+
"x_domain%beg": 0.0,
30+
"x_domain%end": L,
31+
"m": N - 1,
32+
"n": 0,
33+
"p": 0,
34+
"dt": dt,
35+
"t_step_start": 0,
36+
"t_step_stop": Nt,
37+
"t_step_save": Nt,
38+
"num_patches": 2,
39+
"model_eqns": 2,
40+
"num_fluids": 1,
41+
"time_stepper": 3,
42+
"recon_type": "weno",
43+
"weno_order": 5,
44+
"weno_eps": 1.0e-16,
45+
"mapped_weno": "T",
46+
"riemann_solver": 2,
47+
"wave_speeds": 1,
48+
"avg_state": 2,
49+
"bc_x%beg": -3,
50+
"bc_x%end": -3,
51+
"format": 1,
52+
"precision": 2,
53+
"prim_vars_wrt": "T",
54+
"parallel_io": "F",
55+
"fluid_pp(1)%eos": "mie_gruneisen",
56+
"fluid_pp(1)%mg_rho0": rho0,
57+
"fluid_pp(1)%mg_c0": c0,
58+
"fluid_pp(1)%mg_s": s,
59+
"fluid_pp(1)%mg_gruneisen": gruneisen,
60+
**({"fluid_pp(1)%mg_s2": args.s2} if args.s2 else {}),
61+
**({"fluid_pp(1)%mg_s3": args.s3} if args.s3 else {}),
62+
}
63+
for pid, (x_c, vel) in enumerate([(0.25, 0.5 * args.U), (0.75, -0.5 * args.U)], start=1):
64+
case.update(
65+
{
66+
f"patch_icpp({pid})%geometry": 1,
67+
f"patch_icpp({pid})%x_centroid": x_c,
68+
f"patch_icpp({pid})%length_x": 0.5 * L,
69+
f"patch_icpp({pid})%alpha_rho(1)": rho0,
70+
f"patch_icpp({pid})%alpha(1)": 1.0,
71+
f"patch_icpp({pid})%vel(1)": vel,
72+
f"patch_icpp({pid})%pres": p0,
73+
}
74+
)
75+
print(json.dumps(case))

‎src/common/m_constants.fpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,11 @@ module m_constants
118118
!! cannot be auto-generated, so these are hand-written.
119119
integer, parameter :: eos_stiffened_gas = 1
120120
integer, parameter :: eos_ideal_gas = 2
121+
integer, parameter :: eos_mie_gruneisen = 3
122+
integer, parameter :: eos_jwl = 4
123+
integer, parameter :: eos_vinet = 5
124+
integer, parameter :: eos_rk4_steps = 8 !< fixed-step RK4 along a phasic isentrope or a reference temperature
125+
integer, parameter :: ode_isentrope = 1, ode_reference_temperature = 2 !< the two ODEs s_rk4 integrates
121126
integer, parameter :: num_synth_shells_max = 50 !< Max energy shells for synthetic turbulence
122127
integer, parameter :: num_turb_sources_max = 10 !< Max Gaussian forcing zones for synthetic turbulence
123128

‎src/common/m_derived_types.fpp‎

Lines changed: 47 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -386,22 +386,42 @@ module m_derived_types
386386
!> Derived type annexing the physical parameters (PP) of the fluids. These include the specific heat ratio function and liquid
387387
!! stiffness function.
388388
type physical_parameters
389-
real(wp) :: gamma !< Sp. heat ratio
390-
real(wp) :: pi_inf !< Liquid stiffness
391-
real(wp), dimension(2) :: Re !< Reynolds number
392-
real(wp) :: cv !< heat capacity
393-
real(wp) :: qv !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004))
394-
real(wp) :: qvp !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
389+
real(wp) :: gamma !< Sp. heat ratio
390+
real(wp) :: pi_inf !< Liquid stiffness
391+
real(wp), dimension(2) :: Re !< Reynolds number
392+
real(wp) :: cv !< heat capacity
393+
real(wp) :: qv !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004))
394+
real(wp) :: qvp !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
395395
real(wp) :: G
396-
integer :: eos !< Equation of state selector (eos_* in m_constants)
397-
logical :: non_newtonian !< Enable Herschel-Bulkley non-Newtonian viscosity
398-
real(wp) :: K !< HB consistency index
399-
real(wp) :: nn !< HB flow behavior index
400-
real(wp) :: tau0 !< HB yield stress (0 => power-law)
401-
real(wp) :: hb_m !< Papanastasiou regularization parameter
402-
real(wp) :: mu_min !< Lower viscosity clamp (inactive sentinel = dflt_real)
403-
real(wp) :: mu_max !< Upper viscosity clamp (required when non_newtonian)
404-
real(wp) :: mu_bulk !< Bulk viscosity for NN (inactive sentinel = dflt_real)
396+
integer :: eos !< Equation of state selector (eos_* in m_constants)
397+
real(wp) :: mg_rho0 !< Mie-Gruneisen reference density
398+
real(wp) :: mg_c0 !< Mie-Gruneisen bulk sound speed at mg_rho0
399+
real(wp) :: mg_s !< Mie-Gruneisen linear Hugoniot slope, u_s = c0 + s u_p
400+
real(wp) :: mg_gruneisen !< Gruneisen coefficient Gamma_G (not the shear modulus G)
401+
real(wp) :: mg_gruneisen_a !< d(Gamma_G)/d(mu): Gamma_G = Gamma_0 + a mu, zero keeps it constant
402+
real(wp) :: mg_t0 !< temperature at the reference density (for T output)
403+
real(wp) :: mg_s2, mg_s3 !< u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3; zero keeps the fit linear
404+
real(wp) :: jwl_a !< JWL A
405+
real(wp) :: jwl_b !< JWL B
406+
real(wp) :: jwl_r1 !< JWL R1
407+
real(wp) :: jwl_r2 !< JWL R2
408+
real(wp) :: jwl_omega !< JWL omega (its Gruneisen coefficient)
409+
real(wp) :: jwl_rho0 !< JWL reference density
410+
real(wp) :: jwl_t0 !< temperature at the reference density (for T output)
411+
real(wp) :: vinet_k0 !< Vinet bulk modulus at rho0
412+
real(wp) :: vinet_k0p !< Vinet pressure derivative of the bulk modulus
413+
real(wp) :: vinet_rho0 !< Vinet reference density
414+
real(wp) :: vinet_gruneisen !< Gruneisen coefficient at rho0
415+
real(wp) :: vinet_gruneisen_a !< d(Gamma_G)/d(mu)
416+
real(wp) :: vinet_t0 !< temperature at the reference density (for T output)
417+
logical :: non_newtonian !< Enable Herschel-Bulkley non-Newtonian viscosity
418+
real(wp) :: K !< HB consistency index
419+
real(wp) :: nn !< HB flow behavior index
420+
real(wp) :: tau0 !< HB yield stress (0 => power-law)
421+
real(wp) :: hb_m !< Papanastasiou regularization parameter
422+
real(wp) :: mu_min !< Lower viscosity clamp (inactive sentinel = dflt_real)
423+
real(wp) :: mu_max !< Upper viscosity clamp (required when non_newtonian)
424+
real(wp) :: mu_bulk !< Bulk viscosity for NN (inactive sentinel = dflt_real)
405425
end type physical_parameters
406426

407427
!> Derived type annexing the physical parameters required for sub-grid bubble models
@@ -534,6 +554,18 @@ module m_derived_types
534554
integer :: substeps !< Operator-split sub-steps per time step (0 = source added to the flow RHS)
535555
end type reactive_burn_parameters
536556
557+
!> Coefficients of one fluid's equation of state, resolved once at init. Held as a record per fluid rather than as parallel
558+
!! arrays: every read wants several of these for a single fluid, so one base address serves them all, where fifteen arrays cost
559+
!! fifteen live descriptors in the Riemann kernels.
560+
type eos_coefficients
561+
real(wp) :: rho0, t0 !< Reference density [kg/m^3] and temperature [K]
562+
real(wp) :: gruneisen0, gruneisen_a !< Gruneisen closure Gamma_G = Gamma_0 + a mu
563+
real(wp) :: c0, s, s2, s3 !< Mie-Gruneisen Hugoniot u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3
564+
real(wp) :: mu_max !< Compression at which a cubic Hugoniot fit turns over
565+
real(wp) :: a, b, r1, r2 !< JWL principal isentrope p = A exp(-R1 V) + B exp(-R2 V)
566+
real(wp) :: k0, k0p !< Vinet bulk modulus and its pressure derivative
567+
end type eos_coefficients
568+
537569
!> Lagrangian bubble parameters
538570
type bubbles_lagrange_parameters
539571

0 commit comments

Comments
 (0)