Skip to content

Three new examples cases #1060

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Nov 27, 2025
Merged

Three new examples cases #1060

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
831cbdc
Add 3 new examples. (2D_kelvin_helmholtz, 2D_richtmyer_meshkov, 2D_vi…
Nov 25, 2025
22e3756
Remove initial conditions in README's from previous commit
Nov 25, 2025
6f6e3d3
Fix formatting issues.
Nov 25, 2025
d66c7ef
Fix minor precision mismatch
Nov 25, 2025
1d9d0ec
Merge branch 'master' into threeNewExamples
anandrdbz Nov 25, 2025
53b3ae9
Run ./mfc.sh format
Nov 25, 2025
384d0b1
Use WENO3 instead of WENO5
Nov 27, 2025
ee67776
Add golden files
Nov 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions examples/2D_kelvin_helmholtz/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"/>
89 changes: 89 additions & 0 deletions examples/2D_kelvin_helmholtz/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
}
)
)
Binary file added examples/2D_kelvin_helmholtz/figure0.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/2D_kelvin_helmholtz/figure1.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
10 changes: 10 additions & 0 deletions examples/2D_richtmyer_meshkov/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"/>
98 changes: 98 additions & 0 deletions examples/2D_richtmyer_meshkov/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,

}
)
)
Binary file added examples/2D_richtmyer_meshkov/figure0.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/2D_richtmyer_meshkov/figure1.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
10 changes: 10 additions & 0 deletions examples/2D_viscous_shock_tube/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"/>
91 changes: 91 additions & 0 deletions examples/2D_viscous_shock_tube/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
}
)
)
Binary file added examples/2D_viscous_shock_tube/figure0.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/2D_viscous_shock_tube/figure1.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
27 changes: 27 additions & 0 deletions src/common/include/2dHardcodedIC.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 * sigma ** 2))
gauss2 = exp(- (y_cc(j) - 0.25_wp) ** 2 / (2 * 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 * 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 - alpha_air
q_prim_vf(contxb)%sf(i, j, 0) = alpha_sf6 * 5.04
q_prim_vf(contxe)%sf(i, j, 0) = alpha_air * 1.0
q_prim_vf(advxb)%sf(i, j, 0) = alpha_sf6
q_prim_vf(advxe)%sf(i, j, 0) = alpha_air
end if
Comment on lines 143 to 157
Copy link
Contributor

@coderabbitai coderabbitai bot Nov 25, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Case 208 overwrites the shared eps variable and uses inconsistent precision literals.

  1. Line 145 reassigns eps = 1.0e-6_wp, which shadows the module-level eps initialized at Line 12. If another case runs after case 208 in the same execution context, the modified eps value could cause unexpected behavior.

  2. Several literals lack the _wp suffix for working precision: 1.0 on lines 152, 154; 5.04 on line 153. Per coding guidelines, use wp for consistent precision.

Apply this diff to fix precision consistency:

     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 * lam) then
+        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 - alpha_air
-                q_prim_vf(contxb)%sf(i, j, 0) = alpha_sf6 * 5.04
-                q_prim_vf(contxe)%sf(i, j, 0) = alpha_air * 1.0
+                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

Note: If you need a local epsilon distinct from the module-level eps, consider using a differently named variable (e.g., eps_local) to avoid shadowing.

Committable suggestion skipped: line range outside the PR's diff.

🤖 Prompt for AI Agents 
In src/common/include/2dHardcodedIC.fpp around lines 143-157, avoid shadowing
the module-level eps by introducing a locally named variable (e.g., eps_local)
instead of reassigning eps; set eps_local = 1.0e-6_wp and replace uses of eps in
this case with eps_local. Also make all numeric literals use working precision
suffix _wp (change 1.0 -> 1.0_wp and 5.04 -> 5.04_wp and any other plain
literals in these lines), and ensure any derived computations (erf argument,
alpha assignments, q_prim_vf values) use the new eps_local and _wp-suffixed
literals for consistent precision.

Copy link
Contributor

@wilfonba wilfonba Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Cowsreal The comments about missing _wp are worth fixing here. The others aren't important.

Copy link
Contributor

@coderabbitai coderabbitai bot Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like the humans are having a chat. I'll hop back into my burrow for now. If you need me again, just tag @coderabbitai in a new comment, and I'll come hopping out!

Copy link
Contributor

@qodo-code-review qodo-code-review bot Nov 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code Implementation 🛠️

Implementation: In case (208), remove the reassignment of eps and ensure all numeric literals use the working precision suffix. Replace plain literals with _wp variants and adjust the comparison constant 0.7 to 0.7_wp. Keep using the existing eps variable without modifying its value.

Suggested change
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 * 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 - alpha_air
q_prim_vf(contxb)%sf(i, j, 0) = alpha_sf6 * 5.04
q_prim_vf(contxe)%sf(i, j, 0) = alpha_air * 1.0
q_prim_vf(advxb)%sf(i, j, 0) = alpha_sf6
q_prim_vf(advxe)%sf(i, j, 0) = alpha_air
end if
case (208) ! Richtmeyer Meshkov Instability
lam = 1.0_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

See review comment here


case (250) ! MHD Orszag-Tang vortex
! gamma = 5/3
! rho = 25/(36*pi)
Expand Down