Skip to content

refactor: unify equation/range indices into sys_idx; clean BC globals #1365

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
sbryngelson wants to merge 14 commits into
base: master
Choose a base branch
from
Open

refactor: unify equation/range indices into sys_idx; clean BC globals #1365

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
10ca810
refactor: replace implicit bc globals with bc_xyz_info derived type; …
sbryngelson Apr 9, 2026
db99a8a
refactor: consolidate scalar equation indices into eqn_idx derived type
sbryngelson Apr 9, 2026
28b88d0
fix: replace Fypp bc${XYZ}$b/e templates with bc_${XYZ}$%beg/end; fix…
sbryngelson Apr 9, 2026
44592ab
refactor: consolidate all equation/range indices into sys_idx derived…
sbryngelson Apr 10, 2026
a19cae2
refactor: rewrite indices.dat summary using sys_idx; add missing xi, …
sbryngelson Apr 10, 2026
4022b72
refactor: split qbmm_idx out of bub_bounds_info; sys_idx is now alloc…
sbryngelson Apr 10, 2026
141226f
refactor: rename sys_idx -> eqn_idx, sys_idx_info -> eqn_idx_info
sbryngelson Apr 10, 2026
a4c18bd
refactor: replace fragile sys_size bounds with explicit eqn_idx refer…
sbryngelson Apr 10, 2026
6e7e148
fix: address PR review findings from eqn-idx-refactor
sbryngelson Apr 10, 2026
6737414
fix: use then/end if with @:DEALLOCATE macro for qbmm_idx%moms
sbryngelson Apr 10, 2026
5bb8b8d
fix: update QPVF_IDX_VARS in case.py to use eqn_idx references
sbryngelson Apr 10, 2026
8c93ebc
refactor: replace shorthand contxb/momxb/advxb/chemxb with eqn_idx re…
sbryngelson Apr 10, 2026
9d45f6c
refactor: replace shorthand index variables with eqn_idx references t…
sbryngelson Apr 10, 2026
bb715b1
refactor: replace xibeg/xiend shorthands with eqn_idx%xi%beg/end refe…
sbryngelson Apr 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions src/common/include/1dHardcodedIC.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,12 @@
select case (patch_icpp(patch_id)%hcid)
case (150) ! 1D Smooth Alfven Case for MHD
! velocity
q_prim_vf(momxb + 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x_cc(i))
q_prim_vf(momxb + 2)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x_cc(i))
q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x_cc(i))
q_prim_vf(eqn_idx%mom%beg + 2)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x_cc(i))

! magnetic field
q_prim_vf(B_idx%end - 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x_cc(i))
q_prim_vf(B_idx%end)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x_cc(i))
q_prim_vf(eqn_idx%B%end - 1)%sf(i, 0, 0) = 0.1_wp*sin(2._wp*pi*x_cc(i))
q_prim_vf(eqn_idx%B%end)%sf(i, 0, 0) = 0.1_wp*cos(2._wp*pi*x_cc(i))
case (170) ! 1D profile from external data (e.g. Cantera, SDtoolbox)
! This hardcoded case can be used to start a simulation with initial conditions given from a known 1D profile (e.g. Cantera,
! SDtoolbox)
Expand All @@ -21,36 +21,36 @@
! 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))
q_prim_vf(eqn_idx%cont%beg + 0)%sf(i, 0, 0) = 1 + 0.2*sin(5*x_cc(i))
end if
case (181) ! Titarev-Torro problem
! 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)*pi)
q_prim_vf(eqn_idx%cont%beg + 0)%sf(i, 0, 0) = 1 + 0.1*sin(20*x_cc(i)*pi)
case (182) ! Multi-component diffusion
! This patch is a hard-coded for test suite optimization (multiple component diffusion)
x_mid_diffu = 0.05_wp/2.0_wp
width_sq = (2.5_wp*10.0_wp**(-3.0_wp))**2
profile_shape = 1.0_wp - 0.5_wp*exp(-(x_cc(i) - x_mid_diffu)**2/width_sq)
q_prim_vf(momxb)%sf(i, 0, 0) = 0.0_wp
q_prim_vf(E_idx)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5
q_prim_vf(advxb)%sf(i, 0, 0) = 1.0_wp
q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5
q_prim_vf(eqn_idx%adv%beg)%sf(i, 0, 0) = 1.0_wp

y1 = (0.195_wp - 0.142_wp)*profile_shape + 0.142_wp
y2 = (0.0_wp - 0.1_wp)*profile_shape + 0.1_wp
y3 = (0.214_wp - 0.0_wp)*profile_shape + 0.0_wp
y4 = (0.591_wp - 0.758_wp)*profile_shape + 0.758_wp

