mod_small_values.t

Go to the documentation of this file.

!> Module for handling problematic values in simulations, such as negative
!> pressures
module mod_small_values
 
  implicit none
  private
 
  !> How to handle small values
  character(len=20), public :: small_values_method = "error"
 
  !> Average over this many cells in each direction
  integer, public :: small_values_daverage = 1
 
  !> trace small values in the source file using traceback flag of compiler
  logical, public :: trace_small_values=.false.
 
  !> Whether to apply small value fixes to certain variables
  logical, public, allocatable :: small_values_fix_iw(:)
 
  public :: small_values_error
  public :: small_values_average
 
contains
 
  subroutine small_values_error(wprim, x, ixI^L, ixO^L, w_flag, subname)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: wprim(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    logical, intent(in)          :: w_flag(ixi^s,1:nw)
    character(len=*), intent(in) :: subname
 
    double precision :: x^d
    integer :: iw,iiw,ix^d
 
    if (.not.crash) then
      do iw=1,nw
       {do ix^db= ixo^lim^db\}
          if(w_flag(ix^d,iw)) then
            write(*,*) "Error: small value of ", trim(prim_wnames(iw)),wprim(ix^d,iw),&
                 " encountered when call ", subname
            write(*,*) "Iteration: ", it, " Time: ", global_time, "Processor: ",mype
            write(*,*) "Location: ", x(ix^d,:)
            select case (coordinate)
              case (cartesian,cartesian_stretched,cartesian_expansion)
              case (cylindrical)
                {^ifoned
                x1=x(ix^d,1)}
                {^iftwod
                select case (phi_)
                case (2)
                   x1=x(ix^d,1)*cos(x(ix^d,2))
                   x2=x(ix^d,1)*sin(x(ix^d,2))
                case default
                   x1=x(ix^d,1)
                   x2=x(ix^d,2)
                end select}
                {^ifthreed
                x1=x(ix^d,1)*cos(x(ix^d,phi_))
                x2=x(ix^d,1)*sin(x(ix^d,phi_))
                x3=x(ix^d,z_)}
                write(*,*) "Cartesian location: ", x^d
              case (spherical)
                x1=x(ix^d,1){^nooned*sin(x(ix^d,2))}{^ifthreed*cos(x(ix^d,3))}
                {^iftwod
                x2=x(ix^d,1)*cos(x(ix^d,2))}
                {^ifthreed
                x2=x(ix^d,1)*sin(x(ix^d,2))*sin(x(ix^d,3))
                x3=x(ix^d,1)*cos(x(ix^d,2))}
                write(*,*) "Cartesian location: ", x^d
              case default
                write(*,*) "No converter for coordinate=",coordinate
            end select
            write(*,*) "Cell number: ", ix^d
            do iiw=1,nw
              write(*,*) trim(prim_wnames(iiw)),": ",wprim(ix^d,iiw)
            end do
            ! use erroneous arithmetic operation to crash the run
            if(trace_small_values) write(*,*) sqrt(wprim(ix^d,iw)-bigdouble)
            write(*,*) "Saving status at the previous time step"
            crash=.true.
          end if
       {end do^D&\}
      end do
    end if
  end subroutine small_values_error
 
  subroutine small_values_average(ixI^L, ixO^L, w, x, w_flag, windex)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    logical, intent(in)             :: w_flag(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s, 1:nw)
    double precision, intent(in)    :: x(ixi^s, 1:ndim)
    integer, optional, intent(in)   :: windex
    integer                         :: iw, kxo^l, ix^d, i, nwstart, nwend
 
    if(present(windex)) then
      nwstart=windex
      nwend=windex
    else
      nwstart=1
      nwend=nw
    end if
 
    do iw = nwstart, nwend
     {do ix^db= ixo^lim^db\}
      ! point with local failure identified by w_flag
        if (w_flag(ix^d,iw)) then
          ! verify in cube with border width small_values_daverage the presence of
          ! cells where all went ok
          do i = 1, max(small_values_daverage, 1)
            {kxomin^d= max(ix^d-i, iximin^d);
            kxomax^d= min(ix^d+i, iximax^d);\}
            ! in case cells are fine within smaller cube than 
            ! the userset small_values_daverage: use that smaller cube
            if(any(w_flag(kxo^s,iw) .eqv. .false.)) exit
          end do
 
          if(any(w_flag(kxo^s,iw) .eqv. .false.)) then
            ! within surrounding cube, cells without problem were found
 
            ! faulty cells are corrected by averaging here
            ! only average those which were ok and replace faulty cells
            if(small_values_fix_iw(iw)) then
              w(ix^d, iw) = sum(w(kxo^s, iw), w_flag(kxo^s,iw) .eqv. .false.)&
                   / count(w_flag(kxo^s,iw) .eqv. .false.)
            end if
          else
            write(*,*) "no cells without error were found in cube of size", & 
                 small_values_daverage
            write(*,*) "at location:", x(ix^d, 1:ndim)
            write(*,*) "at index:", ix^d
            write(*,*) "w numer:", iw
            !write(*,*) "Saving status at the previous time step"
            !crash=.true.
            write(*,*) "replace with small values"
            if(iw==iw_e) w(ix^d, iw)=small_pressure
            if(iw==iw_rho) w(ix^d, iw)=small_density
          end if
        end if
     {enddo^d&\}
    end do
  end subroutine small_values_average
 
end module mod_small_values