SCC26 Practice Problem Prep

examples/2D_shockdroplet/case.py

@@ -1,30 +1,70 @@
 #!/usr/bin/env python3
+import argparse
 import json
 import math
 
-Ma = 1.4
-ps = 238558
-rho_post_a = 2.18
+parser = argparse.ArgumentParser(prog="2D_shockdroplet", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument(
+    "--mfc",
+    type=json.loads,
+    default="{}",
+    metavar="DICT",
+    help="MFC's toolchain's internal state.",
+)
+parser.add_argument("--mach", type=float, default=1.4, help="Shock Mach number.")
+parser.add_argument("--adap-dt", action="store_true", help="Use CFL-based adaptive time-stepping.")
+parser.add_argument("--cfl", type=float, default=0.2, help="CFL number.")
+parser.add_argument("--tend", type=float, default=1.0, help="Dimensionless end time.")
+parser.add_argument("--res", type=float, default=1.0, help="Resolution scaling factor.")
+args = parser.parse_args()
+
+Ma = args.mach
+p_a = 101325.0
 rho_a = 1.204
 rho_w = 1000
 gam_a = 1.4
 gam_w = 6.12
 pi_w = 3.43e8
-vel = 226
+
+# Post-shock state from the Rankine-Hugoniot conditions (shock moving into still air)
+ps = p_a * (1 + 2 * gam_a / (gam_a + 1) * (Ma**2 - 1))
+rho_post_a = rho_a * (gam_a + 1) * Ma**2 / ((gam_a - 1) * Ma**2 + 2)
+c_a = math.sqrt(gam_a * p_a / rho_a)
+vel = 2 * c_a / (gam_a + 1) * (Ma**2 - 1) / Ma
+
 rho = 1
 c_l = math.sqrt(1.4 * ps / rho)
 eps = 1e-9
 
 D = 0.048
-Ny = 299.0
-Nx = 1199.0
-dx = 0.25 / Nx  # 8.3e-6
+Ny = (300 * args.res) - 1
+Nx = (1200 * args.res) - 1
+dx = 0.25 / Nx
+
+# End time in units of the Ranger & Nicholls (1969) breakup time scale
+# t_breakup = D sqrt(rho_l/rho_g) / u_g, with post-shock gas conditions
+tau_end = args.tend
+time_end = tau_end * D * math.sqrt(rho_w / rho_post_a) / vel
+cfl = args.cfl
 
-time_end = 0.005  # 50us
-cfl = 0.25
+dt = cfl * dx / c_l
+Nt = int(time_end / dt)
 
-dt = cfl * dx / c_l  # 5.3E-9
-Nt = int(time_end / dt)  # 10000
+if args.adap_dt:
+    time_stepping = {
+        "cfl_adap_dt": "T",
+        "cfl_target": cfl,
+        "n_start": 0,
+        "t_stop": time_end,
+        "t_save": time_end / 100,
+    }
+else:
+    time_stepping = {
+        "dt": dt,
+        "t_step_start": 0,
+        "t_step_stop": Nt,
+        "t_step_save": math.ceil(Nt / 100),
+    }
 
 print(
     json.dumps(
@@ -44,10 +84,7 @@
             "m": int(Nx),
             "n": int(Ny),
             "p": 0,
-            "dt": dt,
-            "t_step_start": 0,
-            "t_step_stop": Nt,
-            "t_step_save": Nt,  # math.ceil(Nt/100),
+            **time_stepping,
             # Simulation Algorithm Parameters
             "num_patches": 3,
             "model_eqns": "5eq",
@@ -83,7 +120,7 @@
             "patch_icpp(1)%length_y": 14 * D,
             "patch_icpp(1)%vel(1)": 0.0,
             "patch_icpp(1)%vel(2)": 0.0e00,
-            "patch_icpp(1)%pres": 101325.0,
+            "patch_icpp(1)%pres": p_a,
             "patch_icpp(1)%alpha_rho(1)": eps * 1000,
             "patch_icpp(1)%alpha_rho(2)": (1 - eps) * 1.17,
             "patch_icpp(1)%alpha(1)": eps,
@@ -110,7 +147,7 @@
             "patch_icpp(3)%alter_patch(1)": "T",
             "patch_icpp(3)%vel(1)": 0.0,
             "patch_icpp(3)%vel(2)": 0.0e00,
-            "patch_icpp(3)%pres": 101325.0,
+            "patch_icpp(3)%pres": p_a,
             "patch_icpp(3)%alpha_rho(1)": (1 - eps) * rho_w,
             "patch_icpp(3)%alpha_rho(2)": eps * 1.17,
             "patch_icpp(3)%alpha(1)": 1 - eps,  # 0.95

src/common/m_variables_conversion.fpp

@@ -424,6 +424,8 @@ contains
     !! each model has different variable sets and EOS.
     subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, q_T_sf, qK_prim_vf, ibounds)
 
+        use m_global_parameters_common, only: shear_indices  ! Performance fix with AMDFlang
+
         type(scalar_field), dimension(sys_size), intent(in)    :: qK_cons_vf
         type(scalar_field), intent(inout)                      :: q_T_sf
         type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf