Clamp the cubic Hugoniot at its turnover, and count the relaxation solves that miss

Author: sbryngelson

src/common/m_global_parameters_common.fpp

@@ -59,11 +59,12 @@ module m_global_parameters_common
     !> Reference state and Gruneisen closure Gamma_G = Gamma_0 + a mu of a state-dependent EOS, whatever its family
     real(wp), allocatable, dimension(:) :: rho0s, t0s, gruneisen0s, gruneisen_as
     real(wp), allocatable, dimension(:) :: mg_c0s, mg_ss, mg_s2s, mg_s3s
+    real(wp), allocatable, dimension(:) :: mg_mu_maxs  !< Compression at which a cubic Hugoniot fit turns over
     real(wp), allocatable, dimension(:) :: jwl_as, jwl_bs, jwl_r1s, jwl_r2s
     real(wp), allocatable, dimension(:) :: vinet_k0s, vinet_k0ps
-    logical :: any_state_dependent_eos  !< True when some fluid's coefficients vary with density; set at init
-    $:GPU_DECLARE(create='[eoss, rho0s, t0s, gruneisen0s, gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, jwl_as, jwl_bs, jwl_r1s, &
-                  & jwl_r2s, vinet_k0s, vinet_k0ps, any_state_dependent_eos]')
+    logical :: any_state_dependent_eos                 !< True when some fluid's coefficients vary with density; set at init
+    $:GPU_DECLARE(create='[eoss, rho0s, t0s, gruneisen0s, gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, mg_mu_maxs, jwl_as, &
+                  & jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, vinet_k0ps, any_state_dependent_eos]')
     !> @}
 
     !> @name Fluids participating in shear and bulk viscosity

src/common/m_variables_conversion.fpp

@@ -270,8 +270,8 @@ contains
         @:ALLOCATE(gammas (1:num_fluids))
         @:ALLOCATE(eoss (1:num_fluids), rho0s (1:num_fluids), t0s (1:num_fluids), gruneisen0s (1:num_fluids), &
                    & gruneisen_as (1:num_fluids), mg_c0s (1:num_fluids), mg_ss (1:num_fluids), mg_s2s (1:num_fluids), &
-                   & mg_s3s (1:num_fluids), jwl_as (1:num_fluids), jwl_bs (1:num_fluids), jwl_r1s (1:num_fluids), &
-                   & jwl_r2s (1:num_fluids), vinet_k0s (1:num_fluids), vinet_k0ps (1:num_fluids))
+                   & mg_s3s (1:num_fluids), mg_mu_maxs (1:num_fluids), jwl_as (1:num_fluids), jwl_bs (1:num_fluids), &
+                   & jwl_r1s (1:num_fluids), jwl_r2s (1:num_fluids), vinet_k0s (1:num_fluids), vinet_k0ps (1:num_fluids))
         @:ALLOCATE(isentrope_n (1:num_fluids))
         @:ALLOCATE(pi_infs(1:num_fluids))
         @:ALLOCATE(isentrope_B(1:num_fluids))
@@ -303,6 +303,9 @@ contains
             mg_ss(i) = fluid_pp(i)%mg_s
             mg_s2s(i) = fluid_pp(i)%mg_s2
             mg_s3s(i) = fluid_pp(i)%mg_s3
+            ! Where a cubic Hugoniot fit turns over: mu(u_p) peaks where c0 = s2 u_p^2 + 2 s3 u_p^3, and past it
+            ! no shock state exists, so the Newton below would wander. Solved once here, on the host.
+            mg_mu_maxs(i) = f_hugoniot_compression_limit(fluid_pp(i)%mg_c0, fluid_pp(i)%mg_s, fluid_pp(i)%mg_s2, fluid_pp(i)%mg_s3)
             jwl_as(i) = fluid_pp(i)%jwl_a
             jwl_bs(i) = fluid_pp(i)%jwl_b
             jwl_r1s(i) = fluid_pp(i)%jwl_r1
@@ -335,7 +338,8 @@ contains
             if (f_is_state_dependent(i)) any_state_dependent_eos = .true.
         end do
         $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, rho0s, t0s, gruneisen0s, &