q_prim_vf(chemxb)%sf(i, 0, 0) = y1
q_prim_vf(chemxb + 1)%sf(i, 0, 0) = y2
q_prim_vf(chemxb + 2)%sf(i, 0, 0) = y3
q_prim_vf(chemxb + 3)%sf(i, 0, 0) = y4
q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = y1
q_prim_vf(eqn_idx%species%beg + 1)%sf(i, 0, 0) = y2
q_prim_vf(eqn_idx%species%beg + 2)%sf(i, 0, 0) = y3
q_prim_vf(eqn_idx%species%beg + 3)%sf(i, 0, 0) = y4

temp = (320.0_wp - 1350.0_wp)*profile_shape + 1350.0_wp

molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp

q_prim_vf(contxb)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
case default
call s_int_to_str(patch_id, iStr)
call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))
Expand Down
195 changes: 98 additions & 97 deletions src/common/include/2dHardcodedIC.fpp
View file
Open in desktop

Large diffs are not rendered by default.

68 changes: 34 additions & 34 deletions src/common/include/3dHardcodedIC.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,30 +84,30 @@
if (alph > 1._wp - eps) alph = 1._wp - eps

if (y_cc(j) > intH) then
q_prim_vf(advxb)%sf(i, j, k) = alph
q_prim_vf(advxe)%sf(i, j, k) = 1._wp - alph
q_prim_vf(contxb)%sf(i, j, k) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, k) = (1._wp - alph)*rhoL
q_prim_vf(E_idx)%sf(i, j, k) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j))
q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = alph
q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = alph*rhoH
q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = (1._wp - alph)*rhoL
q_prim_vf(eqn_idx%E)%sf(i, j, k) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j))
else
q_prim_vf(advxb)%sf(i, j, k) = alph
q_prim_vf(advxe)%sf(i, j, k) = 1._wp - alph
q_prim_vf(contxb)%sf(i, j, k) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, k) = (1._wp - alph)*rhoL
q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = alph
q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = alph*rhoH
q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = (1._wp - alph)*rhoL
pInt = pref + rhoH*9.81_wp*(1.2_wp - intH)
q_prim_vf(E_idx)%sf(i, j, k) = pInt + rhoL*9.81_wp*(intH - y_cc(j))
q_prim_vf(eqn_idx%E)%sf(i, j, k) = pInt + rhoL*9.81_wp*(intH - y_cc(j))
end if
case (301) ! (3D lung geometry in X direction, |sin(*)+sin(*)|)
h = 0.0_wp
lam = 1.0_wp
amp = patch_icpp(patch_id)%a(2)
intH = amp*abs((sin(2*pi*y_cc(j)/lam - pi/2) + sin(2*pi*z_cc(k)/lam - pi/2)) + h)
if (x_cc(i) > intH) then
q_prim_vf(contxb)%sf(i, j, k) = patch_icpp(1)%alpha_rho(1)
q_prim_vf(contxe)%sf(i, j, k) = patch_icpp(1)%alpha_rho(2)
q_prim_vf(E_idx)%sf(i, j, k) = patch_icpp(1)%pres
q_prim_vf(advxb)%sf(i, j, k) = patch_icpp(1)%alpha(1)
q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2)
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = patch_icpp(1)%alpha_rho(1)
q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = patch_icpp(1)%alpha_rho(2)
q_prim_vf(eqn_idx%E)%sf(i, j, k) = patch_icpp(1)%pres
q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = patch_icpp(1)%alpha(1)
q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = patch_icpp(1)%alpha(2)
end if
case (302) ! 3D Jet with IGR
ux_th = 10*sqrt(1.4*0.4)
Expand All @@ -126,19 +126,19 @@
rcut = f_cut_on(r - r_th, eps_smooth)
xcut = f_cut_on(x_cc(i), eps_smooth)

q_prim_vf(momxb)%sf(i, j, k) = ux_th*rcut*xcut + ux_am
q_prim_vf(momxb + 1)%sf(i, j, k) = 0._wp
q_prim_vf(momxe)%sf(i, j, k) = 0._wp
q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) = ux_th*rcut*xcut + ux_am
q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, k) = 0._wp
q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = 0._wp

if (num_fluids == 1) then
q_prim_vf(contxb)%sf(i, j, k) = (rho_th - rho_am)*rcut*xcut + rho_am
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = (rho_th - rho_am)*rcut*xcut + rho_am
else
q_prim_vf(advxb)%sf(i, j, k) = (1._wp - 2._wp*eps)*rcut*xcut + eps
q_prim_vf(contxb)%sf(i, j, k) = rho_th*q_prim_vf(advxb)%sf(i, j, k)
q_prim_vf(contxe)%sf(i, j, k) = rho_am*(1._wp - q_prim_vf(advxb)%sf(i, j, k))
q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = (1._wp - 2._wp*eps)*rcut*xcut + eps
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = rho_th*q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k)
q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = rho_am*(1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))
end if

