Evaluate Mie-Gruneisen coefficients per cell along the 5-equation path

examples/1D_mg_acoustic/case.py

@@ -0,0 +1,70 @@
+"""
+Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state.
+The pulse is a simple wave (drho, dp = c^2 drho, du = c drho/rho0) so only the right-going
+characteristic carries it; the harness measures its speed against the general analytic c.
+"""
+
+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.3, 0.05
+N, L, T_end = args.N, 1.0, 0.4
+dt = args.cfl * (L / N) / c
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+pulse = f"{amp}*exp(-((x - {x0})/{width})**2)"
+
+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,
+            "num_patches": 1,
+            "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",
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%alpha_rho(1)": f"{rho0} + {pulse}",
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "patch_icpp(1)%vel(1)": f"{c}/{rho0}*{pulse}",
+            "patch_icpp(1)%pres": f"{p0} + {c}**2*{pulse}",
+            "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,7 @@ 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 :: 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,10 @@ 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)
         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,11 @@ 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
+    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, any_state_dependent_eos]')
     !> @}
 
     !> @name Fluids participating in shear and bulk viscosity

src/common/m_variables_conversion.fpp

@@ -28,8 +28,8 @@ module m_variables_conversion
         & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, f_bulk_modulus, &
         & f_pressure, f_phase_internal_energy, f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, f_pressure_on_isentrope, &
         & s_compute_mixture_coefficients_dt, s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, f_elastic_energy, &
-        & f_hypoelastic_energy, f_relativistic_enthalpy, s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, &
-        & isentrope_B, cvs, qvs, qvps
+        & f_hypoelastic_energy, f_relativistic_enthalpy, s_eos_coefficients, s_finalize_variables_conversion_module, gammas, &
+        & isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps
 
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
@@ -267,6 +267,8 @@ contains
         $:GPU_UPDATE(device='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]')
 
         @:ALLOCATE(gammas (1:num_fluids))
+        @:ALLOCATE(eoss (1:num_fluids), mg_rho0s (1:num_fluids), mg_c0s (1:num_fluids), mg_ss (1:num_fluids), &
+                   & mg_gruneisens (1:num_fluids))
         @:ALLOCATE(isentrope_n (1:num_fluids))
         @:ALLOCATE(pi_infs(1:num_fluids))
         @:ALLOCATE(isentrope_B(1:num_fluids))
@@ -275,6 +277,7 @@ contains
         @:ALLOCATE(qvps    (1:num_fluids))
         @:ALLOCATE(Gs_vc     (1:num_fluids))
 
+        any_state_dependent_eos = .false.
         do i = 1, num_fluids
             gammas(i) = fluid_pp(i)%gamma
             isentrope_n(i) = f_isentrope_exponent(gammas(i))
@@ -292,8 +295,15 @@ contains
             cvs(i) = fluid_pp(i)%cv
             qvs(i) = fluid_pp(i)%qv
             qvps(i) = fluid_pp(i)%qvp
+            eoss(i) = fluid_pp(i)%eos
+            mg_rho0s(i) = fluid_pp(i)%mg_rho0
+            mg_c0s(i) = fluid_pp(i)%mg_c0
+            mg_ss(i) = fluid_pp(i)%mg_s
+            mg_gruneisens(i) = fluid_pp(i)%mg_gruneisen
+            if (fluid_pp(i)%eos == eos_mie_gruneisen) any_state_dependent_eos = .true.
         end do
-        $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc]')
+        $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, &
+                     & mg_gruneisens, any_state_dependent_eos]')
 
         @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max)))
         Res_vc = dflt_real
@@ -1173,7 +1183,7 @@ contains
 
         if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf)
 
-        @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc)
+        @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens)
         if (allocated(bubrs_vc)) then
             @:DEALLOCATE(bubrs_vc)
         end if
@@ -1196,6 +1206,7 @@ contains
             real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
         #:endif
         real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
+        real(wp)              :: gamma_i, pi_inf_i, dpi_i
         integer               :: i  !< Loop iterator over fluids
 
         ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
@@ -1217,8 +1228,14 @@ contains
             $:GPU_LOOP(parallelism='[seq]')
             do i = 1, num_fluids
                 rho_K = rho_K + alpha_rho_K(i)
-                gamma_K = gamma_K + alpha_K(i)*gammas(i)
-                pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i)
+                if (any_state_dependent_eos) then
+                    call s_eos_coefficients(alpha_rho_K(i)/max(alpha_K(i), sgm_eps), i, gamma_i, pi_inf_i, dpi_i)
+                else
+                    gamma_i = gammas(i)
+                    pi_inf_i = pi_infs(i)
+                end if
+                gamma_K = gamma_K + alpha_K(i)*gamma_i
+                pi_inf_K = pi_inf_K + alpha_K(i)*pi_inf_i
                 qv_K = qv_K + alpha_rho_K(i)*qvs(i)
             end do
         end if
@@ -1284,6 +1301,49 @@ contains
 
     end subroutine s_compute_energy
 
+    !> Coefficients of fluid i at density rho in the form rho e = Gamma p + Pi that every operator here consumes, with dPi/drho. Any
+    !! Mie-Gruneisen EOS p = p_ref + rho Gamma_G (e - e_ref) is this form with Gamma = 1/Gamma_G and Pi = rho e_ref - p_ref/Gamma_G;
+    !! a new family adds one case supplying its reference curve. Gamma_G is constant, so dGamma/drho = 0. Stiffened and ideal gas
+    !! keep the constants resolved at init, bit for bit.
+    subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in)  :: rho
+        integer, intent(in)   :: i
+        real(wp), intent(out) :: gamma, pi_inf, dpi
+        real(wp)              :: mu, d, p_ref, e_ref, dp_dmu, de_dmu, G0
+
+        select case (eoss(i))
+        case (eos_mie_gruneisen)
+            ! Linear-Hugoniot reference curve, u_s = c0 + s u_p: p_H = rho0 c0^2 mu (1 + mu)/(1 - (s - 1) mu)^2 on
+            ! compression, extended linearly on release, with the Hugoniot energy e_H = p_H mu/(2 rho0 (1 + mu)).
+            ! Pole at mu = 1/(s - 1); the validator warns when the initial state is near it.
+            mu = rho/mg_rho0s(i) - 1._wp
+            if (mu >= 0._wp) then
+                d = 1._wp - (mg_ss(i) - 1._wp)*mu
+                p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu*(1._wp + mu)/(d*d)
+                dp_dmu = mg_rho0s(i)*mg_c0s(i)**2*((1._wp + 2._wp*mu)*d + 2._wp*(mg_ss(i) - 1._wp)*mu*(1._wp + mu))/(d*d*d)
+            else
+                p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu
+                dp_dmu = mg_rho0s(i)*mg_c0s(i)**2
+            end if
+            e_ref = p_ref*mu/(2._wp*mg_rho0s(i)*(1._wp + mu))
+            de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*mg_rho0s(i)*(1._wp + mu)**2)
+            G0 = mg_gruneisens(i)
+        case default
+            gamma = gammas(i)
+            pi_inf = pi_infs(i)
+            dpi = 0._wp
+            return
+        end select
+
+        gamma = 1._wp/G0
+        pi_inf = rho*e_ref - p_ref/G0
+        dpi = e_ref + (rho*de_dmu - dp_dmu/G0)/mg_rho0s(i)  ! d/drho = (1/rho0) d/dmu
+
+    end subroutine s_eos_coefficients
+
     !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n.
     function f_isentrope_exponent(gamma) result(n)
 
@@ -1420,7 +1480,7 @@ contains
 
     !> Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H, |u|^2 and qv all cancel out of c^2
     !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant.
-    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
@@ -1431,8 +1491,14 @@ contains
             real(wp), dimension(num_fluids), intent(in) :: adv
         #:endif
         real(wp), intent(out) :: c
-        real(wp)              :: alf  !< Subgrid void fraction; dilute by construction
-        integer               :: q
+        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+            real(wp), dimension(3), intent(in), optional :: alpha_rho
+        #:else
+            real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
+        #:endif
+        real(wp) :: alf  !< Subgrid void fraction; dilute by construction
+        real(wp) :: rho_q, gamma_q, pi_inf_q, dpi_q
+        integer  :: q
 
         if (chemistry) then  ! Reacting mixture sound speed
             c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
@@ -1441,7 +1507,16 @@ contains
         else
             ! Every case below is a bulk modulus over a density. The equation of state enters
             ! only through f_bulk_modulus; the cases differ in how the phases are mixed.
-            if (alt_soundspeed) then  ! Wood's law: volume-weighted harmonic mean
+            if (any_state_dependent_eos .and. present(alpha_rho)) then  ! frozen mixing: each phase's modulus at its own density
+                c = 0._wp
+                $:GPU_LOOP(parallelism='[seq]')
+                do q = 1, num_fluids
+                    rho_q = alpha_rho(q)/max(adv(q), sgm_eps)
+                    call s_eos_coefficients(rho_q, q, gamma_q, pi_inf_q, dpi_q)
+                    c = c + adv(q)*(f_bulk_modulus(pres, gamma_q, pi_inf_q) - rho_q*dpi_q/gamma_q)
+                end do
+                c = c/rho
+            else if (alt_soundspeed) then  ! Wood's law: volume-weighted harmonic mean
                 c = 1._wp/(rho*(adv(1)/f_bulk_modulus(pres, gammas(1), pi_infs(1)) + adv(2)/f_bulk_modulus(pres, gammas(2), &
                            & pi_infs(2))))
             else if (model_eqns == model_eqns_6eq) then  ! volume-weighted arithmetic mean
@@ -1475,7 +1550,7 @@ contains
     !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its
     !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from
     !! s_compute_speed_of_sound; keep the condition below in step with the branch list there.
-    subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c)
+    subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
@@ -1486,17 +1561,23 @@ contains
             real(wp), dimension(num_fluids), intent(in) :: adv
         #:endif
         real(wp), intent(out) :: c
+        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+            real(wp), dimension(3), intent(in), optional :: alpha_rho
+        #:else
+            real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
+        #:endif
 
         if (chemistry) then  ! Reacting mixture sound speed
             if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
                 c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
             else
-                call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+                call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
             end if
         else if (relativity) then  ! Relativistic sound speed
             c = sqrt((1._wp + 1._wp/gamma)*pres/rho/H)
-        else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler)) then
-            call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+        else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) &
+                 & .or. any_state_dependent_eos) then
+            call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
         else  ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
             c = (H - 5.e-1*vel_sum - qv/rho)/gamma
 

src/post_process/m_data_output.fpp

@@ -1281,7 +1281,7 @@ contains
 
                     call s_compute_mixture_coefficients(alpha_rho, adv, rho, gamma, pi_inf, qv)
 
-                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
+                    call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
 
                     Ma = maxvel/c
                     if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then