Skip to content

Fix numerical Schlieren for single fluid IGR + add case validation guard for schlieren_alpha #1477

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
Cowsreal wants to merge 5 commits into
base: master
Choose a base branch
from
+31 −6
Open

Fix numerical Schlieren for single fluid IGR + add case validation guard for schlieren_alpha #1477

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
8320e55
Fix numerical Schlieren for single fluid IGR and add guard for schlie…
Cowsreal Jun 1, 2026
7c27856
Add missing if (igr)
Cowsreal Jun 1, 2026
c6bfbe2
Minor adjustments
Cowsreal Jun 1, 2026
8d4a8bb
Formatting
Cowsreal Jun 1, 2026
76be6a6
Merge branch 'master' into minorChanges
Cowsreal Jun 1, 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
2 changes: 1 addition & 1 deletion docs/documentation/case.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -708,7 +708,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu

- ``[variable's name]_wrt`` activates the output of each specified variable into the database.

- `schlieren_alpha(i)` specifies the intensity of the numerical Schlieren of $i$-th component.
- `schlieren_alpha(i)` specifies the intensity of the numerical Schlieren of $i$-th component. It must be specified for every fluid when `schlieren_wrt` is enabled.

- `fd_order` specifies the order of the finite difference scheme used to compute the vorticity from the velocity field and the numerical schlieren from the density field using an integer of 1, 2, and 4.
`fd_order = 1`, `2`, and `4` correspond to the first, second, and fourth-order finite difference schemes.
Expand Down
23 changes: 19 additions & 4 deletions src/post_process/m_derived_variables.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -438,6 +438,7 @@ contains

real(wp) :: drho_dx, drho_dy, drho_dz !< Spatial derivatives of the density in the x-, y- and z-directions
real(wp), dimension(2) :: gm_rho_max !< Global (max gradient magnitude, rank) pair for density
real(wp) :: alpha_last !< Volume fraction of the fluid not explicitly stored (IGR)
integer :: i, j, k, l

do l = -offset_z%beg, p + offset_z%end
Expand Down Expand Up @@ -495,10 +496,24 @@ contains
do j = -offset_x%beg, m + offset_x%end
q_sf(j, k, l) = 0._wp

do i = 1, eqn_idx%adv%end - eqn_idx%E
q_sf(j, k, l) = q_sf(j, k, l) - schlieren_alpha(i)*q_cons_vf(i + eqn_idx%E)%sf(j, k, l)*gm_rho_sf(j, &
& k, l)/gm_rho_max(1)
end do
! Tracks the volume fraction of the fluid not explicitly stored (IGR reconstructs it as 1 - sum)
if (igr) then
Copy link
Copy Markdown
Member

@sbryngelson sbryngelson Jun 1, 2026

Choose a reason for hiding this comment

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

Optional (efficiency/clarity): both if (igr) checks (here and the one inside the i loop at line 504) test a loop-invariant flag evaluated on every (j,k,l) cell. Since igr is constant for the whole run, you could hoist the IGR/non-IGR split outside the triple loop — e.g. keep the inner alpha_last accumulation only when igr, or branch once above the loops. Negligible cost in post-processing, so not blocking, just a readability note.

Also note this new term divides by gm_rho_max(1), which is 0 for a perfectly uniform density field — the same pre-existing risk as line 502, so no regression, just flagging.

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jun 1, 2026

Choose a reason for hiding this comment

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

Will fix.

! IGR stores only num_fluids-1 volume fractions; the last fluid's volume fraction is untracked.
! For a single fluid this is simply 1.0 everywhere. Without this term, the entire single-fluid
! Schlieren field is dropped, leaving exp(0) = 1 everywhere. Compute that term below.
alpha_last = 1._wp
do i = 1, eqn_idx%adv%end - eqn_idx%E
q_sf(j, k, l) = q_sf(j, k, l) - schlieren_alpha(i)*q_cons_vf(i + eqn_idx%E)%sf(j, k, &
& l)*gm_rho_sf(j, k, l)/gm_rho_max(1)
alpha_last = alpha_last - q_cons_vf(i + eqn_idx%E)%sf(j, k, l)
end do
q_sf(j, k, l) = q_sf(j, k, l) - schlieren_alpha(num_fluids)*alpha_last*gm_rho_sf(j, k, l)/gm_rho_max(1)
else
do i = 1, eqn_idx%adv%end - eqn_idx%E
q_sf(j, k, l) = q_sf(j, k, l) - schlieren_alpha(i)*q_cons_vf(i + eqn_idx%E)%sf(j, k, &
& l)*gm_rho_sf(j, k, l)/gm_rho_max(1)
end do
end if
end do
end do
end do
Expand Down
12 changes: 11 additions & 1 deletion toolchain/mfc/case_validator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1558,14 +1558,24 @@ def check_schlieren(self):
n = self.get("n", 0)
fd_order = self.get("fd_order")
num_fluids = self.get("num_fluids")
model_eqns = self.get("model_eqns")

self.prohibit(n is not None and n == 0 and schlieren_wrt, "schlieren_wrt requires n > 0 (at least 2D)")
self.prohibit(schlieren_wrt and fd_order is None, "fd_order must be set for schlieren_wrt")

# The volume-fraction model (model_eqns /= 1) weights the Schlieren field by schlieren_alpha(i) for every fluid; an unset
# value stays at the -1e6 sentinel and produces nonsensical output. Under IGR the last fluid's volume fraction is
# reconstructed (not stored), so schlieren_alpha must still be set for all num_fluids components, including the single
# fluid of a single-fluid IGR case. The gamma/pi_inf model (model_eqns == 1) does not use schlieren_alpha.
if num_fluids is not None:
for i in range(1, num_fluids + 1):
schlieren_alpha = self.get(f"schlieren_alpha({i})")
if schlieren_alpha is not None:
if schlieren_alpha is None:
self.prohibit(
schlieren_wrt and model_eqns is not None and model_eqns != 1,
Copy link
Copy Markdown
Member

@sbryngelson sbryngelson Jun 1, 2026

Choose a reason for hiding this comment

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

The guard correctly excludes model_eqns == 1, which matches the Fortran (the model_eqns == 1 branch never uses schlieren_alpha). One edge case to consider: IGR's adv-index layout is only set up under model_eqns == 2 (m_global_parameters.fpp), but nothing forbids igr = T with model_eqns ∈ {3,4}. In that (unsupported) combo the new if (igr) term in m_derived_variables.fpp would double-count the last fluid. IGR is fundamentally a model_eqns == 2 feature and would likely fail elsewhere, so this is theoretical — but a belt-and-suspenders self.prohibit(igr and model_eqns is not None and model_eqns != 2, ...) in check_igr would make the assumption explicit.

Minor: if model_eqns is unset (None), the guard won't fire even with schlieren_wrt on; fine in practice since model_eqns is virtually always specified.

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jun 1, 2026

Choose a reason for hiding this comment

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

Irrelevant, IGR is blocked to run for model_eqn /= 2

f"schlieren_alpha({i}) must be set for every fluid when schlieren_wrt is enabled (unless model_eqns == 1)",
)
else:
self.prohibit(schlieren_alpha <= 0, f"schlieren_alpha({i}) must be greater than zero")
self.prohibit(not schlieren_wrt, f"schlieren_alpha({i}) should be set only with schlieren_wrt enabled")

Expand Down
Loading