q_prim_vf(E_idx)%sf(i, j, k) = p_th*rcut*xcut + p_am
q_prim_vf(eqn_idx%E)%sf(i, j, k) = p_th*rcut*xcut + p_am
case (303) ! 3D Multijet
eps_smooth = 3.0_wp
ux_th = 10*sqrt(1.4*0.4)
Expand All @@ -152,19 +152,19 @@
rcut = rcut_arr(j, k)
xcut = f_cut_on(x_cc(i), eps_smooth)

q_prim_vf(momxb)%sf(i, j, k) = ux_th*rcut*xcut + ux_am
q_prim_vf(momxb + 1)%sf(i, j, k) = 0._wp
q_prim_vf(momxe)%sf(i, j, k) = 0._wp
q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) = ux_th*rcut*xcut + ux_am
q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, k) = 0._wp
q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = 0._wp

if (num_fluids == 1) then
q_prim_vf(contxb)%sf(i, j, k) = (rho_th - rho_am)*rcut*xcut + rho_am
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = (rho_th - rho_am)*rcut*xcut + rho_am
else
q_prim_vf(advxb)%sf(i, j, k) = (1._wp - 2._wp*eps)*rcut*xcut + eps
q_prim_vf(contxb)%sf(i, j, k) = rho_th*q_prim_vf(advxb)%sf(i, j, k)
q_prim_vf(contxe)%sf(i, j, k) = rho_am*(1._wp - q_prim_vf(advxb)%sf(i, j, k))
q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = (1._wp - 2._wp*eps)*rcut*xcut + eps
q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = rho_th*q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k)
q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = rho_am*(1._wp - q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k))
end if

q_prim_vf(E_idx)%sf(i, j, k) = p_th*rcut*xcut + p_am
q_prim_vf(eqn_idx%E)%sf(i, j, k) = p_th*rcut*xcut + p_am
case (370) ! 3D extrusion of 2D profile from external data
! This hardcoded case extrudes a 2D profile to initialize a 3D simulation domain
@: HardcodedReadValues()
Expand All @@ -173,10 +173,10 @@
! geometry 9
Mach = 0.1
if (patch_id == 1) then
q_prim_vf(E_idx)%sf(i, j, &
q_prim_vf(eqn_idx%E)%sf(i, j, &
& k) = 101325 + (Mach**2*376.636429464809**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) = Mach*376.636429464809*sin(x_cc(i)/1)*cos(y_cc(j)/1)*sin(z_cc(k)/1)
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)
q_prim_vf(eqn_idx%mom%beg + 0)%sf(i, j, k) = Mach*376.636429464809*sin(x_cc(i)/1)*cos(y_cc(j)/1)*sin(z_cc(k)/1)
q_prim_vf(eqn_idx%mom%beg + 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 default
call s_int_to_str(patch_id, iStr)
Expand Down
12 changes: 6 additions & 6 deletions src/common/include/ExtrusionHardcodedIC.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@
!>
!> **Data Assignment:**
!> - Populates q_prim_vf primitive variable arrays with file data
!> - Handles momentum component indexing with special treatment for momxe
!> - Sets momxe component to zero for 2D/3D cases
!> - Handles momentum component indexing with special treatment for eqn_idx%mom%end
!> - Sets eqn_idx%mom%end component to zero for 2D/3D cases
!>
!> **State Management:**
!> - Uses files_loaded flag to prevent redundant file operations
Expand Down Expand Up @@ -177,18 +177,18 @@
case (2)
idx = i + 1 + global_offset_x - index_x
do f = 1, sys_size - 1
jump = merge(1, 0, f >= momxe)
jump = merge(1, 0, f >= eqn_idx%mom%end)
q_prim_vf(f + jump)%sf(i, j, 0) = stored_values(idx, 1, f)
end do
q_prim_vf(momxe)%sf(i, j, 0) = 0.0_wp
q_prim_vf(eqn_idx%mom%end)%sf(i, j, 0) = 0.0_wp
case (3)
idx = i + 1 + global_offset_x - index_x
idy = j + 1 + global_offset_y - index_y
do f = 1, sys_size - 1
jump = merge(1, 0, f >= momxe)
jump = merge(1, 0, f >= eqn_idx%mom%end)
q_prim_vf(f + jump)%sf(i, j, k) = stored_values(idx, idy, f)
end do
q_prim_vf(momxe)%sf(i, j, k) = 0.0_wp
q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = 0.0_wp
end select
#:enddef

Expand Down
Loading
Loading