Merge branch 'master' into fp-stability-dd-confirm-rank

Author: sbryngelson

docs/documentation/case.md

@@ -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.

src/common/m_variables_conversion.fpp

@@ -1332,7 +1332,11 @@ contains
 #ifdef MFC_DEBUG
         if (disc < 0._wp) then
             print *, 'rho, c, Bx, By, Bz, h, term, disc:', rho, c, B(1), B(2), B(3), h, term, disc
+            ! s_mpi_abort is a host routine and cannot be called from device code
+            ! (this is a GPU routine); on GPU builds, emit the diagnostic print only.
+#ifndef MFC_GPU
             call s_mpi_abort('Error: negative discriminant in s_compute_fast_magnetosonic_speed')
+#endif
         end if
 #endif
 

src/post_process/m_derived_variables.fpp

@@ -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
@@ -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
+                            ! 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

tests/coverage_map.json.gz

@@ -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,
+                        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")