Bound chemistry face states with a scaling limiter (NaN at finite-cylinder rims)

src/simulation/m_rhs.fpp

@@ -16,7 +16,7 @@ module m_rhs
     use m_eos
     use m_weno
     use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_6eq, int_comp_mthinc, recon_type_weno, &
-        & recon_type_muscl
+        & recon_type_muscl, T_surface_max, verysmall, sgm_eps
     use m_muscl
     use m_riemann_solvers
     use m_cbc
@@ -33,6 +33,7 @@ module m_rhs
     use m_surface_tension
     use m_body_forces
     use m_chemistry
+    use m_thermochem, only: get_specific_gas_constant
     use m_conduction
     use m_reactive_burn
     use m_igr
@@ -997,6 +998,8 @@ contains
             end if
         end if
 
+        if (chemistry) call s_bound_face_states(id)
+
         ! Reconstruct viscous derivatives for viscosity
         if (weno_Re_flux) then
             iv%beg = eqn_idx%mom%beg; iv%end = eqn_idx%mom%end
@@ -1018,6 +1021,91 @@ contains
 
     end subroutine s_reconstruct_riemann_states
 
+    !> Scales each cell's chemistry face states toward its cell state by one theta (Zhang & Shu, JCP 2010) so that face rho and p
+    !! stay positive and face T <= T_surface_max, past which cp can go negative and the sound speed NaN. WENO overshoots there at
+    !! sharp IB corners, whose hot ghost layers mirror different faces (#1964).
+    subroutine s_bound_face_states(id)
+
+        integer, intent(in)                  :: id
+        type(int_bounds_info), dimension(3)  :: b
+        real(wp), dimension(${NUM_SPECIES}$) :: Ys
+        real(wp), dimension(3)               :: vb, vL, vR  !< (rho, p, R_mix) of the cell and its two faces
+        real(wp)                             :: theta
+        integer                              :: i, j, k, l, polyn
+
+        ! Same range as s_reconstruct_cell_boundary_values
+
+        polyn = merge(weno_polyn, muscl_polyn, recon_type == recon_type_weno)
+        b = idwbuff
+        b(id)%beg = b(id)%beg + polyn; b(id)%end = b(id)%end - polyn
+
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, Ys, vb, vL, vR, theta]', copyin='[b]')
+        do l = b(3)%beg, b(3)%end
+            do k = b(2)%beg, b(2)%end
+                do j = b(1)%beg, b(1)%end
+                    #:for V, Q in [('vb', 'q_prim_qp%vf({})%sf(j, k, l)'), ('vL', 'qL_rsx_vf(j, k, l, {})'), ('vR', &
+                                    & 'qR_rsx_vf(j, k, l, {})')]
+                        $:GPU_LOOP(parallelism='[seq]')
+                        do i = 1, num_species
+                            Ys(i) = ${Q.format('eqn_idx%species%beg + i - 1')}$
+                        end do
+                        call get_specific_gas_constant(Ys, ${V}$(3))
+                        ${V}$(1) = ${Q.format('eqn_idx%cont%beg')}$; ${V}$(2) = ${Q.format('eqn_idx%E')}$
+                    #:endfor
+                    ! rho, p and R_mix are linear in theta, so positivity is closed form and T <= T_max quadratic
+                    theta = 1._wp
+                    $:GPU_LOOP(parallelism='[seq]')
+                    do i = 1, 3
+                        theta = min(theta, f_theta_floor(vb(i), vL(i)), f_theta_floor(vb(i), vR(i)))
+                    end do
+                    theta = min(f_theta_T(theta, vb, vL), f_theta_T(theta, vb, vR))
+                    if (theta < 1._wp) then
+                        $:GPU_LOOP(parallelism='[seq]')
+                        do i = 1, sys_size
+                            qL_rsx_vf(j, k, l, i) = q_prim_qp%vf(i)%sf(j, k, l) + theta*(qL_rsx_vf(j, k, l, &
+                                      & i) - q_prim_qp%vf(i)%sf(j, k, l))
+                            qR_rsx_vf(j, k, l, i) = q_prim_qp%vf(i)%sf(j, k, l) + theta*(qR_rsx_vf(j, k, l, &
+                                      & i) - q_prim_qp%vf(i)%sf(j, k, l))
+                        end do
+                    end if
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+    end subroutine s_bound_face_states
+
+    !> Largest theta in [0, 1] keeping vbar + theta*(v - vbar) >= verysmall*vbar (0 if vbar <= 0).
+    pure function f_theta_floor(vbar, v) result(theta)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in) :: vbar, v
+        real(wp)             :: theta
+
+        theta = 1._wp
+        if (v < verysmall*vbar) theta = max(0._wp, (1._wp - verysmall)*vbar)/max(vbar - v, sgm_eps)
+
+    end function f_theta_floor
+
+    !> Largest theta' <= theta keeping the face (rho, p, R_mix) scaled toward vb by theta' at T <= T_surface_max (0 if vb is above
+    !! it). g = T_max rho R_mix - p is quadratic in theta'; where g(theta) < 0 its first root is the stable smaller-root form.
+    pure function f_theta_T(theta, vb, vf) result(t)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in)               :: theta
+        real(wp), dimension(3), intent(in) :: vb, vf
+        real(wp), dimension(3)             :: d
+        real(wp)                           :: a, b, c, t
+
+        d = vf - vb
+        a = T_surface_max*d(1)*d(3)
+        b = T_surface_max*(vb(1)*d(3) + vb(3)*d(1)) - d(2)
+        c = T_surface_max*vb(1)*vb(3) - vb(2)
+        t = theta
+        if ((a*theta + b)*theta + c < 0._wp) t = max(0._wp, 2._wp*c/max(sqrt(max(b*b - 4._wp*a*c, 0._wp)) - b, sgm_eps))
+
+    end function f_theta_T
+
     !> Computes one sweep direction's contribution to the RHS: the Riemann solve on the reconstructed cell-boundary states followed
     !! by the advection source term. For dual-pass HLLD the fused solve computes BOTH anchored flux sets in one call; this routine
     !! assembles the hat_L (is_hat_L=.true.) partial RHS, and the caller finalizes + assembles the hat_R set afterwards. is_hat_L