Commits

Harald Klimach committed feb4399

Added the tem_condition_module

  • Participants
  • Parent commits 28a861e

Comments (0)

Files changed (1)

File source/tem_condition_module.f90

+! See copyright notice in the COPYRIGHT file.
+!> This module contains type definition and simple routine to load conditions
+!! Condition type has a threshold and an operator against which a quantity
+!! can be compared. This can be used for flow properties based geometry
+!! changes as well as for assessing the convergence 
+!!
+!! \author Kartik Jain
+module tem_condition_module
+  use flu_binding, only: flu_state
+  use aotus_module, only: flu_State, aot_get_val
+  use aot_table_module, only: aot_table_open, aot_table_close,                 &
+    &                         aot_table_length, aot_get_val
+  use env_module, only: rk, labelLen
+  use tem_tools_module, only: tem_write
+  use tem_aux_module, only: tem_abort
+
+  implicit none
+  private
+
+  public :: tem_condition_type, tem_convergence_type
+  public :: tem_load_convergence
+  public :: tem_convergence_evaluate
+  public :: tem_load_cond_single
+  public :: comparator
+
+
+  !> Datatype containing different conditions to be checked
+  !! Currently only threshold and operator are defined
+  type tem_condition_type
+    !> Contains the threshold value defined in lua file
+    real(kind=rk) :: threshold
+    !> Contains the operator defined in lua file
+    character(len=2) :: operation
+    !> The interval over convergence will be assessed i.e. over how 
+    !! many time steps we want to measure convergence and make sure that
+    !! we have got a straight line or a standard deviation from the 
+    !! specified threshold
+    integer :: conv_inter
+  end type tem_condition_type
+  
+  
+  !> The convergence type which contains convergence flag
+  !! and an instance of the condition type 
+  type tem_convergence_type
+    !> The convergence flag which is set to true if convergence is defined
+    logical :: active  = .false.
+    !> norm type
+    character(len=labelLen) :: norm
+    integer :: iNorm
+    !> number of defined conditions
+    integer :: nConditions
+    !> An instance of the condition type
+    type(tem_condition_type), allocatable :: cond(:)
+    !> state field holding the reference values for the 
+    real(kind=rk), allocatable :: lastState(:)
+    !> nEntries in the state array
+    integer :: nEntries_lastState
+    !> number of performed convergence checks
+    !! corresponds to the entry in the lastState array
+    integer :: nChecks
+    !> absolute Error (.true.) or relative Error( .false.)
+    logical :: absoluteError
+  end type tem_convergence_type
+
+  integer, parameter :: norm_simple = 1
+  integer, parameter :: norm_average = 2
+
+contains
+
+
+!******************************************************************************!
+  !> Load the convergence definition table 
+  !! The convergence object must be part of a tracking object, for which the 
+  !! format has been set to format = 'convergence'
+  !! In the convergence table, you then must define a norm:
+  !! - simple: just check against the state value of the last check, and reach
+  !!           convergence if below the defined threshold
+  !! - average: build the average over a defined set of last checks with nvals
+  !!           stop, if the difference to the current state value is below the given
+  !!           threshold
+  !!   - nvals: define, how many last checks should be taken into account for
+  !!            averaging procedure
+  !! The error is by default calculated to be a relative error. If an absolute
+  !! error is desired, choose absolute=true in the convergence object
+  !! 
+  !! The stopping criterion is defined as a general condition object, where the
+  !! threshold and the operator has to be given
+  !! condition = { threshold = 1.E-6, operator = '<' }
+  !! 
+  !! A sample tracking object with a convergence definition  can look as follows
+  !! \begin{verbatim}
+  !! tracking = {
+  !!   { label = 'convergence', 
+  !!     variable = {'density'}, 
+  !!     shape= { kind = 'canoND',
+  !!              object={origin ={0., 0., 0.} } },
+  !!     time = { min = 0, max = tmax, interval = tmax/1000 },
+  !!     format = 'convergence',
+  !!     convergence = {norm ='average', nvals = 100, 
+  !!          condition = {threshold = 2.0e-10, operator = '<='}
+  !!    }
+  !! \end{verbatim}
+  subroutine tem_load_convergence( me, conf, parent, verbose)
+    !--------------------------------------------------------------------------
+    type(tem_convergence_type), intent(inout) :: me
+    type(flu_state), intent(in)             :: conf
+    integer, intent(in)                     :: parent
+    integer, intent(in)                     :: verbose
+    !--------------------------------------------------------------------------
+    integer :: iError
+    integer :: thandle
+    !--------------------------------------------------------------------------
+    ! open the convergence table
+    call aot_table_open(L=conf, parent=parent, thandle=thandle, &
+                        key='convergence')
+    me%active = .true.
+    ! get the kind of the norm
+    call aot_get_val( L=conf, thandle=thandle, val=me%norm,                 &
+      &               ErrCode=iError, key='norm', default = 'simple' )
+    call aot_get_val( L=conf, thandle=thandle, val=me%nEntries_lastState, &
+      &               ErrCode=iError, key='nvals', default = 1 )
+    call aot_get_val( L=conf, thandle=thandle, val=me%absoluteError, &
+      &               ErrCode=iError, key='absolute', default = .false. )
+    call tem_load_condition( me = me%cond,                              &
+      &                      conf = conf, parent = thandle,             &
+      &                      verbose = verbose)
+    me%nConditions = size( me%cond )
+    me%nChecks     = 0
+    call tem_write(' loaded convergence with nConditions', me%nConditions ) 
+    call tem_write('    Norm : '//trim(me%norm )) 
+    if( me%absoluteError ) then
+      call tem_write('    absolute error')
+    else
+      call tem_write('    relative error')
+    end if
+    call tem_write('    nVal : ', me%nEntries_lastState )
+    call aot_table_close(L=conf, thandle=thandle )
+
+    ! Allocate some arrays
+    select case( trim( me%norm ))
+    case( 'simple' )
+      me%iNorm = norm_simple
+    case( 'average')
+      me%iNorm = norm_average
+    case default
+      call tem_write( 'Error: norm name is not correct' )
+      call tem_abort
+    end select
+    allocate( me%lastState( me%nEntries_lastState ))
+    me%lastState(:) = huge( me%lastState(1) )/real( me%nEntries_lastState, kind=rk)
+ 
+  end subroutine tem_load_convergence
+!******************************************************************************!
+
+!******************************************************************************!
+  !> Load the condition table in case of convergence
+  !! \Author: Kartik Jain
+  subroutine tem_load_condition(me, conf, parent, verbose)
+    !--------------------------------------------------------------------------
+    type(tem_condition_type), allocatable, intent(inout) :: me(:)
+    type(flu_state), intent(in)             :: conf
+    integer, intent(in)                     :: parent
+    integer, intent(in), optional           :: verbose
+    !--------------------------------------------------------------------------
+    integer :: iError
+    integer :: cond_handle          !< handle for the condition table
+    integer :: sub_cond_handle      !< handle for subtables inside condition
+    integer :: nCond                !< number of conditions
+    integer :: iCond                !< index for condition loop
+    integer :: verbose_local        !< index for condition loop
+    !--------------------------------------------------------------------------
+    if( present( verbose )) then
+      verbose_local = verbose 
+    else
+      verbose_local = 1
+    end if
+    !! Open the condition table
+    call aot_table_open(L=conf, parent=parent, thandle=cond_handle, &
+                        key='condition')
+    !! The tables inside condition table should be equal to the nVars
+    !! If thats not the case we return an error message
+    nCond = aot_table_length( L=conf, thandle=cond_handle)
+      !! check single or multiple table
+      call aot_table_open( L=conf, parent=cond_handle, thandle=sub_cond_handle, &
+        &                  pos=1 )
+      
+      if (sub_cond_handle .eq. 0) then
+        ! just one table within
+        allocate ( me(1) )
+        call tem_load_cond_single( me(1), conf, cond_handle, verbose)
+      else 
+        ! IF there are more tables within condition
+        allocate ( me(nCond) )
+        do iCond = 1, nCond
+          call aot_table_open( L=conf, parent=cond_handle, thandle=sub_cond_handle, &
+            &                  pos=iCond)
+          call tem_load_cond_single(me(iCond), conf, sub_cond_handle, verbose)
+          call aot_table_close(L=conf, thandle=sub_cond_handle )
+        end do
+      end if ! sub condition check
+      call aot_table_close(L=conf, thandle=cond_handle )
+
+  end subroutine tem_load_condition
+!******************************************************************************!
+
+
+!******************************************************************************!
+  !> Load the conditions for geomIncr and convergence check within tracking 
+  !! conditions mean the operator and threshold against which the macroscopic
+  !! variable has to be compared 
+  subroutine tem_load_cond_single(cond, conf, thandle, verbose)
+    !--------------------------------------------------------------------------
+    type(tem_condition_type), intent(inout) :: cond
+    type(flu_state), intent(in)             :: conf
+    integer, intent(in)                     :: thandle
+    integer, intent(in), optional           :: verbose
+    !--------------------------------------------------------------------------
+    integer :: iError
+    integer :: verbose_local
+    !--------------------------------------------------------------------------
+    if( present( verbose )) then
+      verbose_local = verbose 
+    else
+      verbose_local = 1
+    end if
+
+    call aot_get_val( L=conf, thandle=thandle, val=cond%threshold,     &
+      &               ErrCode=iError, key='threshold' )
+
+    call aot_get_val( L=conf, thandle=thandle, val=cond%operation,          &
+      &               ErrCode=iError, key='operator' )
+
+  end subroutine tem_load_cond_single
+!******************************************************************************!
+
+
+!******************************************************************************!
+  !> return a logical if the input relation holds where the relation is:
+  !! val <operation> threshold
+  !! \Author: Kartik Jain
+  function comparator( val, operation, threshold) result(comp)
+    !-------------------------------------------------------------------------- 
+    real(kind=rk), intent(in)      :: val
+    character(len=2), intent(in)   :: operation
+    real(kind=rk), intent(in)      :: threshold
+    logical                        :: comp !< return value
+    !------------------------------------------------------------------------- 
+    comp = .false.
+
+    select case( trim(operation))
+      case('<') 
+        if (val .lt. threshold) then
+          comp = .true.
+        else
+          comp = .false.
+        end if
+      case('<=')
+        if (val .le. threshold) then
+          comp = .true.
+        else
+          comp = .false.
+        end if
+      case('>')
+        if (val .gt. threshold) then
+          comp = .true.
+        else
+          comp = .false.
+        end if
+      case('>=')
+        if (val .ge. threshold) then
+          comp = .true.
+        else
+          comp = .false.
+        end if
+      case('=')
+        if (val .eq. threshold) then
+          comp = .true.
+        else
+          comp = .false.
+        end if        
+    end select
+
+  end function comparator
+!******************************************************************************!
+
+
+!******************************************************************************!
+  !> Evaluate if the convergence was achieved
+  subroutine tem_convergence_evaluate( me, state, achieved )
+    !--------------------------------------------------------------------------
+    type(tem_convergence_type), intent(inout) :: me
+    real(kind=rk), intent(in)             :: state(:)
+    logical, intent(out)                  :: achieved
+    !--------------------------------------------------------------------------
+    integer :: iCondition
+    logical :: compare
+    real(kind=rk) :: residual
+    !--------------------------------------------------------------------------
+    achieved = .false.
+    do iCondition = 1, me%nConditions
+      ! Compare the results against threshold
+      !norm = abs(state(1) - me%result_prev)
+      residual = evaluate_residual( me = me, state = state(:) )
+      compare = comparator(                                                    &
+        &           val = residual,                                            &
+        &           operation = me%cond( iCondition )%operation,               &
+        &           threshold = me%cond( iCondition )%threshold )
+   
+      ! update the convergence achieved if the last compare was succesful
+      achieved = ( compare .or. achieved )
+    end do
+
+  end subroutine tem_convergence_evaluate
+!******************************************************************************!
+
+
+!******************************************************************************!
+  !> evaluate the residual
+  !! For relative errors (defined in convergence%absoluteError ), the result is
+  !! divided by the current status value
+  function evaluate_residual( me, state ) result( res )
+    !--------------------------------------------------------------------------
+    type(tem_convergence_type), intent(inout) :: me
+    real(kind=rk), intent(in)              :: state(:)
+    real(kind=rk) :: res
+    !--------------------------------------------------------------------------
+    integer :: pos_state, pos_lastState
+    real(kind=rk) :: average, factor
+    !--------------------------------------------------------------------------
+    ! \todo if more entries than the first should be tracked, change this
+    ! pos_state
+    pos_state = 1
+    ! Increase the counter for the checks
+    me%nChecks = me%nChecks + 1
+    ! Reset the result
+    res = huge( res )
+    if( me%absoluteError ) then
+      ! For absolute error, just take the difference between value and reference
+      ! error = ( val - ref )
+      factor = 1._rk
+    else
+      ! For relative errors, divide by the value 
+      ! error = ( val - ref ) / val
+      factor = 1._rk/state(pos_state)
+    end if
+
+    select case (me%iNorm )
+    case( norm_simple )
+      res = abs((state( pos_state ) - me%lastState(1))*factor )
+      ! update the result at t-1 to t as when we arrive at t+1, it will
+      ! be required
+      me%lastState(1) = state( pos_state )
+    case( norm_average )
+      average = sum( me%lastState)/real( me%nEntries_lastState, kind=rk)
+      res = abs( (state( pos_state ) - average )*factor )
+      pos_lastState = mod( me%nChecks -1, me%nEntries_lastState ) + 1
+      me%lastState( pos_lastState ) = state( pos_state )
+      !write(*,*) 'lastState', me%lastState
+      !write(*,*) 'state:', state(pos_state),'average', average, 'res', res, &
+      !  &        'pos last', pos_lastState
+    end select
+
+  end function evaluate_residual
+!******************************************************************************!
+
+end module tem_condition_module