-                     & gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, jwl_as, jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, vinet_k0ps, any_state_dependent_eos]')
+                     & gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, mg_mu_maxs, jwl_as, jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, &
+                     & vinet_k0ps, any_state_dependent_eos]')
 
         @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max)))
         Res_vc = dflt_real
@@ -1219,7 +1223,7 @@ contains
         if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf)
 
         @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, rho0s, t0s, gruneisen0s, &
-                     & gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, jwl_as, jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, vinet_k0ps)
+                     & gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, mg_mu_maxs, jwl_as, jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, vinet_k0ps)
         if (allocated(bubrs_vc)) then
             @:DEALLOCATE(bubrs_vc)
         end if
@@ -1353,6 +1357,9 @@ contains
         integer               :: iter
 
         mu = rho/rho0s(i) - 1._wp
+        ! Past the fit's turnover there is no shock state to find; clamp rather than let the Newton below wander
+        ! off and return a silently wrong pressure. mg_mu_maxs is huge for the linear fit, so this is a no-op there.
+        if (eoss(i) == eos_mie_gruneisen .and. mu > mg_mu_maxs(i)) mu = mg_mu_maxs(i)
         select case (eoss(i))
         case (eos_mie_gruneisen)
             ! Hugoniot reference u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3, with p_H = rho0 u_s u_p and the Hugoniot
@@ -1435,6 +1442,33 @@ contains
 
     end function f_has_isentropic_reference
 
+    !> The largest compression a cubic Hugoniot fit can represent. mu(u_p) = u_p/(u_s - u_p) rises, peaks where c0 = s2 u_p^2 + 2 s3
+    !! u_p^3, and falls after; only the rising branch is a physical shock. Returns a huge value for the linear fit, which never
+    !! turns over. Host-side: called once per fluid at initialization.
+    impure function f_hugoniot_compression_limit(c0, s, s2, s3) result(mu_max)
+
+        real(wp), intent(in) :: c0, s, s2, s3
+        real(wp)             :: mu_max, up, f, df, us
+        integer              :: iter
+
+        if (s2 == 0._wp .and. s3 == 0._wp) then
+            mu_max = huge(1._wp)
+            return
+        end if
+
+        ! Newton on c0 - s2 u^2 - 2 s3 u^3 = 0, from a guess that brackets the physical range
+        up = c0
+        do iter = 1, 100
+            f = c0 - s2*up**2 - 2._wp*s3*up**3
+            df = -2._wp*s2*up - 6._wp*s3*up**2
+            if (abs(df) < verysmall) exit
+            up = max(up - f/df, verysmall)
+        end do
+        us = c0 + up*(s + up*(s2 + up*s3))
+        mu_max = up/max(us - up, verysmall)
+
+    end function f_hugoniot_compression_limit
+
     !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and
     !! ideal gas keep the constants resolved at init, bit for bit.
     subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)

src/simulation/m_pressure_relaxation.fpp

@@ -16,7 +16,12 @@ module m_pressure_relaxation
 
     implicit none
 
-    private; public :: s_pressure_relaxation_procedure
+    private; public :: s_pressure_relaxation_procedure, s_report_pressure_relaxation
+
+    !> Cell updates in which the Newton below stopped on its iteration cap rather than on the tolerance. Counted because the answer
+    !! then depends on that cap, and nothing else in the code says so.
+    integer  :: n_hit_cap = 0
+    real(wp) :: worst_residual = 0._wp  !< Largest |f| left behind when the cap was reached
 
 contains
 
@@ -32,18 +37,27 @@ contains
             real(wp), dimension(num_fluids) :: alpha_rho, alpha
         #:endif
         real(wp) :: rho, gamma, pi_inf, qv_mix
+        integer  :: hit_cap, hit_cap_sum
+        real(wp) :: resid, resid_max
 
         ! Formed here, not one call deeper: CCE OpenACC accepts a num_fluids-sized array passed to a device routine from a
         ! parallel-loop body, and rejects the same call from inside another acc routine seq.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, l, alpha_rho, alpha, rho, gamma, pi_inf, qv_mix]', collapse=3)
