precision.f90 Source File

!> @file precision.f90
!> @brief Working precision and common constants for the solver.

module precision
  use iso_fortran_env, only: real64
  implicit none
  private
  public :: wp
  public :: safe_vel

  !> Working precision kind parameter.
  !! Change to real32 or real128 to switch the entire solver's precision.
  integer, parameter :: wp = real64

contains

  !> FP-trap-safe velocity = momentum / density for output/diagnostics.
  !!
  !! Returns 0 for a vacuum/non-physical cell (rho <= 0) WITHOUT performing the
  !! division, so it cannot raise a zero/overflow/invalid SIGFPE under
  !! -ffpe-trap=invalid,zero,overflow. For valid rho > 0 the result is bit-for-bit
  !! identical to a plain mom/rho, so golden-master output is unchanged.
  !!
  !! NB: deliberately an if/else, not merge() — merge evaluates BOTH arguments
  !! and would still divide by zero in the rho <= 0 branch.
  !!
  !! @param[in] mom  Momentum component (rho*u).
  !! @param[in] rho  Density.
  !! @return Velocity, or 0 when rho <= 0.
  pure elemental function safe_vel(mom, rho) result(v)
    real(wp), intent(in) :: mom, rho
    real(wp) :: v
    v = 0.0_wp
    if (rho > 0.0_wp) v = mom / rho
  end function safe_vel

end module precision