evaluate_step_rejection – LowMachReact-Hex

subroutine evaluate_step_rejection(params, stats, max_dT_chem, max_dY_chem, max_rel_drho_chem, ke_before, ke_after, cfl_after, fields, energy, transport, reject, reason)

Uses

- ieee_arithmetic

Evaluate multi-guard physical and numerical safety gates for step rejection and retry.

Evaluates multiple safety constraints over the trial timestep. If any limits are violated (e.g. NaN/Inf presence, chemistry increments, density vacuum floors, temperature excursions, or PCG linear solver convergence bounds), the trial step is rejected, triggering a temporal step reduction and retry.

Arguments

Type	Intent		Name
type(case_params_t),	intent(in)	::	params	Input case parameters and user namelist settings.
type(solver_stats_t),	intent(in)	::	stats	Cumulative solver and PCG convergence metrics.
real(kind=rk),	intent(in)	::	max_dT_chem	Maximum temperature increment during the chemistry advance [K].
real(kind=rk),	intent(in)	::	max_dY_chem	Maximum temperature increment during the chemistry advance [K]. Maximum species mass fraction increment during the chemistry advance.
real(kind=rk),	intent(in)	::	max_rel_drho_chem	Maximum temperature increment during the chemistry advance [K]. Maximum species mass fraction increment during the chemistry advance. Maximum relative density change during the chemistry advance.
real(kind=rk),	intent(in)	::	ke_before	Kinetic energy of the domain before the trial timestep [J].
real(kind=rk),	intent(in)	::	ke_after	Kinetic energy of the domain before the trial timestep [J]. Trial domain kinetic energy after the trial timestep [J].
real(kind=rk),	intent(in)	::	cfl_after	Kinetic energy of the domain before the trial timestep [J]. Trial domain kinetic energy after the trial timestep [J]. Peak domain CFL value computed at the trial state.
type(flow_fields_t),	intent(in)	::	fields	Struct containing current flow fields (u, p).
type(energy_fields_t),	intent(in)	::	energy	Struct containing current thermodynamic/energy fields (T, h).
type(transport_properties_t),	intent(in)	::	transport	Struct containing current transport/thermo properties (rho, mu).
logical,	intent(out)	::	reject	Output logical indicating if the trial step must be rejected.
character(len=*),	intent(out)	::	reason	Output character string describing the first tripped safety gate.

Source Code

   subroutine evaluate_step_rejection(params, stats, max_dT_chem, max_dY_chem, max_rel_drho_chem, &
                                      ke_before, ke_after, cfl_after, fields, energy, transport, &
                                      reject, reason)
      use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
      type(case_params_t), intent(in) :: params
      type(solver_stats_t), intent(in) :: stats
      real(rk), intent(in) :: max_dT_chem, max_dY_chem, max_rel_drho_chem
      real(rk), intent(in) :: ke_before, ke_after, cfl_after
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      logical, intent(out) :: reject
      character(len=*), intent(out) :: reason

      real(rk) :: ke_growth, cfl_limit

      reject = .false.
      reason = 'accepted'
      if (.not. params%enable_step_rejection) return

      ! Option B.1: NaN and Infinity Sanitation
      if (params%reject_on_nan) then
         if (any(ieee_is_nan(fields%u)) .or. any(abs(fields%u) > 1.0e8_rk)) then
            reject = .true.
            reason = 'NaN or Inf detected in velocity field u'
            return
         end if
         if (any(ieee_is_nan(fields%p)) .or. any(abs(fields%p) > 1.0e12_rk)) then
            reject = .true.
            reason = 'NaN or Inf detected in pressure field p'
            return
         end if
         if (allocated(energy%T)) then
            if (any(ieee_is_nan(energy%T)) .or. any(energy%T > 1.0e5_rk)) then
               reject = .true.
               reason = 'NaN or Inf detected in temperature field T'
               return
            end if
         end if
         if (allocated(transport%rho)) then
            if (any(ieee_is_nan(transport%rho))) then
               reject = .true.
               reason = 'NaN detected in density field rho'
               return
            end if
         end if
      end if

      ! Option B.2: Chemistry Limits
      if (params%chemistry_max_dT_per_step > zero .and. max_dT_chem > params%chemistry_max_dT_per_step) then
         reject = .true.
         reason = 'chemistry dT exceeded committed limit'
         return
      end if
      if (params%chemistry_max_dY_per_step > zero .and. max_dY_chem > params%chemistry_max_dY_per_step) then
         reject = .true.
         reason = 'chemistry dY exceeded committed limit'
         return
      end if
      if (params%chemistry_max_rel_rho_change_per_step > zero .and. &
          max_rel_drho_chem > params%chemistry_max_rel_rho_change_per_step) then
         reject = .true.
         reason = 'chemistry density change exceeded committed limit'
         return
      end if

      ! Option B.3: CFL Overshoot
      if (params%max_cfl_overshoot_factor > zero) then
         cfl_limit = params%max_cfl * params%max_cfl_overshoot_factor
         if (cfl_after > cfl_limit) then
            reject = .true.
            reason = 'post-step CFL overshoot'
            return
         end if
      end if

      ! Option B.4: Kinetic Energy Growth
      if (params%max_KE_growth_per_step > zero .and. ke_before > params%ke_reject_floor) then
         ke_growth = ke_after / ke_before
         if (ke_growth > params%max_KE_growth_per_step) then
            reject = .true.
            reason = 'kinetic energy growth exceeded limit'
            return
         end if
      end if

      ! Option B.5: PCG Pressure Iterations Ceiling
      if (params%reject_on_pressure_max_iter .and. stats%pressure_iterations >= params%pressure_max_iter) then
         reject = .true.
         reason = 'pressure solve reached max iterations'
         return
      end if

      ! Option B.6: PCG Convergence Health Audit
      if (stats%pressure_residual > params%max_pressure_residual_limit .and. &
          stats%pressure_residual_abs > params%pressure_abs_tol) then
         reject = .true.
         reason = 'pressure solver residual exceeds stability tolerance'
         return
      end if

      ! Option B.7: Density Floor/Vacuum Guard
      if (allocated(transport%rho)) then
         if (minval(transport%rho) < params%rho * params%min_density_ratio_limit) then
            reject = .true.
            reason = 'density vacuum limit violated'
            return
         end if
      end if

      ! Option B.8: Temperature Floor and Ceiling Excursions
      if (allocated(energy%T)) then
         if (minval(energy%T) < params%min_temperature_limit) then
            reject = .true.
            reason = 'temperature fell below safety floor'
            return
         end if
         if (maxval(energy%T) > params%max_temperature_limit) then
            reject = .true.
            reason = 'temperature exceeded safety ceiling'
            return
         end if
      end if
   end subroutine evaluate_step_rejection

evaluate_step_rejection Subroutine

Contents

Source Code

subroutine evaluate_step_rejection(params, stats, max_dT_chem, max_dY_chem, max_rel_drho_chem, ke_before, ke_after, cfl_after, fields, energy, transport, reject, reason)

Uses

Arguments

Source Code