Three new examples cases

examples/2D_kelvin_helmholtz/README.md

@@ -0,0 +1,10 @@
+# Kelvin-Helmholtz Instability (2D)
+
+Reference: See Example 4.8.
+> A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021).
+
+### Initial State 
+<img src="figure0.png" height="MAX_HEIGHT"/>
+
+### Evolved State  
+<img src="figure1.png" height="MAX_HEIGHT"/>

examples/2D_kelvin_helmholtz/case.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+import json
+import math
+
+
+eps = 1e-6
+time_end = 1.0
+time_save = time_end / 100.0
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0,
+            "y_domain%end": 1.0,
+            "m": 512,
+            "n": 512,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.2,
+            "n_start": 0,
+            "t_stop": time_end,
+            "t_save": time_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-10,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -1,
+            "bc_y%end": -1,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Background
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%hcid": 207,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%length_y": 1.0,
+            "patch_icpp(1)%vel(1)": -0.5,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 2.5,
+            "patch_icpp(1)%alpha_rho(1)": 1.0 - eps,
+            "patch_icpp(1)%alpha(1)": 1.0 - eps,
+            "patch_icpp(1)%alpha_rho(2)": eps,
+            "patch_icpp(1)%alpha(2)": eps,
+            # Center Strip (0.25 < y <= 0.75)
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%hcid": 207,
+            "patch_icpp(2)%x_centroid": 0.5,
+            "patch_icpp(2)%y_centroid": 0.5,
+            "patch_icpp(2)%length_x": 1.0,
+            "patch_icpp(2)%length_y": 0.5,
+            "patch_icpp(2)%vel(1)": 0.5,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": 2.5,
+            "patch_icpp(2)%alpha_rho(1)": 2.0 * eps,
+            "patch_icpp(2)%alpha(1)": eps,
+            "patch_icpp(2)%alpha_rho(2)": 2.0 * (1.0 - eps),
+            "patch_icpp(2)%alpha(2)": 1.0 - eps,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0e00,
+            "fluid_pp(2)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 0.0e00,
+        }
+    )
+)

examples/2D_richtmyer_meshkov/README.md

@@ -0,0 +1,10 @@
+# Richtmyer-Meshkov Instability (2D)
+
+Reference: See Example 4.18.
+> A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021).
+
+### Initial State
+<img src="figure0.png" height="MAX_HEIGHT"/>
+
+### Evolved State
+<img src="figure1.png" height="MAX_HEIGHT"/>

examples/2D_richtmyer_meshkov/case.py

@@ -0,0 +1,98 @@
+#!/usr/bin/env python3
+import json
+import math
+
+
+mu = 1.0e-4
+lambd = 1.0
+time_end = 15.0
+time_save = time_end / 20.0
+eps = 1.0e-6
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 16.0 * lambd,
+            "y_domain%beg": 0,
+            "y_domain%end": lambd,
+            "m": 4096,
+            "n": 256,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.1,
+            "n_start": 0,
+            "t_stop": time_end,
+            "t_save": time_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "recon_type": 1,
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-9,
+            "null_weights": "F",
+            "mapped_weno": "T",
+            "mp_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -17,
+            "bc_x%end": -17,
+            "bc_y%beg": -15,
+            "bc_y%end": -15,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Fluid #1 = Heavier Fluid
+            # Fluid #2 = Lighter Fluid
+            # Pre Shock
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%hcid": 208,
+            "patch_icpp(1)%x_centroid": 0.5 * 0.7 * lambd,
+            "patch_icpp(1)%y_centroid": 0.5 * lambd,
+            "patch_icpp(1)%length_x": 0.7 * lambd,
+            "patch_icpp(1)%length_y": lambd,
+            "patch_icpp(1)%vel(1)": 1.24,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 1.0 / 1.4,
+            "patch_icpp(1)%alpha_rho(1)": 1.0 * eps,
+            "patch_icpp(1)%alpha(1)": eps,
+            "patch_icpp(1)%alpha_rho(2)": 1.0 * (1.0 - eps),
+            "patch_icpp(1)%alpha(2)": (1.0 - eps),
+            # Post Shock
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%hcid": 208,
+            "patch_icpp(2)%x_centroid": 0.7 * lambd + 0.5 * 15.3 * lambd,
+            "patch_icpp(2)%y_centroid": 0.5 * lambd,
+            "patch_icpp(2)%length_x": 15.3 * lambd,
+            "patch_icpp(2)%length_y": lambd,
+            "patch_icpp(2)%vel(1)": 0.8787,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": 1.6272 / 1.4,
+            "patch_icpp(2)%alpha_rho(1)": 1.4112 * eps,
+            "patch_icpp(2)%alpha(1)": eps,
+            "patch_icpp(2)%alpha_rho(2)": 1.4112 * (1.0 - eps),
+            "patch_icpp(2)%alpha(2)": (1.0 - eps),
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.093 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0e00,
+            "fluid_pp(2)%gamma": 1.0e00 / (1.4 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 0.0e00,
+            "viscous": "T",
+            "fluid_pp(1)%Re(1)": 1 / mu,
+            "fluid_pp(2)%Re(1)": 1 / mu,
+
+        }
+    )
+)

