JWL, and state-dependent EOS through Wood's law, characteristic BCs, probes and post-process output

examples/1D_jwl_release/case.py

@@ -0,0 +1,68 @@
+"""
+Isentropic release of a JWL fluid: a Riemann problem between two states of the same fluid.
+Everything left of the contact keeps the left state's entropy, so the rarefaction fan and the
+left star state must lie on the JWL isentrope through (rho0, p0), which is closed form.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D JWL isentropic release")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=400)
+parser.add_argument("--cfl", type=float, default=0.4)
+args = parser.parse_args()
+
+rho0, p0, rho_r, p_r = 1.0, 1.0, 0.3, 0.1
+jwl = {"a": 6.0, "b": 0.15, "r1": 4.0, "r2": 1.0, "omega": 0.3, "rho0": rho0}
+N, L, T_end = args.N, 1.0, 0.15
+c_max = math.sqrt(((1.0 + jwl["omega"]) * p0 + jwl["a"] + jwl["b"]) / rho0)  # generous bound on c + |u|
+dt = args.cfl * (L / N) / c_max
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+case = {
+    "run_time_info": "F",
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": N - 1,
+    "n": 0,
+    "p": 0,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": Nt,
+    "num_patches": 2,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "time_stepper": 3,
+    "recon_type": "weno",
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -3,
+    "bc_x%end": -3,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "F",
+    "fluid_pp(1)%eos": "jwl",
+    **{f"fluid_pp(1)%jwl_{k}": v for k, v in jwl.items()},
+}
+for pid, (x_c, rho, pres) in enumerate([(0.25, rho0, p0), (0.75, rho_r, p_r)], start=1):
+    case.update(
+        {
+            f"patch_icpp({pid})%geometry": 1,
+            f"patch_icpp({pid})%x_centroid": x_c,
+            f"patch_icpp({pid})%length_x": 0.5 * L,
+            f"patch_icpp({pid})%alpha_rho(1)": rho,
+            f"patch_icpp({pid})%alpha(1)": 1.0,
+            f"patch_icpp({pid})%vel(1)": 0.0,
+            f"patch_icpp({pid})%pres": pres,
+        }
+    )
+print(json.dumps(case))

examples/1D_mg_acoustic/case.py

@@ -0,0 +1,79 @@
+"""
+Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state.
+A rectangle patch carries a simple-wave perturbation (drho, dp = c^2 drho, du = c drho/rho0), so
+only the right-going characteristic is excited; the harness tracks the centroid of drho against
+the general analytic c. Patches rather than an analytic IC: an IC expression is compiled in, and
+every distinct one costs the test suite a full rebuild.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Mie-Gruneisen acoustic pulse")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=200)
+parser.add_argument("--cfl", type=float, default=0.4)
+args = parser.parse_args()
+
+rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4
+c = math.sqrt(c0**2 + (1.0 + gruneisen) * p0 / rho0)  # the frozen speed at the reference state
+amp, x0, width = 1.0e-4, 0.25, 0.1
+N, L, T_end = args.N, 1.0, 0.3
+dt = args.cfl * (L / N) / c
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "F",
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "m": N - 1,
+            "n": 0,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": Nt,
+            "t_step_save": Nt,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "time_stepper": 3,
+            "recon_type": "weno",
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "num_patches": 2,
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%alpha_rho(1)": rho0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": x0,
+            "patch_icpp(2)%length_x": width,
+            "patch_icpp(2)%alpha_rho(1)": rho0 + amp,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            "patch_icpp(2)%vel(1)": c / rho0 * amp,
+            "patch_icpp(2)%pres": p0 + c**2 * amp,
+            "fluid_pp(1)%eos": "mie_gruneisen",
+            "fluid_pp(1)%mg_rho0": rho0,
+            "fluid_pp(1)%mg_c0": c0,
+            "fluid_pp(1)%mg_s": s,
+            "fluid_pp(1)%mg_gruneisen": gruneisen,
+        }
+    )
+)

examples/1D_mg_impact/case.py

@@ -0,0 +1,71 @@
+"""
+Symmetric impact of two Mie-Gruneisen slabs approaching at relative speed U.
+Each slab is brought to rest by a shock with particle-velocity jump U/2, so the shock state
+lies exactly on the Hugoniot u_s = c0 + s u_p; the harness checks the shock speed and the
+plateau density against that relation.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Mie-Gruneisen symmetric impact")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=800)
+parser.add_argument("--U", type=float, default=1.0, help="closing speed of the two slabs")
+parser.add_argument("--cfl", type=float, default=0.4)
+args = parser.parse_args()
+
+rho0, p0, c0, s, gruneisen = 1.0, 1.0e-3, 1.0, 1.5, 0.4
+N, L, T_end = args.N, 1.0, 0.2
+dt = args.cfl * (L / N) / (c0 + (s + 1.0) * args.U)
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+case = {
+    "run_time_info": "F",
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": N - 1,
+    "n": 0,
+    "p": 0,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": Nt,
+    "num_patches": 2,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "time_stepper": 3,
+    "recon_type": "weno",
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -3,
+    "bc_x%end": -3,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "F",
+    "fluid_pp(1)%eos": "mie_gruneisen",
+    "fluid_pp(1)%mg_rho0": rho0,
+    "fluid_pp(1)%mg_c0": c0,
+    "fluid_pp(1)%mg_s": s,
+    "fluid_pp(1)%mg_gruneisen": gruneisen,
+}
+for pid, (x_c, vel) in enumerate([(0.25, 0.5 * args.U), (0.75, -0.5 * args.U)], start=1):
+    case.update(
+        {
+            f"patch_icpp({pid})%geometry": 1,
+            f"patch_icpp({pid})%x_centroid": x_c,
+            f"patch_icpp({pid})%length_x": 0.5 * L,
+            f"patch_icpp({pid})%alpha_rho(1)": rho0,
+            f"patch_icpp({pid})%alpha(1)": 1.0,
+            f"patch_icpp({pid})%vel(1)": vel,
+            f"patch_icpp({pid})%pres": p0,
+        }
+    )
+print(json.dumps(case))

src/common/m_constants.fpp

@@ -118,6 +118,8 @@ module m_constants
     !! cannot be auto-generated, so these are hand-written.
     integer, parameter :: eos_stiffened_gas = 1
     integer, parameter :: eos_ideal_gas = 2
+    integer, parameter :: eos_mie_gruneisen = 3
+    integer, parameter :: eos_jwl = 4
     integer, parameter :: num_synth_shells_max = 50  !< Max energy shells for synthetic turbulence
     integer, parameter :: num_turb_sources_max = 10  !< Max Gaussian forcing zones for synthetic turbulence
 

src/common/m_derived_types.fpp

@@ -394,6 +394,16 @@ module m_derived_types
         real(wp)               :: qvp            !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
         real(wp)               :: G
         integer                :: eos            !< Equation of state selector (eos_* in m_constants)
+        real(wp)               :: mg_rho0        !< Mie-Gruneisen reference density
+        real(wp)               :: mg_c0          !< Mie-Gruneisen bulk sound speed at mg_rho0
+        real(wp)               :: mg_s           !< Mie-Gruneisen linear Hugoniot slope, u_s = c0 + s u_p
+        real(wp)               :: mg_gruneisen   !< Gruneisen coefficient Gamma_G (not the shear modulus G)
+        real(wp)               :: jwl_a          !< JWL A
+        real(wp)               :: jwl_b          !< JWL B
+        real(wp)               :: jwl_r1         !< JWL R1
+        real(wp)               :: jwl_r2         !< JWL R2
+        real(wp)               :: jwl_omega      !< JWL omega (its Gruneisen coefficient)
+        real(wp)               :: jwl_rho0       !< JWL reference density
         logical                :: non_newtonian  !< Enable Herschel-Bulkley non-Newtonian viscosity
         real(wp)               :: K              !< HB consistency index
         real(wp)               :: nn             !< HB flow behavior index

src/common/m_global_parameters_common.fpp

@@ -54,6 +54,12 @@ module m_global_parameters_common
     !! written as p + B = const*rho**n.
     real(wp), allocatable, dimension(:) :: gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps
     $:GPU_DECLARE(create='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps]')
+    !> Per-fluid EOS selector and Mie-Gruneisen reference curve, resolved once at init like the arrays above
+    integer, allocatable, dimension(:) :: eoss
+    real(wp), allocatable, dimension(:) :: mg_rho0s, mg_c0s, mg_ss, mg_gruneisens
+    real(wp), allocatable, dimension(:) :: jwl_as, jwl_bs, jwl_r1s, jwl_r2s, jwl_omegas, jwl_rho0s
+    logical :: any_state_dependent_eos  !< True when some fluid's coefficients vary with density; set at init
+    $:GPU_DECLARE(create='[eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens, jwl_as, jwl_bs, jwl_r1s, jwl_r2s, jwl_omegas, jwl_rho0s, any_state_dependent_eos]')
     !> @}
 
     !> @name Fluids participating in shear and bulk viscosity