From af75f84fcfdd6740c57695f7dd9c48878428bd61 Mon Sep 17 00:00:00 2001 From: Mark Zhang Date: Thu, 18 Dec 2025 18:30:19 -0500 Subject: [PATCH 1/5] Integrate IBM with IGR --- src/simulation/m_ibm.fpp | 236 +++++++++++++++++++++++++------- src/simulation/m_igr.fpp | 28 +++- toolchain/mfc/case_validator.py | 4 +- 3 files changed, 212 insertions(+), 56 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 048032b9ae..c8eb33ce7a 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -150,6 +150,42 @@ contains end subroutine s_populate_ib_buffers + subroutine s_update_igr(jac_sf) + type(scalar_field), dimension(1), intent(inout) :: jac_sf + integer :: j, k, l, r, s, t, i + integer :: j1, j2, k1, k2, l1, l2 + real(wp) :: coeff, jac_IP + type(ghost_point) :: gp + 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 + + $: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_sf(1)%sf(j, k, l) + end do + end do + end do + jac_sf(1)%sf(r, s, t) = jac_IP + end do + $:END_GPU_PARALLEL_LOOP() + end if + end subroutine s_update_igr + !> Subroutine that updates the conservative variables at the ghost points !! @param q_cons_vf Conservative Variables !! @param q_prim_vf Primitive variables @@ -196,18 +232,20 @@ contains type(ghost_point) :: gp type(ghost_point) :: innerp - ! set the Moving IBM interior Pressure Values - $:GPU_PARALLEL_LOOP(private='[i,j,k]', copyin='[E_idx]', collapse=3) - do l = 0, p - do k = 0, n - do j = 0, m - if (ib_markers%sf(j, k, l) /= 0) then - q_prim_vf(E_idx)%sf(j, k, l) = 1._wp - end if + if (.not. igr) then + ! set the Moving IBM interior Pressure Values + $:GPU_PARALLEL_LOOP(private='[i,j,k]', copyin='[E_idx]', collapse=3) + do l = 0, p + do k = 0, n + do j = 0, m + if (ib_markers%sf(j, k, l) /= 0) then + q_prim_vf(E_idx)%sf(j, k, l) = 1._wp + 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]') @@ -239,34 +277,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) 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 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 @@ -322,12 +378,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 @@ -340,6 +398,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) @@ -827,11 +886,14 @@ 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) $: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 @@ -847,6 +909,12 @@ contains integer :: i, j, k, l, q !< Iterator variables integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables real(wp) :: coeff + real(wp) :: alphaSum + 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 @@ -861,6 +929,7 @@ contains alpha_IP = 0._wp pres_IP = 0._wp vel_IP = 0._wp + pres = 0._wp if (surface_tension) c_IP = 0._wp @@ -887,31 +956,94 @@ 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 + alphaSum = 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 +! 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_vc) + 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 + + 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 + $: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) + alphaSum = alphaSum + 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 - alphaSum) + 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 @@ -926,7 +1058,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 @@ -940,9 +1072,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..71449bbdf9 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 + implicit none private; public :: s_initialize_igr_module, & @@ -368,7 +370,31 @@ 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_update_igr(jac_sf) + 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 085754abdc..be7e701047 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -732,8 +732,8 @@ 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(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, From 13404978e4e0e33fc0ae3561fe4eada0f68aaa9a Mon Sep 17 00:00:00 2001 From: Mark Zhang Date: Thu, 18 Dec 2025 19:24:33 -0500 Subject: [PATCH 2/5] Add comments --- src/simulation/m_ibm.fpp | 8 ++++++-- src/simulation/m_igr.fpp | 1 + toolchain/mfc/case_validator.py | 3 --- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index c8eb33ce7a..78ebb7190c 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -156,6 +156,8 @@ contains 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 @@ -195,7 +197,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), & @@ -311,7 +313,7 @@ contains end if ! set the pressure - ! !TEMPORARY, NEED TO FIX FOR IGR + ! !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 @@ -959,6 +961,8 @@ contains coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1) if (igr) then + ! For IGR, we will need to perform operations on + ! the conservative variables instead alphaSum = 0._wp dyn_pres = 0._wp if (num_fluids == 1) then diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp index 71449bbdf9..810e552d90 100644 --- a/src/simulation/m_igr.fpp +++ b/src/simulation/m_igr.fpp @@ -383,6 +383,7 @@ contains end do $:END_GPU_PARALLEL_LOOP() call s_update_igr(jac_sf) + ! 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 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index be7e701047..f7d99ad56e 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -709,7 +709,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' @@ -732,8 +731,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, From b07045b5ded352d6b65f6f9b444904e1fda7e437 Mon Sep 17 00:00:00 2001 From: Mark Zhang Date: Fri, 19 Dec 2025 00:11:43 -0500 Subject: [PATCH 3/5] Minor fixes --- src/simulation/m_ibm.fpp | 18 +++++++++++++----- src/simulation/m_igr.fpp | 2 +- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 78ebb7190c..80f79c91ac 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -150,6 +150,8 @@ contains end subroutine s_populate_ib_buffers + !> Interpolates sigma from the m_igr module at all ghost points + !! @param jac_sf Sigma, Entropic pressure subroutine s_update_igr(jac_sf) type(scalar_field), dimension(1), intent(inout) :: jac_sf integer :: j, k, l, r, s, t, i @@ -171,6 +173,11 @@ contains 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]') @@ -236,7 +243,7 @@ contains if (.not. igr) then ! set the Moving IBM interior Pressure Values - $:GPU_PARALLEL_LOOP(private='[i,j,k]', copyin='[E_idx]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[j,k,l]', copyin='[E_idx]', collapse=3) do l = 0, p do k = 0, n do j = 0, m @@ -911,7 +918,7 @@ contains integer :: i, j, k, l, q !< Iterator variables integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables real(wp) :: coeff - real(wp) :: alphaSum + 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 @@ -963,7 +970,7 @@ contains if (igr) then ! For IGR, we will need to perform operations on ! the conservative variables instead - alphaSum = 0._wp + 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) @@ -1015,14 +1022,15 @@ contains 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) - alphaSum = alphaSum + q_cons_vf(E_idx + l)%sf(i, j, k) + alpha_sum = alpha_sum + 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 - alphaSum) + alpha_IP(num_fluids) = alpha_IP(num_fluids) + coeff*(1._wp - alpha_sum) end if else pres_IP = pres_IP + coeff* & diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp index 810e552d90..d38ea8d2cd 100644 --- a/src/simulation/m_igr.fpp +++ b/src/simulation/m_igr.fpp @@ -15,7 +15,7 @@ module m_igr use m_boundary_common - use m_ibm + use m_ibm, only: ib, ib_markers, s_update_igr implicit none From 489bc7c53f29ec373889541545220f1ab79273c4 Mon Sep 17 00:00:00 2001 From: Mark Zhang Date: Fri, 19 Dec 2025 16:57:30 -0500 Subject: [PATCH 4/5] Change from using jac_sf to jac --- src/simulation/m_ibm.fpp | 12 ++++++------ src/simulation/m_igr.fpp | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 80f79c91ac..542a37a998 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -151,9 +151,9 @@ contains end subroutine s_populate_ib_buffers !> Interpolates sigma from the m_igr module at all ghost points - !! @param jac_sf Sigma, Entropic pressure - subroutine s_update_igr(jac_sf) - type(scalar_field), dimension(1), intent(inout) :: jac_sf + !! @param jac: Sigma, Entropic pressure + subroutine s_interpolate_sigma_igr(jac) + real(wp), dimension(:, :, :), intent(inout) :: jac integer :: j, k, l, r, s, t, i integer :: j1, j2, k1, k2, l1, l2 real(wp) :: coeff, jac_IP @@ -185,15 +185,15 @@ contains $: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_sf(1)%sf(j, k, l) + jac_IP = jac_IP + coeff*jac(j, k, l) end do end do end do - jac_sf(1)%sf(r, s, t) = jac_IP + jac(r, s, t) = jac_IP end do $:END_GPU_PARALLEL_LOOP() end if - end subroutine s_update_igr + end subroutine s_interpolate_sigma_igr !> Subroutine that updates the conservative variables at the ghost points !! @param q_cons_vf Conservative Variables diff --git a/src/simulation/m_igr.fpp b/src/simulation/m_igr.fpp index d38ea8d2cd..4f25ed9850 100644 --- a/src/simulation/m_igr.fpp +++ b/src/simulation/m_igr.fpp @@ -15,7 +15,7 @@ module m_igr use m_boundary_common - use m_ibm, only: ib, ib_markers, s_update_igr + use m_ibm, only: ib, ib_markers, s_interpolate_sigma_igr implicit none @@ -382,7 +382,7 @@ contains end do end do $:END_GPU_PARALLEL_LOOP() - call s_update_igr(jac_sf) + 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) From d914ba3030cca884621803056725577dafbd3e51 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Mon, 29 Dec 2025 13:03:37 -0500 Subject: [PATCH 5/5] add igr forward facing step example --- examples/2D_IGR_forward_facing_step/case.py | 97 +++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 examples/2D_IGR_forward_facing_step/case.py 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..65e47acc20 --- /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: No slip rectangle with rouded corner + "patch_ib(1)%geometry": 3, + "patch_ib(1)%x_centroid": 11.5 * h + rc*h, + "patch_ib(1)%y_centroid": 0 * h, + "patch_ib(1)%length_x": 17 * h, + "patch_ib(1)%length_y": 2 * 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": -rc * h, + "patch_ib(2)%length_x": rc*h, + "patch_ib(2)%length_y": 2 * 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, + ) +)