examples/2D_viscous_shock_tube/README.md

@@ -0,0 +1,10 @@
+# Viscous Shock Tube (2D)
+
+Reference: See Example 4.13.
+> A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021)., see Example 4.13
+
+### Initial State 
+<img src="figure0.png" height="MAX_HEIGHT"/>
+
+### Evolved State  
+<img src="figure1.png" height="MAX_HEIGHT"/>

examples/2D_viscous_shock_tube/case.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+import json
+import math
+
+
+mu = 1.0e-3
+time_end = 1.2
+time_save = time_end / 100.0
+eps = 1.0e-6
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0,
+            "y_domain%end": 0.5,
+            "m": 1280,
+            "n": 640,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.1,
+            "n_start": 0,
+            "t_stop": time_end,
+            "t_save": time_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-10,
+            "null_weights": "F",
+            "mapped_weno": "T",
+            "mp_weno": "T",
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -8,
+            "bc_x%end": -15,
+            "bc_y%beg": -16,
+            "bc_y%end": -15,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Left
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5 * 0.5,
+            "patch_icpp(1)%y_centroid": 0.5 * 0.5,
+            "patch_icpp(1)%length_x": 0.5,
+            "patch_icpp(1)%length_y": 0.5,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 120.0 / (7.0 / 5.0),
+            "patch_icpp(1)%alpha_rho(1)": 120.0 * (1 - eps),
+            "patch_icpp(1)%alpha(1)": 1.0 * (1 - eps),
+            "patch_icpp(1)%alpha_rho(2)": 120.0 * eps,
+            "patch_icpp(1)%alpha(2)": 1.0 * eps,
+            # Right
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%x_centroid": 0.5 + 0.5 * 0.5,
+            "patch_icpp(2)%y_centroid": 0.5 * 0.5,
+            "patch_icpp(2)%length_x": 0.5,
+            "patch_icpp(2)%length_y": 0.5,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": 1.2 / (7.0 / 5.0),
+            "patch_icpp(2)%alpha_rho(1)": 1.2 * eps,
+            "patch_icpp(2)%alpha(1)": 1.0 * eps,
+            "patch_icpp(2)%alpha_rho(2)": 1.2 * (1 - eps),
+            "patch_icpp(2)%alpha(2)": 1.0 * (1 - eps),
+            # Fluids Physical Parameters
+            "viscous": "T",
+            "fluid_pp(1)%gamma": 1.0e00 / ((7.0 / 5.0) - 1.0e00),
+            "fluid_pp(1)%Re(1)": 1 / mu,
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(2)%gamma": 1.0e00 / ((7.0 / 5.0) - 1.0e00),
+            "fluid_pp(2)%Re(1)": 1 / mu,
+            "fluid_pp(2)%pi_inf": 0.0,
+        }
+    )
+)

src/common/include/2dHardcodedIC.fpp

@@ -4,6 +4,10 @@
     real(wp) :: r, rmax, gam, umax, p0
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph
     real(wp) :: factor
+    ! # 207
+    real(wp) :: sigma, gauss1, gauss2
+    ! # 208
+    real(wp) :: ei, d, fsm, alpha_air, alpha_sf6
 
     eps = 1.e-9_wp
 #:enddef
@@ -129,6 +133,29 @@
             q_prim_vf(advxe)%sf(i, j, 0) = patch_icpp(1)%alpha(2)
         end if
 
+    case (207) ! Kelvin Helmholtz Instability
+        sigma = 0.05_wp / sqrt(2.0_wp)
+        gauss1 = exp(- (y_cc(j) - 0.75_wp) ** 2 / (2.0_wp * sigma ** 2))
+        gauss2 = exp(- (y_cc(j) - 0.25_wp) ** 2 / (2.0_wp * sigma ** 2))
+        q_prim_vf(momxb + 1)%sf(i, j, 0) = &
+                0.1_wp * sin(4.0_wp * pi * x_cc(i)) * (gauss1 + gauss2)
+
+    case (208) ! Richtmeyer Meshkov Instability
+        lam = 1.0_wp
+        eps = 1.0e-6_wp
+        ei = 5.0_wp
+        ! Smoothening function to smooth out sharp discontinuity in the interface
+        if(x_cc(i) <= 0.7_wp * lam) then
+                d = x_cc(i) - lam * (0.4_wp - 0.1_wp * sin(2.0_wp * pi * (y_cc(j) / lam + 0.25_wp)))
+                fsm = 0.5_wp * (1.0_wp + erf(d / (ei * sqrt(dx * dy))))
+                alpha_air = eps + (1.0_wp - 2.0_wp * eps) * fsm
+                alpha_sf6 = 1.0_wp - alpha_air
+                q_prim_vf(contxb)%sf(i, j, 0) = alpha_sf6 * 5.04_wp
+                q_prim_vf(contxe)%sf(i, j, 0) = alpha_air * 1.0_wp
+                q_prim_vf(advxb)%sf(i, j, 0) = alpha_sf6
+                q_prim_vf(advxe)%sf(i, j, 0) = alpha_air 
+        end if
+
     case (250) ! MHD Orszag-Tang vortex
         ! gamma = 5/3
         !   rho = 25/(36*pi)