+        hit_cap_sum = 0
+        resid_max = 0._wp
+        $:GPU_PARALLEL_LOOP(private='[i, j, k, l, alpha_rho, alpha, rho, gamma, pi_inf, qv_mix, hit_cap, resid]', collapse=3, &
+                            & reduction='[[hit_cap_sum], [resid_max]]', reductionOp='[+, max]')
         do l = 0, p
             do k = 0, n
                 do j = 0, m
                     if (mpp_lim) call s_correct_volume_fractions(q_cons_vf, j, k, l)
 
+                    hit_cap = 0
+                    resid = 0._wp
                     if (s_needs_pressure_relaxation(q_cons_vf, j, k, l)) then
-                        call s_equilibrate_pressure(q_cons_vf, j, k, l)
+                        call s_equilibrate_pressure(q_cons_vf, j, k, l, hit_cap, resid)
                     end if
+                    hit_cap_sum = hit_cap_sum + hit_cap
+                    resid_max = max(resid_max, resid)
 
                     $:GPU_LOOP(parallelism='[seq]')
                     do i = 1, num_fluids
@@ -58,9 +72,22 @@ contains
             end do
         end do
         $:END_GPU_PARALLEL_LOOP()
+        n_hit_cap = n_hit_cap + hit_cap_sum
+        worst_residual = max(worst_residual, resid_max)
 
     end subroutine s_pressure_relaxation_procedure
 
+    !> One line at the end of a run if the equilibration ever stopped on its iteration cap. Silence means every cell reached the
+    !! tolerance.
+    impure subroutine s_report_pressure_relaxation
+
+        if (n_hit_cap > 0 .and. proc_rank == 0) then
+            print '(A,I0,A,ES10.3,A)', ' Pressure relaxation reached its iteration cap in ', n_hit_cap, &
+                & ' cell updates; worst residual left behind ', worst_residual, ' against a 1e-10 tolerance.'
+        end if
+
+    end subroutine s_report_pressure_relaxation
+
     !> Check if pressure relaxation is needed for this cell
     logical function s_needs_pressure_relaxation(q_cons_vf, j, k, l)
 
@@ -111,12 +138,14 @@ contains
     end subroutine s_correct_volume_fractions
 
     !> Main pressure equilibration using Newton-Raphson
-    subroutine s_equilibrate_pressure(q_cons_vf, j, k, l)
+    subroutine s_equilibrate_pressure(q_cons_vf, j, k, l, hit_cap, resid)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         integer, intent(in)                                    :: j, k, l
+        integer, intent(out)                                   :: hit_cap
+        real(wp), intent(out)                                  :: resid
         real(wp)                                               :: pres_relax, f_pres, df_pres
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
             real(wp), dimension(3) :: pres_K_init, rho_K_init, rho_K_s
@@ -194,6 +223,17 @@ contains
             end if
         end do
 
+        ! Written as .not. (<=) so a NaN residual is reported as a miss, like a diverged one.
+        hit_cap = 0
+        resid = 0._wp
+        if (.not. (abs(f_pres) <= TOLERANCE)) then
+            hit_cap = 1
+            ! A NaN residual has to be mapped, not passed on: max() with a NaN keeps the other operand, so the
+            ! reduction would report the miss as zero.
+            resid = abs(f_pres)
+            if (f_pres /= f_pres) resid = huge(1._wp)
+        end if
+
         ! Update volume fractions. The Newton above often stops on the iteration cap rather than on the
         ! tolerance, and the answer then depends on that cap, so an unconverged-but-physical density is still
         ! used -- the alternative would move every six-equation result. What is refused is a density that is

src/simulation/m_start_up.fpp

@@ -22,6 +22,7 @@ module m_start_up
     use m_boundary_io
     use m_acoustic_src
     use m_rhs
+    use m_pressure_relaxation, only: s_report_pressure_relaxation
     use m_chemistry
     use m_data_output
     use m_time_steppers
@@ -1108,6 +1109,8 @@ contains
     !> Finalize and deallocate all simulation sub-modules in reverse initialization order
     impure subroutine s_finalize_modules
 
+        if (model_eqns == model_eqns_6eq) call s_report_pressure_relaxation()
+
         call s_finalize_time_steppers_module()
         if (hypoelasticity) call s_finalize_hypoelastic_module()
         call s_finalize_derived_variables_module()