diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 7d55e301b3..bae50068ad 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -225,6 +225,7 @@ end if Some patch configurations are not adequately handled with the above analytic variable definitions. In this case, a hard coded patch can be used. Hard coded patches can be added by adding additional hard coded patch identifiers to `src/pre_process/include/1[2,3]dHardcodedIC.fpp`. +When using a hard coded patch, the `patch_icpp(patch_id)%%hcid` must be set to the hard-coded patch id. For example, to add a 2D Hardcoded patch with an id of 200, one would add the following to `src/pre_process/include/2dHardcodedIC.fpp` ```f90 @@ -232,7 +233,7 @@ For example, to add a 2D Hardcoded patch with an id of 200, one would add the fo ! Primitive variables assignment ``` -and use `patch_icpp(i)%%geometry = 7` and `patch_icpp(i)%%hcid = 200` in the input file. +and use `patch_icpp(i)%%hcid = 200` in the input file. Additional variables can be declared in `Hardcoded1[2,3]DVariables` and used in `hardcoded1[2,3]D`. As a convention, any hard coded patches that are part of the MFC master branch should be identified as 1[2,3]xx where the first digit indicates the number of dimensions. @@ -313,7 +314,7 @@ These parameters should be prepended with `patch_ib(j)%` where $j$ is the patch #### Parameter Descriptions - `geometry` defines the type of geometry of a patch with an integer number. -Definitions for currently implemented patch types are list in table [Immersed Boundary Patch Type](#immersed-boundary-patch-types) +Definitions for currently implemented patch types are list in table [Patch Types](#patch-types). - `x[y,z]_centroid` is the centroid location of the patch in the x[y,z]-direction @@ -996,16 +997,16 @@ This boundary condition can be used for subsonic inflow (`bc_[x,y,z]%[beg,end]` | 3 | Rectangle | 2 | N | Coordinate-aligned. Requires `[x,y]_centroid` and `length_[x,y]`. | | 4 | Sweep line | 2 | Y | Not coordinate aligned. Requires `[x,y]_centroid` and `normal(i)`. | | 5 | Ellipse | 2 | Y | Requires `[x,y]_centroid` and `radii(i)`. | -| 6 | N/A | 2 | N | No longer exists. Empty. | -| 7 | 2D Hardcoded | 2 | N | Assigns the primitive variables as analytical functions. | +| 6 | N/A | N/A | N/A | No longer exists. Empty. | +| 7 | N/A | N/A | N/A | No longer exists. Empty. | | 8 | Sphere | 3 | Y | Requires `[x,y,z]_centroid` and `radius` | | 9 | Cuboid | 3 | N | Coordinate-aligned. Requires `[x,y,z]_centroid` and `length_[x,y,z]`. | | 10 | Cylinder | 3 | Y | Requires `[x,y,z]_centroid`, `radius`, and `length_[x,y,z]`. | | 11 | Sweep plane | 3 | Y | Not coordinate-aligned. Requires `x[y,z]_centroid` and `normal(i)`. | | 12 | Ellipsoid | 3 | Y | Requires `[x,y,z]_centroid` and `radii(i)`. | -| 13 | 3D Hardcoded | 3 | N | Assigns the primitive variables as analytical functions | +| 13 | N/A | N/A | N/A | No longer exists. Empty. | | 14 | Spherical Harmonic | 3 | N | Requires `[x,y,z]_centroid`, `radius`, `epsilon`, `beta` | -| 15 | 1D Hardcoded | 1 | N | Assigns the primitive variables as analytical functions | +| 15 | N/A | N/A | N/A | No longer exists. Empty. | | 16 | 1D bubble pulse | 1 | N | Requires `x_centroid`, `length_x` | | 17 | Spiral | 2 | N | Requires `[x,y]_centroid` | | 18 | 2D Varcircle | 2 | Y | Requires `[x,y]_centroid`, `radius`, and `thickness` | diff --git a/examples/1D_shuosher_analytical/case.py b/examples/1D_shuosher_analytical/case.py new file mode 100644 index 0000000000..30e50acff5 --- /dev/null +++ b/examples/1D_shuosher_analytical/case.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 1000 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend, Nt = 1.8, 20000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -5.0, + "x_domain%end": 5.0, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(math.ceil(Nt / 10.0)), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Background to cover whole domain with basic line patch + # Patch 1 Left (-5 < x < -4) + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": -4.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%vel(1)": 2.629369, + "patch_icpp(1)%pres": 10.3333, + "patch_icpp(1)%alpha_rho(1)": 3.957143, + "patch_icpp(1)%alpha(1)": 1.0, + # One anlytic patch to take care of -4 < x < 5 + # Patch 2 Analytic + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.5, + "patch_icpp(2)%length_x": 9.0, + "patch_icpp(2)%vel(1)": 0.0, + "patch_icpp(2)%pres": 1.0, + "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/examples/1D_shuosher_old/case.py b/examples/1D_shuosher_old/case.py index 30e50acff5..24dd240ddf 100644 --- a/examples/1D_shuosher_old/case.py +++ b/examples/1D_shuosher_old/case.py @@ -64,7 +64,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0.0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_shuosher_teno5/case.py b/examples/1D_shuosher_teno5/case.py index 1bbe603ebd..73189b3cec 100644 --- a/examples/1D_shuosher_teno5/case.py +++ b/examples/1D_shuosher_teno5/case.py @@ -66,7 +66,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_shuosher_teno7/case.py b/examples/1D_shuosher_teno7/case.py index c0ebc43436..8a21aa1309 100644 --- a/examples/1D_shuosher_teno7/case.py +++ b/examples/1D_shuosher_teno7/case.py @@ -66,7 +66,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0.0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_shuosher_wenojs5/case.py b/examples/1D_shuosher_wenojs5/case.py index 01f18f1ec0..b99e5ef7e2 100644 --- a/examples/1D_shuosher_wenojs5/case.py +++ b/examples/1D_shuosher_wenojs5/case.py @@ -65,7 +65,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_shuosher_wenom5/case.py b/examples/1D_shuosher_wenom5/case.py index 7549d7a337..8c96fea5b6 100644 --- a/examples/1D_shuosher_wenom5/case.py +++ b/examples/1D_shuosher_wenom5/case.py @@ -65,7 +65,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_shuosher_wenoz5/case.py b/examples/1D_shuosher_wenoz5/case.py index 7750fb4b1c..a7e8e3e3cd 100644 --- a/examples/1D_shuosher_wenoz5/case.py +++ b/examples/1D_shuosher_wenoz5/case.py @@ -65,7 +65,8 @@ "patch_icpp(2)%length_x": 9.0, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)", + "patch_icpp(2)%alpha_rho(1)": 0, + "patch_icpp(2)%hcid": 180, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_titarevtorro/case.py b/examples/1D_titarevtorro/case.py index 3a51bd86b0..63dfe39b79 100644 --- a/examples/1D_titarevtorro/case.py +++ b/examples/1D_titarevtorro/case.py @@ -62,7 +62,8 @@ "patch_icpp(2)%length_x": 9.5, "patch_icpp(2)%vel(1)": 0.0, "patch_icpp(2)%pres": 1.0, - "patch_icpp(2)%alpha_rho(1)": "1 + 0.1*sin(20*x*pi)", + "patch_icpp(2)%alpha_rho(1)": 0.0, + "patch_icpp(2)%hcid": 181, "patch_icpp(2)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/1D_titarevtorro_analytical/case.py b/examples/1D_titarevtorro_analytical/case.py new file mode 100644 index 0000000000..3a51bd86b0 --- /dev/null +++ b/examples/1D_titarevtorro_analytical/case.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 999 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend, Nt = 5, 20000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -5.0, + "x_domain%end": 5.0, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(math.ceil(Nt / 10.0)), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L (-5 < x < -4.5) + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": -4.75, + "patch_icpp(1)%length_x": 0.5, + "patch_icpp(1)%vel(1)": 0.523326, + "patch_icpp(1)%pres": 1.805, + "patch_icpp(1)%alpha_rho(1)": 1.515695, + "patch_icpp(1)%alpha(1)": 1.0, + # Patch 2 R (-4.5 < x < 5) + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.25, + "patch_icpp(2)%length_x": 9.5, + "patch_icpp(2)%vel(1)": 0.0, + "patch_icpp(2)%pres": 1.0, + "patch_icpp(2)%alpha_rho(1)": "1 + 0.1*sin(20*x*pi)", + "patch_icpp(2)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/examples/2D_GreshoVortex/case.py b/examples/2D_GreshoVortex/case.py index dad9435599..96f3638333 100644 --- a/examples/2D_GreshoVortex/case.py +++ b/examples/2D_GreshoVortex/case.py @@ -67,7 +67,7 @@ "fd_order": 1, "omega_wrt(3)": "T", # Patch 1 - "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%geometry": 3, "patch_icpp(1)%hcid": 203, "patch_icpp(1)%x_centroid": 0.5, "patch_icpp(1)%y_centroid": 0.5, diff --git a/examples/2D_acoustic_pulse/case.py b/examples/2D_acoustic_pulse/case.py index 73cec89a66..899242ca48 100644 --- a/examples/2D_acoustic_pulse/case.py +++ b/examples/2D_acoustic_pulse/case.py @@ -81,8 +81,9 @@ "patch_icpp(2)%radius": 1.0, "patch_icpp(2)%vel(1)": 0, "patch_icpp(2)%vel(2)": 0, - "patch_icpp(2)%pres": f"{p_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**({gam} / ({gam} - 1))", - "patch_icpp(2)%alpha_rho(1)": f"{rho_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**(1 / ({gam} - 1))", + "patch_icpp(2)%pres": 0.0, + "patch_icpp(2)%alpha_rho(1)": 0.0, + "patch_icpp(2)%hcid": 281, "patch_icpp(2)%alpha(1)": 1.0, "patch_icpp(2)%alter_patch(1)": "T", # CBC Inflow / Outflow diff --git a/examples/2D_acoustic_pulse_analytical/case.py b/examples/2D_acoustic_pulse_analytical/case.py new file mode 100644 index 0000000000..73cec89a66 --- /dev/null +++ b/examples/2D_acoustic_pulse_analytical/case.py @@ -0,0 +1,102 @@ +import math +import json + +# Numerical setup +Nx = 99 +Ny = 99 +dx = 8.0 / (1.0 * (Nx + 1)) + +alf_st = 0.4 + +p_inf = 101325 +rho_inf = 1 +gam = 1.4 + +c = math.sqrt(gam * (p_inf) / rho_inf) +cfl = 0.3 +mydt = cfl * dx / c +Tfinal = 80 * 1 / c +Nt = int(Tfinal / mydt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -4, + "x_domain%end": 4, + "y_domain%beg": -4, + "y_domain%end": 4, + "m": Nx, + "n": Ny, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": int(Nt / 100), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -8, + "bc_x%end": -8, + "bc_y%beg": -8, + "bc_y%end": -8, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "T", + "omega_wrt(3)": "T", + "fd_order": 2, + # Patch 1 + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0, + "patch_icpp(1)%y_centroid": 0, + "patch_icpp(1)%length_x": 8.0, + "patch_icpp(1)%length_y": 8.0, + "patch_icpp(1)%vel(1)": 0, + "patch_icpp(1)%vel(2)": 0, + "patch_icpp(1)%pres": p_inf, + "patch_icpp(1)%alpha_rho(1)": rho_inf, + "patch_icpp(1)%alpha(1)": 1.0, + # Patch 2 + "patch_icpp(2)%geometry": 2, + "patch_icpp(2)%x_centroid": 0, + "patch_icpp(2)%y_centroid": 0, + "patch_icpp(2)%radius": 1.0, + "patch_icpp(2)%vel(1)": 0, + "patch_icpp(2)%vel(2)": 0, + "patch_icpp(2)%pres": f"{p_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**({gam} / ({gam} - 1))", + "patch_icpp(2)%alpha_rho(1)": f"{rho_inf}*(1 - 0.5*({gam} - 1)*({alf_st})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**(1 / ({gam} - 1))", + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%alter_patch(1)": "T", + # CBC Inflow / Outflow + "bc_x%grcbc_in": "F", + "bc_x%grcbc_out": "T", + "bc_x%grcbc_vel_out": "F", + "bc_x%pres_out": p_inf, + "bc_y%grcbc_in": "F", + "bc_y%grcbc_out": "T", + "bc_y%grcbc_vel_out": "F", + "bc_y%pres_out": p_inf, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (gam - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/examples/2D_hardcodied_ic/case.py b/examples/2D_hardcodied_ic/case.py index 52132bf5ef..e26e01f8ac 100644 --- a/examples/2D_hardcodied_ic/case.py +++ b/examples/2D_hardcodied_ic/case.py @@ -55,7 +55,7 @@ "prim_vars_wrt": "T", "parallel_io": "T", # Patch 1: Base - "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%geometry": 3, "patch_icpp(1)%hcid": 200, "patch_icpp(1)%x_centroid": 4.0, "patch_icpp(1)%y_centroid": 4.0, diff --git a/examples/2D_isentropicvortex/case.py b/examples/2D_isentropicvortex/case.py index 96acf97eee..a3b2b63e78 100644 --- a/examples/2D_isentropicvortex/case.py +++ b/examples/2D_isentropicvortex/case.py @@ -101,10 +101,11 @@ "patch_icpp(1)%y_centroid": 0, "patch_icpp(1)%length_x": 10.0, "patch_icpp(1)%length_y": 10.0, - "patch_icpp(1)%vel(1)": vel1, - "patch_icpp(1)%vel(2)": vel2, - "patch_icpp(1)%pres": pres, - "patch_icpp(1)%alpha_rho(1)": alpha_rho1, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": 0.0, + "patch_icpp(1)%alpha_rho(1)": 0.0, + "patch_icpp(1)%hcid": 280, "patch_icpp(1)%alpha(1)": 1.0, # Fluids Physical Parameters "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), diff --git a/examples/2D_isentropicvortex_analytical/case.py b/examples/2D_isentropicvortex_analytical/case.py new file mode 100644 index 0000000000..f25773e1a5 --- /dev/null +++ b/examples/2D_isentropicvortex_analytical/case.py @@ -0,0 +1,115 @@ +import math +import json + +# Parameters +epsilon = "5d0" +alpha = "1d0" +gamma = "1.4" + +# Initial conditions +vel1_i = "0d0" +vel2_i = "0d0" +T_i = "1d0" +alpha_rho1_i = "1d0" +pres_i = "1d0" + +# Perturbations +vel1 = f"{vel1_i} + (y - yc)*({epsilon}/(2d0*pi))*" + f"exp({alpha}*(1d0 - (x - xc)**2d0 - (y - yc)**2d0))" +vel2 = f"{vel2_i} - (x - xc)*({epsilon}/(2d0*pi))*" + f"exp({alpha}*(1d0 - (x - xc)**2d0 - (y - yc)**2d0))" +alpha_rho1 = ( + f"{alpha_rho1_i}*(1d0 - ({alpha_rho1_i}/{pres_i})*({epsilon}/(2d0*pi))*" + + f"({epsilon}/(8d0*{alpha}*({gamma} + 1d0)*pi))*" + + f"exp(2d0*{alpha}*(1d0 - (x - xc)**2d0" + + f"- (y - yc)**2d0))" + + f")**{gamma}" +) +pres = ( + f"{pres_i}*(1d0 - ({alpha_rho1_i}/{pres_i})*({epsilon}/(2d0*pi))*" + + f"({epsilon}/(8d0*{alpha}*({gamma} + 1d0)*pi))*" + + f"exp(2d0*{alpha}*(1d0 - (x - xc)**2d0" + + f"- (y - yc)**2d0))" + + f")**({gamma} + 1d0)" +) + +# Numerical setup +Nx = 399 +dx = 10.0 / (1.0 * (Nx + 1)) + +c = 1.4**0.5 +C = 0.3 +mydt = C * dx / c * 0.01 +Nt = 1 / mydt + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -3, + "x_domain%end": 3, + "y_domain%beg": -3, + "y_domain%end": 3, + "stretch_x": "T", + "stretch_y": "T", + "loops_x": 2, + "loops_y": 2, + "a_x": 1.03, + "a_y": 1.03, + "x_a": -1.5, + "y_a": -1.5, + "x_b": 1.5, + "y_b": 1.5, + "m": Nx, + "n": Nx, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": 10, + "t_step_save": 1, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 3, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -4, + "bc_x%end": -4, + "bc_y%beg": -4, + "bc_y%end": -4, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "omega_wrt(3)": "T", + "fd_order": 2, + # Patch 1 + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0, + "patch_icpp(1)%y_centroid": 0, + "patch_icpp(1)%length_x": 10.0, + "patch_icpp(1)%length_y": 10.0, + "patch_icpp(1)%vel(1)": vel1, + "patch_icpp(1)%vel(2)": vel2, + "patch_icpp(1)%pres": pres, + "patch_icpp(1)%alpha_rho(1)": alpha_rho1, + "patch_icpp(1)%hcid": 280, + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/examples/2D_lungwave/case.py b/examples/2D_lungwave/case.py index 5a60492463..85b0387919 100644 --- a/examples/2D_lungwave/case.py +++ b/examples/2D_lungwave/case.py @@ -153,7 +153,7 @@ "patch_icpp(1)%alpha(1)": alphal_back, "patch_icpp(1)%alpha(2)": alphag_back, # Patch 2: Lung - "patch_icpp(2)%geometry": 7, + "patch_icpp(2)%geometry": 3, "patch_icpp(2)%hcid": 205, "patch_icpp(2)%alter_patch(1)": "T", "patch_icpp(2)%x_centroid": dlengx / 2.0, diff --git a/examples/2D_lungwave_horizontal/case.py b/examples/2D_lungwave_horizontal/case.py index 235de7154f..8fea0a0285 100644 --- a/examples/2D_lungwave_horizontal/case.py +++ b/examples/2D_lungwave_horizontal/case.py @@ -149,7 +149,7 @@ "patch_icpp(1)%alpha(1)": alphal_back, "patch_icpp(1)%alpha(2)": alphag_back, # Patch 2: Lung - "patch_icpp(2)%geometry": 7, + "patch_icpp(2)%geometry": 3, "patch_icpp(2)%hcid": 206, "patch_icpp(2)%alter_patch(1)": "T", "patch_icpp(2)%x_centroid": -dlengx / 4.0, diff --git a/examples/2D_orszag_tang/case.py b/examples/2D_orszag_tang/case.py index 5f19ced93d..ae2680144c 100644 --- a/examples/2D_orszag_tang/case.py +++ b/examples/2D_orszag_tang/case.py @@ -54,7 +54,7 @@ # v = (-sin(2*pi*y), sin(2*pi*x), 0) # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) "patch_icpp(1)%hcid": 250, - "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%geometry": 3, "patch_icpp(1)%x_centroid": 0.5, "patch_icpp(1)%y_centroid": 0.5, "patch_icpp(1)%length_x": 1.0, diff --git a/examples/2D_orszag_tang_powell/case.py b/examples/2D_orszag_tang_powell/case.py index 174a80d191..6c8aa09069 100644 --- a/examples/2D_orszag_tang_powell/case.py +++ b/examples/2D_orszag_tang_powell/case.py @@ -56,7 +56,7 @@ # v = (-sin(2*pi*y), sin(2*pi*x), 0) # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) "patch_icpp(1)%hcid": 250, - "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%geometry": 3, "patch_icpp(1)%x_centroid": 0.5, "patch_icpp(1)%y_centroid": 0.5, "patch_icpp(1)%length_x": 1.0, diff --git a/examples/2D_rayleigh_taylor/case.py b/examples/2D_rayleigh_taylor/case.py index 635af8b6e2..029178732a 100644 --- a/examples/2D_rayleigh_taylor/case.py +++ b/examples/2D_rayleigh_taylor/case.py @@ -85,7 +85,7 @@ "p_y": 0.0, "g_y": -9.81, # Water Patch - "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%geometry": 3, "patch_icpp(1)%hcid": 204, "patch_icpp(1)%x_centroid": 0, "patch_icpp(1)%y_centroid": h / 2, diff --git a/examples/2D_zero_circ_vortex/case.py b/examples/2D_zero_circ_vortex/case.py index 793f265b7c..2eb8b24643 100644 --- a/examples/2D_zero_circ_vortex/case.py +++ b/examples/2D_zero_circ_vortex/case.py @@ -81,10 +81,9 @@ "patch_icpp(2)%x_centroid": 0, "patch_icpp(2)%y_centroid": 0, "patch_icpp(2)%radius": 1.0, - "patch_icpp(2)%vel(1)": f"{u_inf}*(1 - ({Mv} / {M_inf}))*y*exp(0.5*(1 - sqrt(x**2 + y**2))) ", - "patch_icpp(2)%vel(2)": f"{u_inf}*(({Mv} / {M_inf}))*x*exp(0.5*(1 - sqrt(x**2 + y**2)))", - "patch_icpp(2)%pres": f"{p_inf}*(1 - 0.5*({gam} - 1)*({Mv} / {M_inf})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**({gam} / ({gam} - 1))", - "patch_icpp(2)%alpha_rho(1)": f"{rho_inf}*(1 - 0.5*({gam} - 1)*({Mv} / {M_inf})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**(1 / ({gam} - 1))", + "patch_icpp(2)%vel(1)": 0.0, + "patch_icpp(2)%vel(2)": 0.0, + "patch_icpp(2)%hcid": 282, "patch_icpp(2)%alpha(1)": 1.0, "patch_icpp(2)%alter_patch(1)": "T", # CBC Inflow / Outflow diff --git a/examples/2D_zero_circ_vortex_analytical/case.py b/examples/2D_zero_circ_vortex_analytical/case.py new file mode 100644 index 0000000000..793f265b7c --- /dev/null +++ b/examples/2D_zero_circ_vortex_analytical/case.py @@ -0,0 +1,109 @@ +import math +import json + +# Numerical setup +Nx = 250 +Ny = 125 +dx = 40.0 / (1.0 * (Nx + 1)) + +M_inf = 0.3 +Mv = 0.1 + +p_inf = 101325 +rho_inf = 1 +gam = 1.4 + +c = math.sqrt(gam * (p_inf) / rho_inf) +u_inf = M_inf * c +cfl = 0.3 +mydt = cfl * dx / c +Tfinal = 150 * 1 / u_inf +Nt = int(Tfinal / mydt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -20, + "x_domain%end": 20, + "y_domain%beg": -10, + "y_domain%end": 10, + "m": Nx, + "n": Ny, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": int(Nt / 100), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -7, + "bc_x%end": -8, + "bc_y%beg": -6, + "bc_y%end": -6, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "T", + "omega_wrt(3)": "T", + "fd_order": 2, + # Patch 1 + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0, + "patch_icpp(1)%y_centroid": 0, + "patch_icpp(1)%length_x": 40.0, + "patch_icpp(1)%length_y": 20.0, + "patch_icpp(1)%vel(1)": u_inf, + "patch_icpp(1)%vel(2)": 0, + "patch_icpp(1)%pres": 101325, + "patch_icpp(1)%alpha_rho(1)": 1, + "patch_icpp(1)%alpha(1)": 1.0, + # Patch 2 + "patch_icpp(2)%geometry": 2, + "patch_icpp(2)%x_centroid": 0, + "patch_icpp(2)%y_centroid": 0, + "patch_icpp(2)%radius": 1.0, + "patch_icpp(2)%vel(1)": f"{u_inf}*(1 - ({Mv} / {M_inf}))*y*exp(0.5*(1 - sqrt(x**2 + y**2))) ", + "patch_icpp(2)%vel(2)": f"{u_inf}*(({Mv} / {M_inf}))*x*exp(0.5*(1 - sqrt(x**2 + y**2)))", + "patch_icpp(2)%pres": f"{p_inf}*(1 - 0.5*({gam} - 1)*({Mv} / {M_inf})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**({gam} / ({gam} - 1))", + "patch_icpp(2)%alpha_rho(1)": f"{rho_inf}*(1 - 0.5*({gam} - 1)*({Mv} / {M_inf})**2*exp(0.5*(1 - sqrt(x**2 + y**2))))**(1 / ({gam} - 1))", + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%alter_patch(1)": "T", + # CBC Inflow / Outflow + "bc_x%grcbc_in": "T", + "bc_x%grcbc_out": "T", + "bc_x%grcbc_vel_out": "T", + "bc_x%vel_in(1)": u_inf, + "bc_x%vel_in(2)": 0, + "bc_x%vel_in(3)": 0, + "bc_x%pres_in": p_inf, + "bc_x%alpha_rho_in(1)": rho_inf, + "bc_x%alpha_in(1)": 1, + "bc_x%vel_out(1)": u_inf, + "bc_x%vel_out(2)": 0, + "bc_x%vel_out(3)": 0, + "bc_x%pres_out": p_inf, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (gam - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/examples/3D_TaylorGreenVortex/case.py b/examples/3D_TaylorGreenVortex/case.py index de440980cb..8b48c2990f 100644 --- a/examples/3D_TaylorGreenVortex/case.py +++ b/examples/3D_TaylorGreenVortex/case.py @@ -85,10 +85,11 @@ "patch_icpp(1)%length_x": 2 * math.pi * L, "patch_icpp(1)%length_y": 2 * math.pi * L, "patch_icpp(1)%length_z": 2 * math.pi * L, - "patch_icpp(1)%vel(1)": f"{V0}*sin(x/{L})*cos(y/{L})*sin(z/{L})", - "patch_icpp(1)%vel(2)": f"-{V0}*cos(x/{L})*sin(y/{L})*sin(z/{L})", + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, "patch_icpp(1)%vel(3)": 0, - "patch_icpp(1)%pres": f"{P0} + ({rho0}*{V0}**2/16)*(cos(2*x/{L}) + cos(2*y/{L}))*(cos(2*z/{L}) + 2)", + "patch_icpp(1)%pres": 0.0, + "patch_icpp(1)%hcid": 380, "patch_icpp(1)%alpha_rho(1)": 1, "patch_icpp(1)%alpha(1)": 1, # Fluids Physical Parameters diff --git a/examples/3D_TaylorGreenVortex_analytical/case.py b/examples/3D_TaylorGreenVortex_analytical/case.py new file mode 100644 index 0000000000..de440980cb --- /dev/null +++ b/examples/3D_TaylorGreenVortex_analytical/case.py @@ -0,0 +1,100 @@ +#!/usr/bin/env python3 +import math +import json + +N = 256 + +Re = 1600 +L = 1 +P0 = 101325 +rho0 = 1 +C0 = math.sqrt(1.4 * P0) +V0 = 0.1 * C0 +mu = V0 * L / Re + +cfl = 0.5 +dx = 2 * math.pi * L / (N + 1) + +dt = cfl * dx / (C0) + +tC = L / V0 +tEnd = 20 * tC + +Nt = int(tEnd / dt) + + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -math.pi * L, + "x_domain%end": math.pi * L, + "y_domain%beg": -math.pi * L, + "y_domain%end": math.pi * L, + "z_domain%beg": -math.pi * L, + "z_domain%end": math.pi * L, + "m": N, + "n": N, + "p": N, + "cyl_coord": "F", + "dt": dt, + "t_step_start": 13529, + "t_step_stop": Nt, + "t_step_save": int(Nt / 100), + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mixture_err": "T", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + "bc_z%beg": -1, + "bc_z%end": -1, + "viscous": "T", + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + # 'prim_vars_wrt' :'T', + "omega_wrt(1)": "T", + "omega_wrt(2)": "T", + "omega_wrt(3)": "T", + "qm_wrt": "T", + "fd_order": 4, + "parallel_io": "T", + # I will use 1 for WATER properties, and 2 for AIR properties + # Patch 1: Background (AIR - 2) + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0, + "patch_icpp(1)%y_centroid": 0, + "patch_icpp(1)%z_centroid": 0, + "patch_icpp(1)%length_x": 2 * math.pi * L, + "patch_icpp(1)%length_y": 2 * math.pi * L, + "patch_icpp(1)%length_z": 2 * math.pi * L, + "patch_icpp(1)%vel(1)": f"{V0}*sin(x/{L})*cos(y/{L})*sin(z/{L})", + "patch_icpp(1)%vel(2)": f"-{V0}*cos(x/{L})*sin(y/{L})*sin(z/{L})", + "patch_icpp(1)%vel(3)": 0, + "patch_icpp(1)%pres": f"{P0} + ({rho0}*{V0}**2/16)*(cos(2*x/{L}) + cos(2*y/{L}))*(cos(2*z/{L}) + 2)", + "patch_icpp(1)%alpha_rho(1)": 1, + "patch_icpp(1)%alpha(1)": 1, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1), + "fluid_pp(1)%pi_inf": 0, + "fluid_pp(1)%Re(1)": 1 / mu, + } + ) +) diff --git a/examples/3D_rayleigh_taylor/case.py b/examples/3D_rayleigh_taylor/case.py index 35b0ecef24..73a9c9a88b 100644 --- a/examples/3D_rayleigh_taylor/case.py +++ b/examples/3D_rayleigh_taylor/case.py @@ -92,7 +92,7 @@ "p_y": 0.0, "g_y": -98.1, # Water Patch - "patch_icpp(1)%geometry": 13, + "patch_icpp(1)%geometry": 9, "patch_icpp(1)%hcid": 300, "patch_icpp(1)%x_centroid": 0, "patch_icpp(1)%y_centroid": h / 2, diff --git a/src/pre_process/include/1dHardcodedIC.fpp b/src/pre_process/include/1dHardcodedIC.fpp index 390409bfc0..e7896fd9f6 100644 --- a/src/pre_process/include/1dHardcodedIC.fpp +++ b/src/pre_process/include/1dHardcodedIC.fpp @@ -3,12 +3,23 @@ #:enddef #:def Hardcoded1D() - select case (patch_icpp(patch_id)%hcid) case (170) ! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera, SDtoolbox) @: HardcodedReadValues() + case (180) + ! This is patch is hard-coded for test suite optimization used in the + ! 1D_shuoser cases: "patch_icpp(2)%alpha_rho(1)": "1 + 0.2*sin(5*x)" + if (patch_id == 2) then + q_prim_vf(contxb + 0)%sf(i, 0, 0) = 1 + 0.2*sin(5*x_cc(i)) + end if + + case (181) + ! This is patch is hard-coded for test suite optimization used in the + ! 1D_titarevtorro cases: "patch_icpp(2)%alpha_rho(1)": "1 + 0.1*sin(20*x*pi)" + q_prim_vf(contxb + 0)%sf(i, 0, 0) = 1 + 0.1*sin(20*x_cc(i)*3.141592653589793) + case default call s_int_to_str(patch_id, iStr) call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr)) diff --git a/src/pre_process/include/2dHardcodedIC.fpp b/src/pre_process/include/2dHardcodedIC.fpp index e9d81789f4..7136b1d072 100644 --- a/src/pre_process/include/2dHardcodedIC.fpp +++ b/src/pre_process/include/2dHardcodedIC.fpp @@ -143,7 +143,6 @@ q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = sin(4._wp*pi*x_cc(i))/sqrt(4._wp*pi) case (251) ! RMHD Cylindrical Blast Wave [Mignone, 2006: Section 4.3.1] - if (x_cc(i)**2 + y_cc(j)**2 < 0.08_wp**2) then q_prim_vf(contxb)%sf(i, j, 0) = 0.01 q_prim_vf(E_idx)%sf(i, j, 0) = 1.0 @@ -161,6 +160,37 @@ ! This hardcoded case extrudes a 1D profile to initialize a 2D simulation domain @: HardcodedReadValues() + case (280) + ! This is patch is hard-coded for test suite optimization used in the + ! 2D_isentropicvortex case: + ! This analytic patch uses geometry 2 + if (patch_id == 1) then + q_prim_vf(E_idx)%sf(i, j, 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*3.141592653589793))*(5.0/(8.0*1.0*(1.4 + 1.0)*3.141592653589793))*exp(2.0*1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**(1.4 + 1.0) + q_prim_vf(contxb + 0)%sf(i, j, 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*3.141592653589793))*(5.0/(8.0*1.0*(1.4 + 1.0)*3.141592653589793))*exp(2.0*1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**1.4 + q_prim_vf(momxb + 0)%sf(i, j, 0) = 0.0 + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*3.141592653589793))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + q_prim_vf(momxb + 1)%sf(i, j, 0) = 0.0 - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*3.141592653589793))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + end if + + case (281) + ! This is patch is hard-coded for test suite optimization used in the + ! 2D_acoustic_pulse case: + ! This analytic patch uses geometry 2 + if (patch_id == 2) then + q_prim_vf(E_idx)%sf(i, j, 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1)) + q_prim_vf(contxb + 0)%sf(i, j, 0) = 1*(1 - 0.5*(1.4 - 1)*(0.4)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1)) + end if + + case (282) + ! This is patch is hard-coded for test suite optimization used in the + ! 2D_zero_circ_vortex case: + ! This analytic patch uses geometry 2 + if (patch_id == 2) then + q_prim_vf(E_idx)%sf(i, j, 0) = 101325*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1.4/(1.4 - 1)) + q_prim_vf(contxb + 0)%sf(i, j, 0) = 1*(1 - 0.5*(1.4 - 1)*(0.1/0.3)**2*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))))**(1/(1.4 - 1)) + q_prim_vf(momxb + 0)%sf(i, j, 0) = 112.99092883944267*(1 - (0.1/0.3))*y_cc(j)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))) + q_prim_vf(momxb + 1)%sf(i, j, 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))) + end if + case default if (proc_rank == 0) then call s_int_to_str(patch_id, iStr) diff --git a/src/pre_process/include/3dHardcodedIC.fpp b/src/pre_process/include/3dHardcodedIC.fpp index 804b89ab73..797272eadf 100644 --- a/src/pre_process/include/3dHardcodedIC.fpp +++ b/src/pre_process/include/3dHardcodedIC.fpp @@ -59,6 +59,16 @@ ! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain @: HardcodedReadValues() + case (380) + ! This is patch is hard-coded for test suite optimization used in the + ! 3D_TaylorGreenVortex case: + ! This analytic patch used geometry 9 + if (patch_id == 1) then + q_prim_vf(E_idx)%sf(i, j, k) = 101325 + (1*37.6636429464809**2/16)*(cos(2*x_cc(i)/1) + cos(2*y_cc(j)/1))*(cos(2*z_cc(k)/1) + 2) + q_prim_vf(momxb + 0)%sf(i, j, k) = 37.6636429464809*sin(x_cc(i)/1)*cos(y_cc(j)/1)*sin(z_cc(k)/1) + q_prim_vf(momxb + 1)%sf(i, j, k) = -37.6636429464809*cos(x_cc(i)/1)*sin(y_cc(j)/1)*sin(z_cc(k)/1) + end if + case default call s_int_to_str(patch_id, iStr) call s_mpi_abort("Invalid hcid specified for patch "//trim(iStr)) diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp index 969b6cd471..681376e0aa 100644 --- a/src/pre_process/m_check_patches.fpp +++ b/src/pre_process/m_check_patches.fpp @@ -62,7 +62,9 @@ contains 'is no longer supported for patch '//trim(iStr)// & '. Exiting.') elseif (patch_icpp(i)%geometry == 7) then - call s_check_2D_analytical_patch_geometry(i) + call s_mpi_abort('geometry 7 (formerly "2D Hardcoded")'// & + 'is no longer supported for patch '//trim(iStr)// & + '. Exiting.') elseif (patch_icpp(i)%geometry == 8) then call s_check_sphere_patch_geometry(i) elseif (patch_icpp(i)%geometry == 9) then @@ -74,11 +76,15 @@ contains elseif (patch_icpp(i)%geometry == 12) then call s_check_ellipsoid_patch_geometry(i) elseif (patch_icpp(i)%geometry == 13) then - call s_check_3D_analytical_patch_geometry(i) + call s_mpi_abort('geometry 13 (formerly "3D Hardcoded")'// & + 'is no longer supported for patch '//trim(iStr)// & + '. Exiting.') elseif (patch_icpp(i)%geometry == 14) then call s_check_spherical_harmonic_patch_geometry(i) elseif (patch_icpp(i)%geometry == 15) then - call s_check_1d_analytical_patch_geometry(i) + call s_mpi_abort('geometry 15 (formerly "1D Hardcoded")'// & + 'is no longer supported for patch '//trim(iStr)// & + '. Exiting.') elseif (patch_icpp(i)%geometry == 16) then print *, '1d pressure sinusoid' elseif (patch_icpp(i)%geometry == 17) then @@ -243,54 +249,6 @@ contains end subroutine s_check_2D_TaylorGreen_vortex_patch_geometry - !> This subroutine checks the model patch input - !! @param patch_id Patch identifier - impure subroutine s_check_1D_analytical_patch_geometry(patch_id) - - integer, intent(in) :: patch_id - call s_int_to_str(patch_id, iStr) - - @:PROHIBIT(n > 0, "1D analytical patch "//trim(iStr)//": n must be zero") - @:PROHIBIT(p > 0, "1D analytical patch "//trim(iStr)//": p must be zero") - @:PROHIBIT(model_eqns /= 4 .and. model_eqns /= 2, "1D analytical patch "//trim(iStr)//": model_eqns must be either 4 or 2") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "1D analytical patch "//trim(iStr)//": x_centroid must be set") - @:PROHIBIT(patch_icpp(patch_id)%length_x <= 0._wp, "1D analytical patch "//trim(iStr)//": length_x must be greater than zero") - - end subroutine s_check_1D_analytical_patch_geometry - - !> This subroutine checks the model patch input - !! @param patch_id Patch identifier - impure subroutine s_check_2D_analytical_patch_geometry(patch_id) - - integer, intent(in) :: patch_id - call s_int_to_str(patch_id, iStr) - - @:PROHIBIT(n == 0, "2D analytical patch "//trim(iStr)//": n must be greater than zero") - @:PROHIBIT(p > 0, "2D analytical patch "//trim(iStr)//": p must be zero") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "2D analytical patch "//trim(iStr)//": x_centroid must be set") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "2D analytical patch "//trim(iStr)//": y_centroid must be set") - @:PROHIBIT(patch_icpp(patch_id)%length_x <= 0._wp, "2D analytical patch "//trim(iStr)//": length_x must be greater than zero") - @:PROHIBIT(patch_icpp(patch_id)%length_y <= 0._wp, "2D analytical patch "//trim(iStr)//": length_y must be greater than zero") - - end subroutine s_check_2D_analytical_patch_geometry - - !> This subroutine checks the model patch input - !! @param patch_id Patch identifier - impure subroutine s_check_3D_analytical_patch_geometry(patch_id) - - integer, intent(in) :: patch_id - call s_int_to_str(patch_id, iStr) - - @:PROHIBIT(p == 0, "3D analytical patch "//trim(iStr)//": p must be greater than zero") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "3D analytical patch "//trim(iStr)//": x_centroid must be set") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "3D analytical patch "//trim(iStr)//": y_centroid must be set") - @:PROHIBIT(f_is_default(patch_icpp(patch_id)%z_centroid), "3D analytical patch "//trim(iStr)//": z_centroid must be set") - @:PROHIBIT(patch_icpp(patch_id)%length_x <= 0._wp, "3D analytical patch "//trim(iStr)//": length_x must be greater than zero") - @:PROHIBIT(patch_icpp(patch_id)%length_y <= 0._wp, "3D analytical patch "//trim(iStr)//": length_y must be greater than zero") - @:PROHIBIT(patch_icpp(patch_id)%length_z <= 0._wp, "3D analytical patch "//trim(iStr)//": length_z must be greater than zero") - - end subroutine s_check_3D_analytical_patch_geometry - !> This subroutine checks the model patch input !! @param patch_id Patch identifier impure subroutine s_check_sphere_patch_geometry(patch_id) diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp index 7423d8ed64..15307cb948 100644 --- a/src/pre_process/m_checker.fpp +++ b/src/pre_process/m_checker.fpp @@ -48,6 +48,7 @@ contains impure subroutine s_check_inputs_restart logical :: skip_check !< Flag to skip the check when iterating over !! x, y, and z directions, for special treatment of cylindrical coordinates + integer :: i @:PROHIBIT((.not. old_grid) .and. old_ic, & "old_ic can only be enabled with old_grid enabled") diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index f7700e84ac..010d1e9d4f 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -99,9 +99,6 @@ contains ! Ellipsoidal patch elseif (patch_icpp(i)%geometry == 12) then call s_ellipsoid(i, patch_id_fp, q_prim_vf) - ! Analytical function patch for testing purposes - elseif (patch_icpp(i)%geometry == 13) then - call s_3D_analytical(i, patch_id_fp, q_prim_vf) ! Spherical harmonic patch elseif (patch_icpp(i)%geometry == 14) then call s_spherical_harmonic(i, patch_id_fp, q_prim_vf) @@ -169,9 +166,6 @@ contains elseif (patch_icpp(i)%geometry == 6) then call s_mpi_abort('This used to be the isentropic vortex patch, '// & 'which no longer exists. See Examples. Exiting.') - ! Analytical function patch for testing purposes - elseif (patch_icpp(i)%geometry == 7) then - call s_2D_analytical(i, patch_id_fp, q_prim_vf) ! Spherical Harmonic Patch elseif (patch_icpp(i)%geometry == 14) then call s_spherical_harmonic(i, patch_id_fp, q_prim_vf) @@ -226,9 +220,6 @@ contains if (patch_icpp(i)%geometry == 1) then call s_line_segment(i, patch_id_fp, q_prim_vf) ! 1d analytical - elseif (patch_icpp(i)%geometry == 15) then - call s_1d_analytical(i, patch_id_fp, q_prim_vf) - ! 1d bubble screen with sinusoidal pressure pulse elseif (patch_icpp(i)%geometry == 16) then call s_1d_bubble_pulse(i, patch_id_fp, q_prim_vf) end if @@ -253,13 +244,19 @@ contains integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - real(wp) :: pi_inf, gamma, lit_gamma + ! Generic loop iterators + integer :: i, j, k - integer :: i, j, k !< Generic loop operators + ! Placeholders for the cell boundary values + real(wp) :: pi_inf, gamma, lit_gamma + @:HardcodedDimensionsExtrusion() + @:Hardcoded1DVariables() pi_inf = fluid_pp(1)%pi_inf gamma = fluid_pp(1)%gamma lit_gamma = (1._wp + gamma)/gamma + j = 0 + k = 0 ! Transferring the line segment's centroid and length information x_centroid = patch_icpp(patch_id)%x_centroid @@ -290,11 +287,17 @@ contains @:analytical() + ! check if this should load a hardcoded patch + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded1D() + end if + ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, 0, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, 0, 0) = patch_id end if end do + @:HardcodedDellacation() end subroutine s_line_segment @@ -314,6 +317,8 @@ contains integer :: i, j, k !< Generic loop iterators real(wp) :: th, thickness, nturns, mya real(wp) :: spiral_x_min, spiral_x_max, spiral_y_min, spiral_y_max + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information @@ -353,12 +358,16 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if end do end do + @:HardcodedDellacation() end subroutine s_spiral @@ -381,6 +390,8 @@ contains real(wp) :: radius integer :: i, j, k !< Generic loop iterators + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information @@ -437,13 +448,19 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if + end if end if end do end do + @:HardcodedDellacation() end subroutine s_circle + ! TODO :: Determine if we want analytical and hardcoded patches for the airfoil geometries !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables @@ -792,6 +809,8 @@ contains ! Generic loop iterators integer :: i, j, k real(wp) :: radius, myr, thickness + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information @@ -824,9 +843,12 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id q_prim_vf(alf_idx)%sf(i, j, 0) = patch_icpp(patch_id)%alpha(1)* & exp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp) @@ -834,6 +856,7 @@ contains end do end do + @:HardcodedDellacation() end subroutine s_varcircle @@ -850,6 +873,8 @@ contains ! Generic loop iterators integer :: i, j, k real(wp) :: radius, myr, thickness + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information @@ -887,9 +912,12 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id q_prim_vf(alf_idx)%sf(i, j, k) = patch_icpp(patch_id)%alpha(1)* & exp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp) @@ -898,6 +926,7 @@ contains end do end do end do + @:HardcodedDellacation() end subroutine s_3dvarcircle @@ -916,6 +945,8 @@ contains integer :: i, j, k !< Generic loop operators real(wp) :: a, b + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() ! Transferring the elliptical patch's radii, centroid, smearing ! patch identity, and smearing coefficient information @@ -958,12 +989,16 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if end do end do + @:HardcodedDellacation() end subroutine s_ellipse @@ -984,6 +1019,8 @@ contains ! Generic loop iterators integer :: i, j, k real(wp) :: a, b, c + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the ellipsoidal patch's radii, centroid, smearing ! patch identity, and smearing coefficient information @@ -1038,13 +1075,17 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id end if end do end do end do + @:HardcodedDellacation() end subroutine s_ellipsoid @@ -1069,6 +1110,8 @@ contains integer :: i, j, k !< generic loop iterators real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() pi_inf = fluid_pp(1)%pi_inf gamma = fluid_pp(1)%gamma @@ -1122,6 +1165,10 @@ contains @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if + if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == 4)) then !zero density, reassign according to Tait EOS q_prim_vf(1)%sf(i, j, 0) = & @@ -1130,13 +1177,14 @@ contains end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if end if end if end do end do + @:HardcodedDellacation() end subroutine s_rectangle @@ -1158,6 +1206,8 @@ contains integer :: i, j, k !< Generic loop operators real(wp) :: a, b, c + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the centroid information of the line to be swept x_centroid = patch_icpp(patch_id)%x_centroid @@ -1198,13 +1248,17 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if end do end do + @:HardcodedDellacation() end subroutine s_sweep_line @@ -1224,6 +1278,8 @@ contains integer :: i, j, k !< generic loop iterators real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters real(wp) :: L0, U0 !< Taylor Green Vortex parameters + @:HardcodedDimensionsExtrusion() + @:Hardcoded2DVariables() pi_inf = fluid_pp(1)%pi_inf gamma = fluid_pp(1)%gamma @@ -1268,9 +1324,12 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded2D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id ! Assign Parameters q_prim_vf(mom_idx%beg)%sf(i, j, 0) = U0*sin(x_cc(i)/L0)*cos(y_cc(j)/L0) @@ -1281,73 +1340,9 @@ contains end if end do end do - - end subroutine s_2D_TaylorGreen_Vortex - - !> This patch assigns the primitive variables as analytical - !! functions such that the code can be verified. - !! @param patch_id is the patch identifier - !! @param patch_id_fp Array to track patch ids - !! @param q_prim_vf Array of primitive variables - subroutine s_1D_analytical(patch_id, patch_id_fp, q_prim_vf) - - ! Patch identifier - integer, intent(in) :: patch_id - integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - - ! Generic loop iterators - integer :: i, j, k - - ! Placeholders for the cell boundary values - real(wp) :: pi_inf, gamma, lit_gamma - @:HardcodedDimensionsExtrusion() - @:Hardcoded1DVariables() - - pi_inf = fluid_pp(1)%pi_inf - gamma = fluid_pp(1)%gamma - lit_gamma = (1._wp + gamma)/gamma - j = 0.0_wp - k = 0.0_wp - - ! Transferring the patch's centroid and length information - x_centroid = patch_icpp(patch_id)%x_centroid - length_x = patch_icpp(patch_id)%length_x - - ! Computing the beginning and the end x- and y-coordinates - ! of the patch based on its centroid and lengths - x_boundary%beg = x_centroid - 0.5_wp*length_x - x_boundary%end = x_centroid + 0.5_wp*length_x - - ! Since the patch doesn't allow for its boundaries to be - ! smoothed out, the pseudo volume fraction is set to 1 to - ! ensure that only the current patch contributes to the fluid - ! state in the cells that this patch covers. - eta = 1._wp - - ! Checking whether the line segment covers a particular cell in the - ! domain and verifying whether the current patch has the permission - ! to write to that cell. If both queries check out, the primitive - ! variables of the current patch are assigned to this cell. - do i = 0, m - if (x_boundary%beg <= x_cc(i) .and. & - x_boundary%end >= x_cc(i) .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, 0, 0))) then - - call s_assign_patch_primitive_variables(patch_id, i, 0, 0, & - eta, q_prim_vf, patch_id_fp) - - @:Hardcoded1D() - - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, 0, 0) = patch_id - - end if - end do - @:HardcodedDellacation() - end subroutine s_1D_analytical + end subroutine s_2D_TaylorGreen_Vortex !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids @@ -1365,6 +1360,8 @@ contains integer :: i, j, k ! Placeholders for the cell boundary values real(wp) :: pi_inf, gamma, lit_gamma + @:HardcodedDimensionsExtrusion() + @:Hardcoded1DVariables() pi_inf = fluid_pp(1)%pi_inf gamma = fluid_pp(1)%gamma @@ -1398,169 +1395,15 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded1D() + end if end if end do - - end subroutine s_1D_bubble_pulse - - !> This patch assigns the primitive variables as analytical - !! functions such that the code can be verified. - !! @param patch_id is the patch identifier - !! @param patch_id_fp Array to track patch ids - !! @param q_prim_vf Array of primitive variables - subroutine s_2D_analytical(patch_id, patch_id_fp, q_prim_vf) - - integer, intent(in) :: patch_id - integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - - integer :: i, j, k !< generic loop iterators - - real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters - real(wp) :: l, U0 !< Taylor Green Vortex parameters - @:HardcodedDimensionsExtrusion() - @:Hardcoded2DVariables() - - pi_inf = fluid_pp(1)%pi_inf - gamma = fluid_pp(1)%gamma - lit_gamma = (1._wp + gamma)/gamma - - k = 0.0_wp - - ! Transferring the patch's centroid and length information - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - length_x = patch_icpp(patch_id)%length_x - length_y = patch_icpp(patch_id)%length_y - - ! Computing the beginning and the end x- and y-coordinates - ! of the patch based on its centroid and lengths - x_boundary%beg = x_centroid - 0.5_wp*length_x - x_boundary%end = x_centroid + 0.5_wp*length_x - y_boundary%beg = y_centroid - 0.5_wp*length_y - y_boundary%end = y_centroid + 0.5_wp*length_y - - ! Since the patch doesn't allow for its boundaries to be - ! smoothed out, the pseudo volume fraction is set to 1 to - ! ensure that only the current patch contributes to the fluid - ! state in the cells that this patch covers. - eta = 1._wp - l = 1._wp - U0 = 0.1_wp - ! Checking whether the patch covers a particular cell in the - ! domain and verifying whether the current patch has the - ! permission to write to that cell. If both queries check out, - ! the primitive variables of the current patch are assigned - ! to this cell. - - do j = 0, n - do i = 0, m - if (x_boundary%beg <= x_cc(i) .and. & - x_boundary%end >= x_cc(i) .and. & - y_boundary%beg <= y_cc(j) .and. & - y_boundary%end >= y_cc(j) .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) then - - call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - eta, q_prim_vf, patch_id_fp) - - @:Hardcoded2D() - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, 0) = patch_id - - end if - end do - end do - @:HardcodedDellacation() - end subroutine s_2D_analytical - - !> This patch assigns the primitive variables as analytical - !! functions such that the code can be verified. - !! @param patch_id is the patch identifier - !! @param patch_id_fp Array to track patch ids - !! @param q_prim_vf Array of primitive variables - subroutine s_3D_analytical(patch_id, patch_id_fp, q_prim_vf) - - integer, intent(in) :: patch_id - integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - - integer :: i, j, k !< generic loop iterators - real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters - @:HardcodedDimensionsExtrusion() - @:Hardcoded3DVariables() - - pi_inf = fluid_pp(1)%pi_inf - gamma = fluid_pp(1)%gamma - lit_gamma = (1._wp + gamma)/gamma - - ! Transferring the patch's centroid and length information - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - z_centroid = patch_icpp(patch_id)%z_centroid - length_x = patch_icpp(patch_id)%length_x - length_y = patch_icpp(patch_id)%length_y - length_z = patch_icpp(patch_id)%length_z - - ! Computing the beginning and the end x-, y- and z-coordinates of - ! the patch based on its centroid and lengths - x_boundary%beg = x_centroid - 0.5_wp*length_x - x_boundary%end = x_centroid + 0.5_wp*length_x - y_boundary%beg = y_centroid - 0.5_wp*length_y - y_boundary%end = y_centroid + 0.5_wp*length_y - z_boundary%beg = z_centroid - 0.5_wp*length_z - z_boundary%end = z_centroid + 0.5_wp*length_z - - ! Since the analytical patch does not allow for its boundaries to get - ! smoothed out, the pseudo volume fraction is set to 1 to make sure - ! that only the current patch contributes to the fluid state in the - ! cells that this patch covers. - eta = 1._wp - - ! Checking whether the patch covers a particular cell in the domain - ! and verifying whether the current patch has permission to write to - ! to that cell. If both queries check out, the primitive variables - ! of the current patch are assigned to this cell. - do k = 0, p - do j = 0, n - do i = 0, m - - if (grid_geometry == 3) then - call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) - else - cart_y = y_cc(j) - cart_z = z_cc(k) - end if - - if (x_boundary%beg <= x_cc(i) .and. & - x_boundary%end >= x_cc(i) .and. & - y_boundary%beg <= cart_y .and. & - y_boundary%end >= cart_y .and. & - z_boundary%beg <= cart_z .and. & - z_boundary%end >= cart_z & - .and. & - patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) & - then - - call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) - - @:Hardcoded3D() - - ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id - - end if - - end do - end do - end do - - @:HardcodedDellacation() - end subroutine s_3D_analytical + end subroutine s_1D_bubble_pulse !> This patch generates the shape of the spherical harmonics !! as a perturbation to a perfect sphere @@ -1723,8 +1566,10 @@ contains logical, optional, intent(in) :: ib !< True if this patch is an immersed boundary ! Generic loop iterators - integer :: i, j, k !< generic loop iterators + integer :: i, j, k real(wp) :: radius + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() !! Variables to initialize the pressure field that corresponds to the !! bubble-collapse test case found in Tiwari et al. (2013) @@ -1791,11 +1636,16 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if + end if end if end do end do end do + @:HardcodedDellacation() end subroutine s_sphere @@ -1818,6 +1668,8 @@ contains type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf integer :: i, j, k !< Generic loop iterators + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the cuboid's centroid and length information if (present(ib)) then @@ -1883,9 +1735,12 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id end if end if @@ -1893,6 +1748,7 @@ contains end do end do end do + @:HardcodedDellacation() end subroutine s_cuboid @@ -1917,6 +1773,8 @@ contains integer :: i, j, k !< Generic loop iterators real(wp) :: radius + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the cylindrical patch's centroid, length, radius, ! smoothing patch identity and smoothing coefficient information @@ -2037,14 +1895,18 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id end if end if end do end do end do + @:HardcodedDellacation() end subroutine s_cylinder @@ -2066,6 +1928,8 @@ contains integer :: i, j, k !< Generic loop iterators real(wp) :: a, b, c, d + @:HardcodedDimensionsExtrusion() + @:Hardcoded3DVariables() ! Transferring the centroid information of the plane to be swept x_centroid = patch_icpp(patch_id)%x_centroid @@ -2119,14 +1983,18 @@ contains eta, q_prim_vf, patch_id_fp) @:analytical() + if (patch_icpp(patch_id)%hcid /= dflt_int) then + @:Hardcoded3D() + end if ! Updating the patch identities bookkeeping variable - if (1._wp - eta < 1.e-16_wp) patch_id_fp(i, j, k) = patch_id + if (1._wp - eta < sgm_eps) patch_id_fp(i, j, k) = patch_id end if end do end do end do + @:HardcodedDellacation() end subroutine s_sweep_plane diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 34c422f02a..a006e25291 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -96,6 +96,8 @@ def __is_ic_analytical(self, key: str, val: str) -> bool: # pylint: disable=too-many-locals def __get_pre_fpp(self, print: bool) -> str: + # generates the content of an FFP file that will hold the functions for + # some initial condition DATA = { 1: {'ptypes': [1, 15, 16], 'sf_idx': 'i, 0, 0'}, 2: {'ptypes': [2, 3, 4, 5, 6, 7, 17, 18, 21], 'sf_idx': 'i, j, 0'}, @@ -104,6 +106,8 @@ def __get_pre_fpp(self, print: bool) -> str: patches = {} + # iterates over the parameters and checks if they are defined as an + # analytical function. If so, append it to the `patches`` object for key, val in self.params.items(): if not self.__is_ic_analytical(key, val): continue @@ -117,12 +121,16 @@ def __get_pre_fpp(self, print: bool) -> str: srcs = [] + # for each analytical patch that is required to be added, generate + # the string that contains that function. for pid, items in patches.items(): ptype = self.params[f"patch_icpp({pid})%geometry"] if ptype not in DATA['ptypes']: raise common.MFCException(f"Patch #{pid} of type {ptype} cannot be analytically defined.") + # function that defines how we will replace variable names with + # values from the case file def rhs_replace(match): return { 'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)', @@ -137,6 +145,8 @@ def rhs_replace(match): }.get(match.group(), match.group()) lines = [] + # perform the replacement of strings for each analytic function + # to generate some fortran string representing the code passed in for attribute, expr in items: if print: cons.print(f"* Codegen: {attribute} = {expr}") @@ -153,6 +163,9 @@ def rhs_replace(match): lines.append(f" {lhs} = {rhs}") + # concatenates all of the analytic lines into a single string with + # each element separated by new line characters. Then write those + # new lines as a fully concatenated string with fortran syntax srcs.append(f"""\ if (patch_id == {pid}) then {f'{chr(10)}'.join(lines)} @@ -168,7 +181,6 @@ def rhs_replace(match): {f'{chr(10)}{chr(10)}'.join(srcs)} #:enddef """ - return content def __get_sim_fpp(self, print: bool) -> str: diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 6e1027d9c7..0ff6ccf490 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -935,9 +935,26 @@ def foreach_example(): if path == "scaling": continue - # # List of currently broken examples -> currently attempting to fix! - brokenCases = ["2D_ibm_cfl_dt", "1D_sodHypo", "2D_viscous", "2D_laplace_pressure_jump", "2D_bubbly_steady_shock", "2D_advection", "2D_hardcodied_ic", "2D_ibm_multiphase", "2D_acoustic_broadband", "1D_inert_shocktube", "1D_reactive_shocktube", "2D_ibm_steady_shock", "3D_performance_test", "3D_ibm_stl_ellipsoid", "3D_sphbubcollapse", "2D_ibm_stl_wedge", "3D_ibm_stl_pyramid", "3D_ibm_bowshock", "3D_turb_mixing", "2D_mixing_artificial_Ma", "2D_lagrange_bubblescreen", "3D_lagrange_bubblescreen", "2D_triple_point"] - if path in brokenCases: + # # List of all example cases that will be skipped during testing + casesToSkip = ["2D_ibm_cfl_dt", "1D_sodHypo", "2D_viscous", + "2D_laplace_pressure_jump", "2D_bubbly_steady_shock", + "2D_advection", "2D_hardcodied_ic", + "2D_ibm_multiphase", "2D_acoustic_broadband", + "1D_inert_shocktube", "1D_reactive_shocktube", + "2D_ibm_steady_shock", "3D_performance_test", + "3D_ibm_stl_ellipsoid", "3D_sphbubcollapse", + "2D_ibm_stl_wedge", "3D_ibm_stl_pyramid", + "3D_ibm_bowshock", "3D_turb_mixing", + "2D_mixing_artificial_Ma", + "2D_lagrange_bubblescreen", + "3D_lagrange_bubblescreen", "2D_triple_point", + "1D_shuosher_analytical", + "1D_titarevtorro_analytical", + "2D_acoustic_pulse_analytical", + "2D_isentropicvortex_analytical", + "2D_zero_circ_vortex_analytical", + "3D_TaylorGreenVortex_analytical"] + if path in casesToSkip: continue name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}" path = os.path.join(common.MFC_EXAMPLE_DIRPATH, path, "case.py")