diff --git a/examples/2D_IGR_forward_facing_step/README.md b/examples/2D_IGR_forward_facing_step/README.md
new file mode 100644
index 0000000000..3b85c28d3a
--- /dev/null
+++ b/examples/2D_IGR_forward_facing_step/README.md
@@ -0,0 +1,7 @@
+# Forward Facing Step With IGR (2D)
+
+Reference: See Section IV, b.
+> Woodward, P. *(1984). The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54(1), 115–173. https://doi.org/10.1016/0021-9991(84)90140-2*
+
+## Evolved State
+
diff --git a/examples/2D_IGR_forward_facing_step/case.py b/examples/2D_IGR_forward_facing_step/case.py
new file mode 100644
index 0000000000..1b8585f0a6
--- /dev/null
+++ b/examples/2D_IGR_forward_facing_step/case.py
@@ -0,0 +1,97 @@
+import json
+import math
+
+h = 0.2
+
+# Radius as a percentage of height (h)
+rc = 0.2
+
+gam_a = 1.4
+p0 = 1
+rho0 = 1.4
+c0 = math.sqrt(gam_a * p0 / rho0)
+v0 = 3 * c0
+mu = rho0 * v0 * h / 2e5
+
+# Configuring case dictionary
+print(
+ json.dumps(
+ {
+ # Logistics
+ "run_time_info": "T",
+ "x_domain%beg": 0,
+ "x_domain%end": 15 * h,
+ "y_domain%beg": 0,
+ "y_domain%end": 5 * h,
+ "cyl_coord": "F",
+ "m": 1499,
+ "n": 499,
+ "p": 0,
+ "cfl_adap_dt": "T",
+ "cfl_target": 0.6,
+ "n_start": 0,
+ "t_save": 0.04,
+ "t_stop": 4,
+ # Simulation Algorithm Parameters
+ "num_patches": 1,
+ "model_eqns": 2,
+ "alt_soundspeed": "F",
+ "num_fluids": 1,
+ "mpp_lim": "F",
+ "mixture_err": "F",
+ "time_stepper": 3,
+ "igr": "T",
+ "igr_pres_lim": "T",
+ "igr_order": 3,
+ "igr_iter_solver": 1,
+ "num_igr_iters": 5,
+ "num_igr_warm_start_iters": 50,
+ "bc_x%beg": -3,
+ "bc_x%end": -3,
+ "bc_y%beg": -2,
+ "bc_y%end": -2,
+ "ib": "T",
+ "num_ibs": 3,
+ # Formatted Database Files Structure Parameters
+ "format": 1,
+ "precision": 2,
+ "prim_vars_wrt": "T",
+ "parallel_io": "T",
+ # Patch 1 Background
+ "patch_icpp(1)%geometry": 3,
+ "patch_icpp(1)%x_centroid": 7.5 * h,
+ "patch_icpp(1)%y_centroid": 2.5 * h,
+ "patch_icpp(1)%length_x": 15 * h,
+ "patch_icpp(1)%length_y": 5 * h,
+ "patch_icpp(1)%vel(1)": v0,
+ "patch_icpp(1)%vel(2)": 0.0,
+ "patch_icpp(1)%pres": p0,
+ "patch_icpp(1)%alpha_rho(1)": rho0,
+ "patch_icpp(1)%alpha(1)": 1.0,
+ # Patch: Slip rectangle with rounded corner
+ "patch_ib(1)%geometry": 3,
+ "patch_ib(1)%x_centroid": (9 + rc / 2) * h,
+ "patch_ib(1)%y_centroid": 0.5 * h,
+ "patch_ib(1)%length_x": (12 - rc) * h,
+ "patch_ib(1)%length_y": h,
+ "patch_ib(1)%slip": "T",
+ "patch_ib(2)%geometry": 3,
+ "patch_ib(2)%x_centroid": (3 + rc / 2) * h,
+ "patch_ib(2)%y_centroid": (0.5 - rc / 2) * h,
+ "patch_ib(2)%length_x": rc * h,
+ "patch_ib(2)%length_y": (1 - rc) * h,
+ "patch_ib(2)%slip": "T",
+ "patch_ib(3)%geometry": 2,
+ "patch_ib(3)%x_centroid": (3 + rc) * h,
+ "patch_ib(3)%y_centroid": (1 - rc) * h,
+ "patch_ib(3)%radius": rc * h,
+ "patch_ib(3)%slip": "T",
+ # Fluids Physical Parameters
+ "fluid_pp(1)%gamma": 1.0 / (gam_a - 1.0),
+ "fluid_pp(1)%pi_inf": 0.0,
+ "viscous": "T",
+ "fluid_pp(1)%Re(1)": 1 / mu,
+ },
+ indent=4,
+ )
+)
diff --git a/examples/2D_IGR_forward_facing_step/figure.png b/examples/2D_IGR_forward_facing_step/figure.png
new file mode 100644
index 0000000000..b2bde6dc15
Binary files /dev/null and b/examples/2D_IGR_forward_facing_step/figure.png differ
diff --git a/examples/3D_IGR_bowshock/README.md b/examples/3D_IGR_bowshock/README.md
new file mode 100644
index 0000000000..5cc7ea5916
--- /dev/null
+++ b/examples/3D_IGR_bowshock/README.md
@@ -0,0 +1,7 @@
+# Mach 2 Flow Over Sphere (bowshock) With IGR (3D)
+
+Reference: See Section 5.4
+> Uddin, H. *(2014). A Cartesian-based embedded geometry technique with adaptive high-order finite differences for compressible flow around complex geometries. Journal of Computational Physics, 262, 379–407. https://doi.org/10.1016/j.jcp.2014.01.004*
+
+## Figure
+
diff --git a/examples/3D_IGR_bowshock/case.py b/examples/3D_IGR_bowshock/case.py
new file mode 100644
index 0000000000..da30ead932
--- /dev/null
+++ b/examples/3D_IGR_bowshock/case.py
@@ -0,0 +1,101 @@
+import json
+import math
+
+gam = 1.4
+
+D = 0.1
+N = 400
+
+tf = 0.01
+saveFreq = tf / 200
+
+# Configuring case dictionary
+print(
+ json.dumps(
+ {
+ # Logistics
+ "run_time_info": "T",
+ # Computational Domain Parameters
+ # x direction
+ "x_domain%beg": -5.0 * D,
+ "x_domain%end": 5.0 * D,
+ # y direction
+ "y_domain%beg": -2.5 * D,
+ "y_domain%end": 2.5 * D,
+ # z direction
+ "z_domain%beg": -2.5 * D,
+ "z_domain%end": 2.5 * D,
+ "cyl_coord": "F",
+ "m": 2 * N,
+ "n": N,
+ "p": N,
+ "cfl_adap_dt": "T",
+ "cfl_target": 0.4,
+ "n_start": 0,
+ "t_save": saveFreq,
+ "t_stop": tf,
+ # Simulation Algorithm Parameters
+ "num_patches": 1,
+ # Use the 5 equation model
+ "model_eqns": 2,
+ "alt_soundspeed": "F",
+ "num_fluids": 1,
+ "mpp_lim": "F",
+ "mixture_err": "T",
+ "time_stepper": 3,
+ # Use IGR5
+ "igr": "T",
+ "igr_pres_lim": "T",
+ "igr_order": 5,
+ "alf_factor": 10,
+ "igr_iter_solver": 2,
+ "num_igr_iters": 15,
+ "num_igr_warm_start_iters": 50,
+ # We use ghost-cell extrapolation
+ "bc_x%beg": -3,
+ "bc_x%end": -3,
+ "bc_y%beg": -3,
+ "bc_y%end": -3,
+ "bc_z%beg": -3,
+ "bc_z%end": -3,
+ # Set IB to True
+ "ib": "T",
+ "num_ibs": 1,
+ "viscous": "T",
+ # Formatted Database Files Structure Parameters
+ "format": 1,
+ "precision": 2,
+ "prim_vars_wrt": "T",
+ "E_wrt": "T",
+ "parallel_io": "T",
+ # Main Patch: 10D:5D:5D rectangular prism centered at the (0, 0, 0)
+ # HCID 390 smooths the x-velocity profile in a small region around
+ # the sphere using a tanh profile.
+ "patch_icpp(1)%geometry": 9,
+ "patch_icpp(1)%hcid": 390,
+ "patch_icpp(1)%x_centroid": 0.0,
+ "patch_icpp(1)%y_centroid": 0.0,
+ "patch_icpp(1)%z_centroid": 0.0,
+ "patch_icpp(1)%length_x": 10 * D,
+ "patch_icpp(1)%length_y": 5 * D,
+ "patch_icpp(1)%length_z": 5 * D,
+ "patch_icpp(1)%vel(1)": 527.2e00,
+ "patch_icpp(1)%vel(2)": 0.0e00,
+ "patch_icpp(1)%vel(3)": 0.0e00,
+ "patch_icpp(1)%pres": 10918.2549,
+ "patch_icpp(1)%alpha_rho(1)": 0.2199,
+ "patch_icpp(1)%alpha(1)": 1.0e00,
+ # Patch: Sphere Immersed Boundary
+ "patch_ib(1)%geometry": 8,
+ "patch_ib(1)%x_centroid": -3.0 * D,
+ "patch_ib(1)%y_centroid": 0.0,
+ "patch_ib(1)%z_centroid": 0.0,
+ "patch_ib(1)%radius": D / 2,
+ "patch_ib(1)%slip": "T",
+ # Fluids Physical Parameters
+ "fluid_pp(1)%gamma": 1.0e00 / (gam - 1.0e00),
+ "fluid_pp(1)%pi_inf": 0,
+ "fluid_pp(1)%Re(1)": 650000,
+ }
+ )
+)
diff --git a/examples/3D_IGR_bowshock/figure.png b/examples/3D_IGR_bowshock/figure.png
new file mode 100644
index 0000000000..75076cc7c8
Binary files /dev/null and b/examples/3D_IGR_bowshock/figure.png differ
diff --git a/src/common/include/3dHardcodedIC.fpp b/src/common/include/3dHardcodedIC.fpp
index 40d5f03ad3..6e7a9607de 100644
--- a/src/common/include/3dHardcodedIC.fpp
+++ b/src/common/include/3dHardcodedIC.fpp
@@ -189,6 +189,27 @@
q_prim_vf(momxb + 1)%sf(i, j, k) = -Mach*376.636429464809*cos(x_cc(i)/1)*sin(y_cc(j)/1)*sin(z_cc(k)/1)
end if
+ case (390)
+ ! This is for smoothing the x-velocity in a small neighborhood around the sphere for the bowshock case
+ r = (x_cc(i) + 0.3_wp)*(x_cc(i) + 0.3_wp) + y_cc(j)*y_cc(j) + z_cc(k)*z_cc(k)
+ r = sqrt(r)
+ ! Size of smoothing region
+ rcut = 0.025_wp
+
+ if (r < 0.05 + rcut + 0.01) then
+ if (r < 0.05_wp - 0.005) then
+ ! If inside sphere, set velocity to 0
+ q_prim_vf(momxb)%sf(i, j, k) = 0.0_wp
+ else if (r > 0.05_wp + rcut + 0.005) then
+ ! If outside the smoothing region, do nothing
+ else
+ ! If inside the smoothing region, then interpolate values via tanh function
+ r = (r - 0.05_wp)/rcut
+ q_prim_vf(momxb)%sf(i, j, k) = 0.0_wp + (527.2_wp - 0.0_wp)* &
+ (0.5_wp*(1.0_wp + tanh(5.0_wp*(r - 0.5_wp))))
+ end if
+ 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_ib_patches.fpp b/src/pre_process/m_check_ib_patches.fpp
index 1b09929dda..213e966a02 100644
--- a/src/pre_process/m_check_ib_patches.fpp
+++ b/src/pre_process/m_check_ib_patches.fpp
@@ -37,6 +37,9 @@ contains
integer :: i
+ @:PROHIBIT(igr .and. any(patch_ib(:)%moving_ibm > 0), "Cannot use &
+ moving immersed boundary with IGR. All patch_ib(:)%moving_ibm must be 0.")
+
do i = 1, num_patches_max
if (i <= num_ibs) then
! call s_check_patch_geometry(i)
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index b634ec0078..940e9866cc 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -38,6 +38,7 @@ module m_ibm
; public :: s_initialize_ibm_module, &
s_ibm_setup, &
s_ibm_correct_state, &
+ s_interpolate_sigma_igr, &
s_finalize_ibm_module
integer, allocatable, dimension(:, :, :) :: patch_id_fp
@@ -154,6 +155,51 @@ contains
end subroutine s_populate_ib_buffers
+ !> Interpolates sigma from the m_igr module at all ghost points
+ !! @param jac: Sigma, Entropic pressure
+ subroutine s_interpolate_sigma_igr(jac)
+ real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:), intent(inout) :: jac
+ integer :: j, k, l, r, s, t, i
+ integer :: j1, j2, k1, k2, l1, l2
+ real(wp) :: coeff, jac_IP
+ type(ghost_point) :: gp
+
+ ! At all ghost points, use its image point to interpolate sigma
+ if (num_gps > 0) then
+ $:GPU_PARALLEL_LOOP(private='[i, j, k, l, j1, j2, k1, k2, l1, l2, r, s, t, gp, coeff, jac_IP]')
+ do i = 1, num_gps
+ jac_IP = 0._wp
+ gp = ghost_points(i)
+ r = gp%loc(1)
+ s = gp%loc(2)
+ t = gp%loc(3)
+
+ j1 = gp%ip_grid(1); j2 = j1 + 1
+ k1 = gp%ip_grid(2); k2 = k1 + 1
+ l1 = gp%ip_grid(3); l2 = l1 + 1
+
+ if (p == 0) then
+ l1 = 0
+ l2 = 0
+ end if
+
+ $:GPU_LOOP(parallelism='[seq]')
+ do l = l1, l2
+ $:GPU_LOOP(parallelism='[seq]')
+ do k = k1, k2
+ $:GPU_LOOP(parallelism='[seq]')
+ do j = j1, j2
+ coeff = gp%interp_coeffs(j - j1 + 1, k - k1 + 1, l - l1 + 1)
+ jac_IP = jac_IP + coeff*jac(j, k, l)
+ end do
+ end do
+ end do
+ jac(r, s, t) = jac_IP
+ end do
+ $:END_GPU_PARALLEL_LOOP()
+ end if
+ end subroutine s_interpolate_sigma_igr
+
!> Subroutine that updates the conservative variables at the ghost points
!! @param q_cons_vf Conservative Variables
!! @param q_prim_vf Primitive variables
@@ -163,7 +209,7 @@ contains
type(scalar_field), &
dimension(sys_size), &
- intent(INOUT) :: q_cons_vf !< Primitive Variables
+ intent(INOUT) :: q_cons_vf !< Conservative Variables
type(scalar_field), &
dimension(sys_size), &
@@ -200,31 +246,33 @@ contains
type(ghost_point) :: gp
type(ghost_point) :: innerp
- ! set the Moving IBM interior Pressure Values
- $:GPU_PARALLEL_LOOP(private='[i,j,k,patch_id,rho]', copyin='[E_idx,momxb]', collapse=3)
- do l = 0, p
- do k = 0, n
- do j = 0, m
- patch_id = ib_markers%sf(j, k, l)
- if (patch_id /= 0) then
- q_prim_vf(E_idx)%sf(j, k, l) = 1._wp
- if (patch_ib(patch_id)%moving_ibm > 0) then
- rho = 0._wp
- do i = 1, num_fluids
- rho = rho + q_prim_vf(contxb + i - 1)%sf(j, k, l)
- end do
+ if (.not. igr) then
+ ! set the Moving IBM interior Pressure Values
+ $:GPU_PARALLEL_LOOP(private='[i,j,k,patch_id,rho]', copyin='[E_idx,momxb]', collapse=3)
+ do l = 0, p
+ do k = 0, n
+ do j = 0, m
+ patch_id = ib_markers%sf(j, k, l)
+ if (patch_id /= 0) then
+ q_prim_vf(E_idx)%sf(j, k, l) = 1._wp
+ if (patch_ib(patch_id)%moving_ibm > 0) then
+ rho = 0._wp
+ do i = 1, num_fluids
+ rho = rho + q_prim_vf(contxb + i - 1)%sf(j, k, l)
+ end do
- ! Sets the momentum
- do i = 1, num_dims
- q_cons_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
- q_prim_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
- end do
+ ! Sets the momentum
+ do i = 1, num_dims
+ q_cons_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
+ q_prim_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
+ end do
+ end if
end if
- end if
+ end do
end do
end do
- end do
- $:END_GPU_PARALLEL_LOOP()
+ $:END_GPU_PARALLEL_LOOP()
+ end if
if (num_gps > 0) then
$:GPU_PARALLEL_LOOP(private='[i,physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, radial_vector, rotation_velocity, j,k,l,q,qv_K,c_IP,nbub,patch_id]')
@@ -256,34 +304,52 @@ contains
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
+ else if (igr) then
+ call s_interpolate_image_point(q_prim_vf, gp, &
+ alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, q_cons_vf=q_cons_vf, G_K=G_K, Gs=Gs)
else
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
end if
dyn_pres = 0._wp
-
- ! Set q_prim_vf params at GP so that mixture vars calculated properly
- $:GPU_LOOP(parallelism='[seq]')
- do q = 1, num_fluids
- q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
- q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
- end do
+ if (igr) then
+ if (num_fluids == 1) then
+ q_cons_vf(1)%sf(j, k, l) = alpha_rho_IP(1)
+ else
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = 1, num_fluids - 1
+ q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
+ q_cons_vf(E_idx + q)%sf(j, k, l) = alpha_IP(q)
+ end do
+ q_cons_vf(num_fluids)%sf(j, k, l) = alpha_rho_IP(num_fluids)
+ end if
+ else
+ ! Set q_prim_vf params at GP so that mixture vars calculated properly
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = 1, num_fluids
+ q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
+ q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
+ end do
+ end if
if (surface_tension) then
q_prim_vf(c_idx)%sf(j, k, l) = c_IP
end if
! set the pressure
- if (patch_ib(patch_id)%moving_ibm <= 1) then
- q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
- else
- q_prim_vf(E_idx)%sf(j, k, l) = 0._wp
- $:GPU_LOOP(parallelism='[seq]')
- do q = 1, num_fluids
- ! Se the pressure inside a moving immersed boundary based upon the pressure of the image point. acceleration, and normal vector direction
- q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) + pres_IP/(1._wp - 2._wp*abs(levelset%sf(j, k, l, patch_id)*alpha_rho_IP(q)/pres_IP)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j, k, l, patch_id, :)))
- end do
+ ! !TEMPORARY, NEED TO FIX FOR MOVING IGR
+ if (.not. igr) then
+ if (patch_ib(patch_id)%moving_ibm <= 1) then
+ q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
+ else
+ q_prim_vf(E_idx)%sf(j, k, l) = 0._wp
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = 1, num_fluids
+ ! Se the pressure inside a moving immersed boundary based upon the pressure of the image point. acceleration, and normal vector direction
+ q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) + pres_IP/(1._wp - 2._wp*abs(levelset%sf(j, k, l, patch_id)*alpha_rho_IP(q)/pres_IP)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j, k, l, patch_id, :)))
+ end do
+ end if
end if
if (model_eqns /= 4) then
@@ -339,12 +405,14 @@ contains
vel_g(q - momxb + 1)/2._wp
end do
- ! Set continuity and adv vars
- $:GPU_LOOP(parallelism='[seq]')
- do q = 1, num_fluids
- q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
- q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
- end do
+ if (.not. igr) then
+ ! Set continuity and adv vars
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = 1, num_fluids
+ q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
+ q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
+ end do
+ end if
! Set color function
if (surface_tension) then
@@ -357,6 +425,7 @@ contains
else
q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
end if
+
! Set bubble vars
if (bubbles_euler .and. .not. qbmm) then
call s_comp_n_from_prim(alpha_IP(1), r_IP, nbub, weight)
@@ -844,15 +913,21 @@ contains
!! at the cell centers in order to estimate the state at the image point
subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, &
pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, &
- mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
+ mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP, q_cons_vf, G_K, Gs)
$:GPU_ROUTINE(parallelism='[seq]')
type(scalar_field), &
dimension(sys_size), &
intent(IN) :: q_prim_vf !< Primitive Variables
+ type(scalar_field), optional, &
+ dimension(sys_size), &
+ intent(IN) :: q_cons_vf !< Conservative Variables
real(stp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(IN) :: pb_in, mv_in
type(ghost_point), intent(IN) :: gp
+ real(wp), optional, intent(INOUT) :: G_K
+ real(wp), optional, dimension(num_fluids), intent(IN) :: Gs
+
real(wp), intent(INOUT) :: pres_IP
real(wp), dimension(3), intent(INOUT) :: vel_IP
real(wp), intent(INOUT) :: c_IP
@@ -864,6 +939,12 @@ contains
integer :: i, j, k, l, q !< Iterator variables
integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
real(wp) :: coeff
+ real(wp) :: alpha_sum
+ real(wp) :: pres, dyn_pres, pres_mag, T
+ real(wp) :: rhoYks(1:num_species)
+ real(wp) :: rho_K, gamma_K, pi_inf_K, qv_K
+ real(wp), dimension(num_fluids) :: alpha_K, alpha_rho_K
+ real(wp), dimension(2) :: Re_K
i1 = gp%ip_grid(1); i2 = i1 + 1
j1 = gp%ip_grid(2); j2 = j1 + 1
@@ -878,6 +959,7 @@ contains
alpha_IP = 0._wp
pres_IP = 0._wp
vel_IP = 0._wp
+ pres = 0._wp
if (surface_tension) c_IP = 0._wp
@@ -904,31 +986,100 @@ contains
do j = j1, j2
$:GPU_LOOP(parallelism='[seq]')
do k = k1, k2
-
coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1)
- pres_IP = pres_IP + coeff* &
- q_prim_vf(E_idx)%sf(i, j, k)
+ if (igr) then
+ ! For IGR, we will need to perform operations on
+ ! the conservative variables instead
+ alpha_sum = 0._wp
+ dyn_pres = 0._wp
+ if (num_fluids == 1) then
+ alpha_rho_K(1) = q_cons_vf(1)%sf(i, j, k)
+ alpha_K(1) = 1._wp
+ else
+ $:GPU_LOOP(parallelism='[seq]')
+ do l = 1, num_fluids - 1
+ alpha_rho_K(l) = q_cons_vf(l)%sf(i, j, k)
+ alpha_K(l) = q_cons_vf(E_idx + l)%sf(i, j, k)
+ end do
+ alpha_rho_K(num_fluids) = q_cons_vf(num_fluids)%sf(i, j, k)
+ alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1))
+ end if
+ if (model_eqns /= 4) then
+ if (elasticity) then
+ ! NOTE: NEED TO VERIFY ELASTICITY W/ IGR
+ call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &
+ alpha_rho_K, Re_K, G_K, Gs)
+ else
+ call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
+ alpha_K, alpha_rho_K, Re_K)
+ end if
+ end if
- $:GPU_LOOP(parallelism='[seq]')
- do q = momxb, momxe
- vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
- q_prim_vf(q)%sf(i, j, k)
- end do
+ rho_K = max(rho_K, sgm_eps)
- $:GPU_LOOP(parallelism='[seq]')
- do l = contxb, contxe
- alpha_rho_IP(l) = alpha_rho_IP(l) + coeff* &
- q_prim_vf(l)%sf(i, j, k)
- alpha_IP(l) = alpha_IP(l) + coeff* &
- q_prim_vf(advxb + l - 1)%sf(i, j, k)
- end do
+ $:GPU_LOOP(parallelism='[seq]')
+ do l = momxb, momxe
+ if (model_eqns /= 4) then
+ dyn_pres = dyn_pres + 5.e-1_wp*q_cons_vf(l)%sf(i, j, k) &
+ *q_cons_vf(l)%sf(i, j, k)/rho_K
+ end if
+ end do
+ pres_mag = 0._wp
+ T = 0._wp
+ rhoYks = 0._wp
+
+ call s_compute_pressure(q_cons_vf(E_idx)%sf(i, j, k), &
+ q_cons_vf(alf_idx)%sf(i, j, k), &
+ dyn_pres, pi_inf_K, gamma_K, rho_K, &
+ qv_K, rhoYks, pres, T, pres_mag=pres_mag)
+
+ pres_IP = pres_IP + coeff*pres
+
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = momxb, momxe
+ vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
+ q_cons_vf(q)%sf(i, j, k)/rho_K
+ end do
+
+ if (num_fluids == 1) then
+ alpha_rho_IP(1) = alpha_rho_IP(1) + coeff*q_cons_vf(contxb)%sf(i, j, k)
+ alpha_IP(1) = alpha_IP(1) + coeff*1._wp
+ else
+ alpha_sum = 0._wp
+ $:GPU_LOOP(parallelism='[seq]')
+ do l = 1, num_fluids - 1
+ alpha_rho_IP(l) = alpha_rho_IP(l) + coeff*q_cons_vf(l)%sf(i, j, k)
+ alpha_IP(l) = alpha_IP(l) + coeff*q_cons_vf(E_idx + l)%sf(i, j, k)
+ alpha_sum = alpha_sum + coeff*q_cons_vf(E_idx + l)%sf(i, j, k)
+ end do
+ alpha_rho_IP(num_fluids) = alpha_rho_IP(num_fluids) + coeff*q_cons_vf(num_fluids)%sf(i, j, k)
+ alpha_IP(num_fluids) = alpha_IP(num_fluids) + (coeff*1._wp - alpha_sum)
+ end if
+ else
+ pres_IP = pres_IP + coeff* &
+ q_prim_vf(E_idx)%sf(i, j, k)
+
+ $:GPU_LOOP(parallelism='[seq]')
+ do q = momxb, momxe
+ vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
+ q_prim_vf(q)%sf(i, j, k)
+ end do
- if (surface_tension) then
+ $:GPU_LOOP(parallelism='[seq]')
+ do l = contxb, contxe
+ alpha_rho_IP(l) = alpha_rho_IP(l) + coeff* &
+ q_prim_vf(l)%sf(i, j, k)
+ alpha_IP(l) = alpha_IP(l) + coeff* &
+ q_prim_vf(advxb + l - 1)%sf(i, j, k)
+ end do
+ end if
+
+ if (surface_tension .and. .not. igr) then
c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
end if
- if (bubbles_euler .and. .not. qbmm) then
+ if (bubbles_euler .and. .not. qbmm .and. .not. igr) then
$:GPU_LOOP(parallelism='[seq]')
do l = 1, nb
if (polytropic) then
@@ -943,7 +1094,7 @@ contains
end do
end if
- if (qbmm) then
+ if (qbmm .and. .not. igr) then
do l = 1, nb*nmom
nmom_IP(l) = nmom_IP(l) + coeff*q_prim_vf(bubxb - 1 + l)%sf(i, j, k)
end do
@@ -957,9 +1108,7 @@ contains
end do
end do
end if
-
end if
-
end do
end do
end do
diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp
index 60555b554b..301a3054a7 100644
--- a/src/simulation/m_igr.fpp
+++ b/src/simulation/m_igr.fpp
@@ -15,6 +15,8 @@ module m_igr
use m_boundary_common
+ use m_ibm, only: ib_markers, s_interpolate_sigma_igr
+
implicit none
private; public :: s_initialize_igr_module, &
@@ -368,7 +370,32 @@ contains
$:END_GPU_PARALLEL_LOOP()
end if
end do
-
+ if (ib) then
+ $:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
+ do l = 0, p
+ do k = 0, n
+ do j = 0, m
+ if (ib_markers%sf(j, k, l) /= 0) then
+ jac(j, k, l) = 0._wp
+ end if
+ end do
+ end do
+ end do
+ $:END_GPU_PARALLEL_LOOP()
+ call s_interpolate_sigma_igr(jac)
+ ! If using Jacobi Iter, we need to update jac_old again
+ if (igr_iter_solver == 1) then
+ $:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
+ do l = idwbuff(3)%beg, idwbuff(3)%end
+ do k = idwbuff(2)%beg, idwbuff(2)%end
+ do j = idwbuff(1)%beg, idwbuff(1)%end
+ jac_old(j, k, l) = jac(j, k, l)
+ end do
+ end do
+ end do
+ $:END_GPU_PARALLEL_LOOP()
+ end if
+ end if
end subroutine s_igr_iterative_solve
subroutine s_igr_sigma_x(q_cons_vf, rhs_vf)
diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py
index d52cce412a..7f062f2f59 100644
--- a/toolchain/mfc/case_validator.py
+++ b/toolchain/mfc/case_validator.py
@@ -706,7 +706,6 @@ def check_igr_simulation(self): # pylint: disable=too-many-locals
igr_iter_solver = self.get('igr_iter_solver')
alf_factor = self.get('alf_factor')
model_eqns = self.get('model_eqns')
- ib = self.get('ib', 'F') == 'T'
bubbles_euler = self.get('bubbles_euler', 'F') == 'T'
bubbles_lagrange = self.get('bubbles_lagrange', 'F') == 'T'
alt_soundspeed = self.get('alt_soundspeed', 'F') == 'T'
@@ -718,6 +717,7 @@ def check_igr_simulation(self): # pylint: disable=too-many-locals
hyperelasticity = self.get('hyperelasticity', 'F') == 'T'
cyl_coord = self.get('cyl_coord', 'F') == 'T'
probe_wrt = self.get('probe_wrt', 'F') == 'T'
+ chemistry = self.get('chemistry', 'F') == 'T'
self.prohibit(num_igr_iters is not None and num_igr_iters < 0,
"num_igr_iters must be greater than or equal to 0")
@@ -729,8 +729,6 @@ def check_igr_simulation(self): # pylint: disable=too-many-locals
"alf_factor must be non-negative")
self.prohibit(model_eqns is not None and model_eqns != 2,
"IGR only supports model_eqns = 2")
- self.prohibit(ib,
- "IGR does not support the immersed boundary method")
self.prohibit(bubbles_euler,
"IGR does not support Euler-Euler bubble models")
self.prohibit(bubbles_lagrange,
@@ -753,6 +751,9 @@ def check_igr_simulation(self): # pylint: disable=too-many-locals
"IGR does not support cylindrical or axisymmetric coordinates")
self.prohibit(probe_wrt,
"IGR does not support probe writes")
+ self.prohibit(chemistry,
+ "IGR does not support chemistry")
+
# Check BCs - IGR does not support characteristic BCs
# Characteristic BCs are BC_CHAR_SLIP_WALL (-5) through BC_CHAR_SUP_OUTFLOW (-12)
diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py
index 228c24ac1c..eb62b23057 100644
--- a/toolchain/mfc/test/cases.py
+++ b/toolchain/mfc/test/cases.py
@@ -1030,8 +1030,10 @@ def foreach_example():
"3D_IGR_TaylorGreenVortex_nvidia",
"2D_backward_facing_step",
"2D_forward_facing_step",
+ "2D_IGR_forward_facing_step",
"1D_convergence",
- "3D_IGR_33jet","1D_multispecies_diffusion"]
+ "3D_IGR_33jet","1D_multispecies_diffusion",
+ "3D_IGR_bowshock"]
if path in casesToSkip:
continue
name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}"