! See copyright notice in the COPYRIGHT file.
!> author: Kartik Jain, Kannan Masilamani
!! This module contains type definition and routines for convergence to set
!! steady state status
!!
!! Convergence is loaded when steady_state is set to true in abort_criteria.
!! Convergence is achieved if all variables defined per convergence table
!! are converged. To check for convergence reduction and condition
!! must be defined for each variable.
!!
module tem_convergence_module

  ! include treelm modules
  use env_module,              only: rk, labelLen, io_buffer_size
  use tem_aux_module,          only: tem_abort
  use tem_logging_module,      only: logUnit
  use tem_tools_module,        only: tem_horizontalSpacer, upper_to_lower
  use tem_comm_env_module,     only: tem_comm_env_type
  use tem_status_module,       only: tem_stat_steady_state, tem_status_type
  use tem_time_module,         only: tem_time_type
  use tem_timeControl_module,  only: tem_timeControl_type,                     &
    &                                tem_timeControl_load,                     &
    &                                tem_timeControl_dump,                     &
    &                                tem_timeControl_out,                      &
    &                                tem_timeControl_check,                    &
    &                                tem_timeControl_reset_trigger
  use tem_varSys_module,       only: tem_varSys_type, tem_get_element_chunk, &
    &                                tem_get_point_chunk
  use tem_varMap_module,       only: tem_varMap_type, &
    &                                tem_create_varMap
  use treelmesh_module,        only: treelmesh_type
  use tem_bc_prop_module,      only: tem_bc_prop_type
  use tem_reduction_spatial_module, only: tem_reduction_spatial_type,        &
    &                                     tem_reduction_spatial_close,       &
    &                                     tem_reduction_spatial_config_type, &
    &                                     tem_load_reduction_spatial,        &
    &                                     tem_reduction_spatial_init,        &
    &                                     tem_reduction_spatial_open,        &
    &                                     tem_reduction_spatial_append,      &
    &                                     tem_reduction_spatial_out
  use tem_shape_module,        only: tem_shape_type, tem_load_shape,           &
    &                                tem_shape_out
  use tem_subTree_module,      only: tem_create_subTree_of
  use tem_subTree_type_module, only: tem_subTree_type
  use tem_solveHead_module,    only: tem_solveHead_type
  use tem_stencil_module,      only: tem_stencilHeader_type
  use tem_condition_module,    only: tem_condition_type, tem_load_condition,   &
    &                                tem_comparator, tem_condition_out
  use tem_geometry_module,     only: tem_BaryOfId

  ! include aotus modules
  use aotus_module,     only: flu_State, aot_get_val, aoterr_Fatal
  use aot_table_module, only: aot_table_open, aot_table_close,                 &
    &                         aot_table_length, aot_get_val
  use aot_out_module,   only: aot_out_type, aot_out_open,  aot_out_close,      &
    &                         aot_out_val, aot_out_open_table,                 &
    &                         aot_out_close_table

  implicit none
  private

  public :: tem_condition_type, tem_convergence_type
  public :: tem_convergence_load
  public :: tem_convergence_dump
  public :: tem_convergence_out
  public :: tem_convergence_check
  public :: tem_init_convergence
  public :: tem_convergence_reset

  !> Convergence description loaded from config file
  type tem_convergenceHeader_type
    !> convergence kind
    character(len=labelLen) :: norm
    !> number of defined conditions
    integer :: nConditions
    !> An instance of the condition type for each variable
    type(tem_condition_type), allocatable :: cond(:)
    !> Number of last values to check for convergence
    integer :: nLastVals
    !> absolute Error (.true.) or relative Error( .false.)
    logical :: absoluteError
    !> array of variable labels to check for convergence
    character(len=labelLen), allocatable :: varName(:)
    !> number of variables to check for convergence
    !! i.e size(varName)
    integer :: nRequestedVars
    !> stores time control parameters
    type(tem_timeControl_type) :: timeControl
    !> convergence shapes
    type(tem_shape_type), allocatable  :: geometry(:)

    !> reduction config
    type(tem_reduction_spatial_config_type) :: redSpatial_config

    !> Logic to decide to use get_point or get_element to dump data
    logical :: useGetPoint

    !> Number of dofs to check for convergence if useGetPoint = .false.
    integer :: nDofs

  end type tem_convergenceHeader_type

  !> The convergence type which contains convergence flag
  !! and an instance of the condition type
  type tem_convergence_type
    !> Convergence header info
    type(tem_convergenceHeader_type) :: header
    !> norm kind
    integer :: norm_kind
    !> state field holding the reference values for the nScalars
    !! size: nLastVals, nScalars
    real(kind=rk), allocatable :: lastState(:,:)
    !> number of performed convergence checks
    !! corresponds to the entry in the lastState array
    integer :: nChecks
    !> Process description to use for the output.
    !! Might be only a subset of the global communicator
    type(tem_comm_env_type)  :: proc
    !> Contains name and position of variables in global varSys
    type(tem_varMap_type) :: varMap
    !> sub-tree resulting from the elements within the convergence shape
    !! The sub-tree also holds the sub-communicator
    type(tem_subTree_type) :: subTree
    !> number of elements that fit in the buffer
    integer :: chunkSize
    !> number of chunks per output
    integer :: nChunks

    !> The number of dofs for each scalar variable of the equation system
    integer :: nDofs

    !> spatial reduction for each variable
    type(tem_reduction_spatial_type), allocatable :: redSpatial(:)

  end type tem_convergence_type

  !> Convergence norm kinds
  integer, parameter :: norm_simple = 1
  integer, parameter :: norm_average = 2

  interface assignment(=)
    module procedure Copy_convergence
  end interface

  interface tem_convergence_dump
    module procedure tem_convergence_dump_vector
    module procedure tem_convergence_dump_single
  end interface tem_convergence_dump

  interface tem_convergence_out
    module procedure tem_convergence_out_vector
    module procedure tem_convergence_out_single
  end interface tem_convergence_out

contains


! **************************************************************************** !
  !> Load the convergence definition table
  !! The convergence object must be part of a convergence 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
  !!~~~~~~~~~~~~~~~~~~~~~{.lua}
  !! condition = { threshold = 1.E-6, operator = '<' }
  !!~~~~~~~~~~~~~~~~~~~~~
  !!
  !! A sample convergence object with a convergence definition  can look as follows
  !! (within time_control table)
  !!~~~~~~~~~~~~~~~~~~~~~{.lua}
  !! abort_criteria = {
  !!  stop_file = 'stop',
  !!  steady_state = false,
  !!   convergence = {
  !!     variable = {'pressure','velocity'},
  !!     shape = {kind = 'all'},
  !!     time_control = {min = {iter=0}, max = {iter=tmax}, interval = {iter=10*dt}},--},
  !!     reduction = {'average','average'},
  !!     norm='average', nvals = 100, absolute = true,
  !!     condition = {
  !!        { threshold = 1.e-15, operator = '<=' },
  !!        { threshold = 1.e-12, operator = '<=' }
  !!     }
  !!   }
  !! }
  !!~~~~~~~~~~~~~~~~~~~~~
  !!
  !! Or another sample:
  !!~~~~~~~~~~~~~~~~~~~~~{.lua}
  !! abort_criteria = {
  !!   stop_file     = 'stop',
  !!   steady_state  = true,
  !!   convergence   = {
  !!     variable = {'pressure_phy'},
  !!     shape = {
  !!       {kind = 'canoND', object = {origin ={0.15-dx,0.2,zpos} }},
  !!       {kind = 'canoND', object = {origin ={0.25+dx,0.2,zpos} }}
  !!     },
  !!     time_control = {min = 0, max = tmax, interval = 10*dt},
  !!     reduction = {'average'},
  !!     norm      = 'average',
  !!     nvals     = 50,
  !!     absolute  = true,
  !!     condition = { threshold = 1.e-10, operator = '<=' }
  !!   }
  !! },
  !!~~~~~~~~~~~~~~~~~~~~~
  !!
  subroutine tem_convergence_load(me, conf, parent, steady_state)
    ! ---------------------------------------------------------------------------
    !> list of the convergence entities to create
    type( tem_convergence_type ), allocatable, intent(out) :: me(:)
    !> general control parameters
    !> handle of the lua config file
    type( flu_state ) :: conf
    !> if the convergence table is a child-table of some other table,
    !! use the parent as a reference
    integer, optional :: parent
    !> Steady flag in abort_criteria to check for convergence
    logical, intent(inout) :: steady_state
    ! ---------------------------------------------------------------------------
    integer :: conv_handle, sub_handle
    integer :: iConv, nConv
    ! ---------------------------------------------------------------------------

    ! Read the number of convergences in the lua file
    call aot_table_open( L       = conf,          &
      &                  thandle = conv_handle,   &
      &                  key     = 'convergence', &
      &                  parent  = parent         )

    if (conv_handle == 0) then
      write(logUnit(1),*) 'WARNING: Abort criteria, steady state is true but'
      write(logUnit(1),*) '         No Convergence table is defined with '
      write(logUnit(1),*) '         conditions to check for steady state'
      write(logUnit(1),*) 'NOTE: Steady state is deactivated'
      steady_state = .false.
      call aot_table_close(L=conf, thandle=conv_handle)
      call tem_horizontalSpacer(fUnit=logUnit(1))
      return
    end if

    write(logUnit(1),*) 'Loading convergence for steady state...'
    ! Check whether convergence had a subtable
    ! If no, then it is a single table, load single convergence entry
    ! else load multiple tables, open convergence subtable
    call aot_table_open( L       = conf,        &
      &                  parent  = conv_handle, &
      &                  thandle = sub_handle,  &
      &                  pos     = 1            )

    ! Only single table
    if (sub_handle == 0) then
      nConv = 1
      write(logUnit(1),*) 'Convergence is a single table'
      allocate( me( nConv ) )
      call tem_load_convergenceHeader( conf       = conf,        &
        &                              sub_handle = conv_handle, &
        &                              me         = me(1)        )
      call aot_table_close(L=conf, thandle=sub_handle)
    else ! Multiple table
      call aot_table_close(L=conf, thandle=sub_handle)
      nConv = aot_table_length(L=conf, thandle=conv_handle)
      ! Allocate the defined number of convergence entities
      allocate( me( nConv ))
      write(logUnit(1),*) 'Number of Convergence entities: ', nConv

      ! Loop over all the definitions and assign the variables from the lua
      ! file on the tem_convergence_type.
      ! Inside this routine it will open convergence subtable. Each subtable
      ! contains one or more convergence variables the stuff is done in the
      ! routine tem_load_convergenceHeader
      do iConv = 1, nConv
        write(logUnit(3),*) 'Loading convergence ', iConv
        call aot_table_open( L       = conf,        &
          &                  parent  = conv_handle, &
          &                  thandle = sub_handle,  &
          &                  pos     = iConv      )
        call tem_load_convergenceHeader( conf           = conf,       &
          &                              sub_handle     = sub_handle, &
          &                              me             = me(iConv)   )
        call aot_table_close(L=conf, thandle=sub_handle)
        write(logUnit(3),*) 'Done'
      end do
    end if ! sub_handle

    call aot_table_close(L=conf, thandle=conv_handle) ! close convergence table
    call tem_horizontalSpacer(fUnit=logUnit(1))

  end subroutine tem_convergence_load
! **************************************************************************** !


! **************************************************************************** !
  !> Read the convergence variables from convergence subtables defined in
  !! configuration from the main lua file
  !!
  !! If convergence is just a single table with single convergence entry
  !! then load only one convergence log exists with
  !! one or more variables using tem_load_convergenceHeader_single.
  !! Else if convergence is table of many log then allocate log and load each
  !! log type using tem_load_convergenceHeader_single
  !! Setup the values for the convergence entities
  !!
  subroutine tem_load_convergenceHeader(me, conf, sub_handle)
    ! --------------------------------------------------------------------------
    !> list of the convergence entities to create
    type( tem_convergence_type ), intent(out) :: me
    !> handle of the lua config file
    type( flu_state ) :: conf
    !> table sub-handle for the convergence table
    integer, intent(in) :: sub_handle
    ! --------------------------------------------------------------------------
    integer :: iError            ! error flag handle
    integer, allocatable :: vError(:)
    character(len=labelLen) :: norm
    ! --------------------------------------------------------------------------

    call aot_get_val( val       = me%header%varname, &
      &               ErrCode   = vError,            &
      &               maxLength = 100,               &
      &               L         = conf,              &
      &               thandle   = sub_handle,        &
      &               key       = "variable"         )

    if ( any(btest(vError, aoterr_Fatal)) ) then
      write(logUnit(1),*) 'FATAL Error occured, while retrieving'
      write(logUnit(1),*) 'list of variables to use in convergence'
      call tem_abort()
    end if

    me%header%nRequestedVars = size(me%header%varName)

    ! load time control to output convergence
    call tem_timeControl_load( conf           = conf,                          &
      &                        parent         = sub_handle,                    &
      &                        me             = me%header%timeControl )
    call tem_timeControl_dump(me%header%timeControl, logUnit(3))

    ! load convergence object shapes like point, line, plane
    call tem_load_shape( conf = conf, parent = sub_handle,                     &
                         me = me%header%geometry )

    if( size( me%header%geometry) < 1) then
      write(logUnit(1),*)'The geometrical objects for the convergence are'//  &
        &            ' not defined correctly.'
      call tem_abort()
    end if

    ! load reductions
    call tem_load_reduction_spatial( conf   = conf,                       &
      &                              parent = sub_handle,                 &
      &                   redSpatial_config = me%header%redSpatial_config )

    if( me%header%redSpatial_config%active ) then
      ! Check if the number of reductions correspond to the number of variables
      ! in the system
      if( size( me%header%redSpatial_config%reduceType ) &
        & /= me%header%nRequestedVars ) then
        write(logUnit(1),*) 'Error: In convergence.'
        write(logUnit(1),*) 'The number of defined reductions does not '// &
          &            'correspond to the '
        write(logUnit(1),*)'number of variables in the system. '
        call tem_abort()
      end if
    else
      write(logUnit(1),*) 'Error: No spatial reduction defined.'
      write(logUnit(1),*) 'NOTE: Convergence requires spatial reduction '//&
        &                 'for each variable'
      call tem_abort()
    end if

    ! get the kind of the convergence norm
    call aot_get_val( L       = conf,                                          &
      &               thandle = sub_handle,                                    &
      &               val     = norm,                                          &
      &               ErrCode = iError,                                        &
      &               key     = 'norm',                                        &
      &               default = 'simple' )
    norm = adjustl(norm)
    norm = upper_to_lower(norm)
    me%header%norm = trim(norm)
    select case( trim( norm ))
    case( 'simple' )
      me%norm_kind = norm_simple
    case( 'average')
      me%norm_kind = norm_average
    case default
      write(logUnit(1),*) 'Error: Unknown convergence norm'
      write(logUnit(1),*) 'Solution: Supported norms '
      write(logUnit(1),*) '        * simple'
      write(logUnit(1),*) '        * average'
      call tem_abort
    end select

    ! Get number of last values to check for convergence
    call aot_get_val( L       = conf,                                          &
      &               thandle = sub_handle,                                    &
      &               val     = me%header%nLastVals,                           &
      &               ErrCode = iError,                                        &
      &               key     = 'nvals',                                       &
      &               default = 1 )

    ! type of convergence error: relative or absolute
    call aot_get_val( L       = conf,                                          &
      &               thandle = sub_handle,                                    &
      &               val     = me%header%absoluteError,                       &
      &               ErrCode = iError,                                        &
      &               key     = 'absolute',                                    &
      &               default = .false. )

    ! To decide whether to use get_point or get_element
    call aot_get_val( L       = conf,                  &
      &               thandle = sub_handle,            &
      &               val     = me%header%useGetPoint, &
      &               ErrCode = iError,                &
      &               default = .false.,               &
      &               key     = 'use_get_point'        )

    ! Get the number of Dofs to be written in the output
    ! The default is set to -1. If the dofs are not specified,
    ! all the dofs should be dumped
    call aot_get_val( L       = conf,            &
      &               thandle = sub_handle,      &
      &               val     = me%header%nDofs, &
      &               ErrCode = iError,          &
      &               default = -1,              &
      &               key     = 'ndofs'          )

    ! load condition for convergence
    call tem_load_condition( me      = me%header%cond,                         &
      &                      conf    = conf,                                   &
      &                      parent  = sub_handle )

    me%header%nConditions = size(me%header%cond)
    ! check if there is condition for each variable
    if (me%header%nConditions /= me%header%nRequestedVars) then
      write(logUnit(1),*) 'Error: Nr. of conditions \= Nr. of variables '
      write(logUnit(1),"(2(A,I0))") 'nCond: ', me%header%nConditions, &
        &                           'nVars: ', me%header%nRequestedVars
      call tem_abort()
    end if

    write(logUnit(1),"(A,I0)") ' loaded convergence with nConditions=', &
      &                 me%header%nConditions
    write(logUnit(1),*) '    Norm : '//trim(norm)
    if( me%header%absoluteError ) then
      write(logUnit(1),*)'    absolute error'
    else
      write(logUnit(1),*)'    relative error'
    end if
    write(logUnit(1),"(A,I0)")'    nVal : ', me%header%nLastVals
    write(logUnit(7),*) '  Use get_point: ', me%header%useGetPoint
    call aot_table_close(L=conf, thandle=sub_handle )

  end subroutine tem_load_convergenceHeader
! **************************************************************************** !


! **************************************************************************** !
  !> Initialize the convergence subtreee
  !!
  !! Identify, how many and which elements exist on my local process and are
  !! requested from the convergences
  !! Empty convergence  are removed, so the convergence(:) might be re-allocated
  !! ~~~~~~~~~~~~~~~~~~~~~{.lua}
  !!
  subroutine tem_init_convergence( me, tree, varSys, bc_prop, globProc,&
    &                              stencil, nDofs )
    ! ---------------------------------------------------------------------------
    !> convergence descriptions
    type(tem_convergence_type),intent(inout), allocatable :: me(:)
    !> Global mesh from which the elements are identified and then stored to
    !! sub-meshes inside the convergences
    type(treelmesh_type), intent(in) :: tree
    !> bc property that used to identify elements of certain BCs
    type( tem_bc_prop_type ), intent(in) :: bc_prop
    !> solver-provided variable systems
    type(tem_varSys_type), intent(in) :: varSys
    !> Process description to use.
    type(tem_comm_env_type), intent(in) :: globProc
    !> stencil used to create subTree of boundary type
    type(tem_stencilHeader_type), optional, intent(in) :: stencil
    !> The number of dofs for each scalar variable of the equation system
    integer, intent(in), optional :: nDofs
    ! ---------------------------------------------------------------------------
    integer :: iConv, nConv, nVars
    integer :: nChunks, chunkSize, nElems, maxComponents, nPoints
    ! temporary convergence array
    type( tem_convergence_type ),allocatable :: tempConv(:)
    ! ---------------------------------------------------------------------------

    call tem_horizontalSpacer(fUnit=logUnit(1))
    write(logUnit(1),*) 'Setting up the convergence infrastructure'
    nConv = 0


    ! Allocate the temporary convergence
    allocate(tempConv(size(me)))

    do iConv = 1, size(me)
      write(logUnit(10),"(A,I0)") 'Initializing convergence object ', iConv

      ! map variables
      ! create convergence variable position in the global varSys
      call tem_create_varMap( varname = me(iConv)%header%varname, &
        &                     varSys  = varSys,                   &
        &                     varMap  = me(iConv)%varMap          )

      ! @todo KM: If variable not found in varSys then skip that variable
      ! reduction and condition info while copying convergence
      !
      ! Terminate if number of variables to check for convergence does not
      ! match with variables found in varMap
      if (me(iConv)%varMap%varPos%nVals /= me(iConv)%header%nRequestedVars) then
        write(logUnit(1),*) 'Error: Mapping Convergences variables'
        write(logUnit(1),*) 'Variables defined in convergence '// &
          &                 'table ', iConv, ' are not found in varSys'
        call tem_abort()
      end if

      ! -----------------------------------------------------------------------
      ! identify convergence elements
      ! -----------------------------------------------------------------------
      call tem_create_subTree_of( inTree    = tree,                           &
        &                         bc_prop   = bc_prop,                        &
        &                         stencil   = stencil,                        &
        &                         subTree   = me( iConv )%subTree,            &
        &                         storePnts = me (iConv )%header%useGetPoint, &
        &                         inShape   = me( iConv )%header%geometry     )

      ! get rid of the empty convergence in order
      ! to avoid empty writes to disk
      if ( me(iConv)%subTree%useGlobalMesh ) then

        ! set convergence communicator
        me(iConv)%proc = globProc

        nConv = nConv + 1
        tempConv( nConv ) = me( iConv )

      else if ( me(iConv)%subTree%nElems > 0 ) then

        nConv = nConv + 1 ! Increase number of log

        ! set convergence communicator
        me(iConv)%proc%comm      = me(iConv)%subTree%global%comm
        me(iConv)%proc%rank      = me(iConv)%subTree%global%myPart
        me(iConv)%proc%comm_size = me(iConv)%subTree%global%nParts
        me(iConv)%proc%root      = 0

        ! Copy all entries from the log derived type to the temporary one.
        ! this might not work on all compilers!!
        ! This assignment is realized by operator overloader Copy_convergence
        tempConv( nConv ) = me( iConv )
      end if ! useGlobalMesh
    end do  ! nConv

    deallocate(me)
    allocate( me(nConv) )

    do iConv = 1, nConv
      ! Copy the stuff from the temporary track
      me(iConv) = tempConv(iConv)

      ! number of variables in varMap
      nVars = me(iConv)%varMap%varPos%nVals

      ! nDofs is valid only for get_element
      if (me(iConv)%header%useGetPoint) then
        me(iConv)%nDofs = 1
      else
        if (present(nDofs)) then
          ! out_config%nDofs is set to -1 if unspecied
          ! in the config file. In this case all the dof's
          ! should be dumped
          if (me(iConv)%header%nDofs < 0) then
            me(iConv)%nDofs = nDofs
          else
            ! Otherwise the number of dofs dumped should
            ! be what's specified in the config
            me(iConv)%nDofs = me(iConv)%header%nDofs
          end if
        else
          me(iConv)%nDofs = 1
        end if
      end if

      if (me(iConv)%subTree%useGlobalMesh) then
        nElems = tree%nElems
        nPoints = tree%nElems
      else
        nElems = me(iConv)%subTree%nElems
        nPoints = me(iConv)%subTree%nPoints
      end if

      ! max nComponent in current convergence variables
      maxComponents = maxval(varSys%method%val(me(iConv)%varMap               &
        &                                     %varPos%val(:nVars))%nComponents)

      if (me(iConv)%header%useGetPoint) then
        chunkSize = min(io_buffer_size/maxComponents, nPoints)
      else
        chunkSize = min(io_buffer_size/(maxComponents*me(iConv)%nDofs), nElems)
      end if
      if ( (nElems > 0) .and. (chunkSize == 0) ) then
        write(logUnit(0),*)'Error in init_convergence: '
        write(logUnit(0),*) 'The chosen io_buffer_size of ', io_buffer_size
        write(logUnit(0),*) 'is too small for evaluating ', maxComponents
        write(logUnit(0),*) 'scalar values'
        write(logUnit(0),*) 'Please increase the io_buffer_size to at &
          & least ', real(maxComponents*me(iConv)%nDofs) / real(131072), ' MB!'
        call tem_abort()
      end if

      nChunks = 0
      if (chunkSize>0) then
        if (me(iConv)%header%useGetPoint) then
          nChunks = ceiling(real(nPoints, kind=rk)/real(chunkSize, kind=rk))
        else
          nChunks = ceiling(real(nElems, kind=rk)/real(chunkSize, kind=rk))
        end if
      else
        nChunks = 0
      end if

      me(iConv)%nChunks = nChunks
      me(iConv)%chunkSize = chunkSize
      me(iConv)%nChecks = 0

      ! Initialize reduction
      call tem_reduction_spatial_init(                                 &
        &      me                = me(iConv)%redSpatial,               &
        &      redSpatial_config = me(iConv)%header%redSpatial_config, &
        &      varSys            = varSys,                             &
        &      varPos            = me(iConv)%varMap%varPos%val(:nVars) )

      ! Allocate some arrays
      allocate( me(iConv)%lastState( me(iConv)%header%nLastVals, &
        &                           me(iConv)%varMap%nScalars ) )
      me(iConv)%lastState = huge( me(iConv)%lastState(1,1) )          &
        &                 / real( me(iConv)%header%nLastVals, kind=rk )

    end do

    deallocate(tempConv)

    call tem_horizontalSpacer(fUnit=logUnit(1))

  end subroutine tem_init_convergence
! **************************************************************************** !


! **************************************************************************** !
  !> This routine resets convergence lastState and nChecks
  subroutine tem_convergence_reset( me )
    ! ---------------------------------------------------------------------------
    !> convergence descriptions
    type(tem_convergence_type),intent(inout) :: me(:)
    ! ---------------------------------------------------------------------------
    integer :: iConv
    ! ---------------------------------------------------------------------------
    do iConv = 1, size(me)
      me(iConv)%nChecks = 0
      me(iConv)%lastState = huge( me(iConv)%lastState(1,1) )          &
        &                 / real( me(iConv)%header%nLastVals, kind=rk )
      call tem_timeControl_reset_trigger(me(iConv)%header%timeControl)
    end do

  end subroutine tem_convergence_reset
! **************************************************************************** !

! **************************************************************************** !
  !> This routine runs over each convergence object and check for convergence
  !! of each requested quantities timeControl is active on current time
  subroutine tem_convergence_check( me, time, status, varSys, tree )
    ! ---------------------------------------------------------------------------
    !> convergence descriptions
    type(tem_convergence_type), target, intent(inout)   :: me(:)
    !> current simulation time
    type(tem_time_type), intent(in)     :: time
    !> Status bits
    type(tem_status_type), intent(inout) :: status
    !> global variable system
    type(tem_varSys_type), intent(in)         :: varSys
    !> global tree
    type(treelmesh_type ), intent(in)         :: tree
    ! ---------------------------------------------------------------------------
    integer :: iConv
    logical :: triggered
    real(kind=rk), allocatable :: res(:)
    ! ---------------------------------------------------------------------------
    allocate(res(io_buffer_size))

    ! Run over all convergence objects
    do iConv = 1,size(me)

      ! Skip this convergence object, if there are no entries in the
      ! variable system
      if (me( iConv )%varMap%nScalars < 1) cycle

      ! Determine, if it is now the time to perform the convergence for each
      ! convergence object
      ! check time control and update it if triggered
      call tem_timeControl_check( me        = me( iConv )%header%timeControl, &
        &                         now       = time,                           &
        &                         comm      = me( iConv )%proc%comm,          &
        &                         triggered = triggered                       )

      ! check convergence if timeControl is triggered
      if( triggered )then
        if (me(iConv)%header%useGetPoint) then
          call tem_convergence_check_point( me     = me(iConv), &
            &                               time   = time,      &
            &                               status = status,    &
            &                               varSys = varSys,    &
            &                               tree   = tree,      &
            &                               res    = res        )
        else
          call tem_convergence_check_element( me     = me(iConv), &
            &                                 time   = time,      &
            &                                 status = status,    &
            &                                 varSys = varSys,    &
            &                                 tree   = tree,      &
            &                                 res    = res        )
        end if
      end if  ! do convergence? interval, tmin, tmax check
    end do ! iConv

    deallocate(res)
  end subroutine tem_convergence_check
! **************************************************************************** !


! **************************************************************************** !
  !> This routine runs over convergence check using get_element interface
  subroutine tem_convergence_check_element( me, time, status, varSys, tree, &
    &                                       res )
    ! ---------------------------------------------------------------------------
    !> convergence descriptions
    type(tem_convergence_type), target, intent(inout)   :: me
    !> current simulation time
    type(tem_time_type), intent(in)     :: time
    !> Status bits
    type(tem_status_type), intent(inout) :: status
    !> global variable system
    type(tem_varSys_type), intent(in)         :: varSys
    !> global tree
    type(treelmesh_type ), intent(in)         :: tree
    !> Output data
    !! size: io_buffer_size
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------------
    integer :: nVars, nElems, nScalars, elemOff, nChunkElems
    integer :: iElem, iChunk
    integer :: buf_start, buf_end
    integer, allocatable :: elemPos(:)
    logical :: isConverged
    ! ---------------------------------------------------------------------------
    ! number of variables in varMap
    nVars = me%varMap%varPos%nVals

    ! Number of scalars in current output
    nScalars = me%varMap%nScalars

    if (me%subTree%useGlobalMesh) then
      nElems = tree%nElems
    else
      nElems = me%subTree%nElems
    end if

    ! open spatial reduction
    call tem_reduction_spatial_open( me     = me%redSpatial,               &
      &                              varSys = varSys,                      &
      &                              varPos = me%varMap%varPos%val(:nVars) )

    ! allocate elemPos to size of chunkSize
    allocate(elemPos(me%chunkSize))

    ! Process all chunks to derive the quantities defined in the tracking object
    do iChunk = 1, me%nChunks
      ! Number of elements read so far in previous chunks.
      elemOff = ((iChunk-1)*me%chunkSize)

      ! number of elements written to THIS chunk
      nChunkElems = min(me%chunkSize, nElems-elemOff)

      ! Compute the element lower and upper bound for the current chunk
      buf_start = elemOff + 1
      buf_end = elemOff + nChunkElems

      if (me%subTree%useGlobalMesh) then
        elemPos(1:nChunkElems) = (/ (iElem, iElem=buf_start, buf_end) /)
      else
        elemPos(1:nChunkElems) = me%subTree                     &
          &                        %map2Global(buf_start:buf_end)
      end if

      ! evaluate all variables on current chunk
      call tem_get_element_chunk(varSys  = varSys,                 &
        &                        varPos  = me%varMap%varPos        &
        &                                           %val(:nVars),  &
        &                        elemPos = elemPos(1:nChunkElems), &
        &                        time    = time,                   &
        &                        tree    = tree,                   &
        &                        nElems  = nChunkElems,            &
        &                        nDofs   = me%nDofs,               &
        &                        res     = res                     )

      ! preform spatial reduction
      call tem_reduction_spatial_append( me     = me%redSpatial,            &
        &                                chunk  = res,                      &
        &                                nElems = nChunkElems,              &
        &                                treeID = tree%treeID(              &
        &                                         elemPos(1:nChunkElems) ), &
        &                                tree   = tree,                     &
        &                                varSys = varSys,                   &
        &                                nDofs  = me%nDofs,                 &
        &                                varPos = me%varMap%varPos          &
        &                                                  %val(:nVars)     )

    end do ! iChunk

    deallocate(elemPos)

    ! Close reduction for current convergence
    call tem_reduction_spatial_close( me   = me%redSpatial, &
      &                               proc = me%proc        )

    ! evaluate residual and check convergence for each scalar
    ! and set steady state flag only in root process of convergence
    ! communicator
    if (me%proc%rank == 0) then
      ! Now check for convergence for each reduced scalars
      call tem_convergence_evaluate( me       = me,         &
        &                            achieved = isConverged )

      ! set steady state if any convergence object is converged
      status%bits(tem_stat_steady_state) = status%bits(tem_stat_steady_state)&
        &                                  .or. isConverged
      if (isConverged) then
        write(logUnit(5),*) 'Reached steady state ', time%iter, isConverged
      end if
    end if

  end subroutine tem_convergence_check_element
! **************************************************************************** !


! **************************************************************************** !
  !> This routine runs over convergence check using get_point interface
  subroutine tem_convergence_check_point( me, time, status, varSys, tree, res )
    ! ---------------------------------------------------------------------------
    !> convergence descriptions
    type(tem_convergence_type), target, intent(inout)   :: me
    !> current simulation time
    type(tem_time_type), intent(in)     :: time
    !> Status bits
    type(tem_status_type), intent(inout) :: status
    !> global variable system
    type(tem_varSys_type), intent(in)         :: varSys
    !> global tree
    type(treelmesh_type ), intent(in)         :: tree
    !> Output data
    !! size: io_buffer_size
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------------
    integer :: nVars, nPoints, nScalars, pointsOff, nChunkPoints
    integer :: iPoint, iChunk, counter
    integer :: buf_start, buf_end
    real(kind=rk), allocatable :: points(:,:)
    logical :: isConverged
    ! ---------------------------------------------------------------------------
    ! number of variables in varMap
    nVars = me%varMap%varPos%nVals

    ! Number of scalars in current output
    nScalars = me%varMap%nScalars

    if (me%subTree%useGlobalMesh) then
      nPoints = tree%nElems
    else
      nPoints = me%subTree%nPoints
    end if

    ! open spatial reduction
    call tem_reduction_spatial_open( me     = me%redSpatial,               &
      &                              varSys = varSys,                      &
      &                              varPos = me%varMap%varPos%val(:nVars) )

    ! allocate points to size of chunkSize
    allocate(points(me%chunkSize,3))

    ! Process all chunks to derive the quantities defined in the tracking object
    do iChunk = 1, me%nChunks
      ! Number of points read so far in previous chunks.
      pointsOff = ((iChunk-1)*me%chunkSize)

      ! number of points written to THIS chunk
      nChunkPoints = min(me%chunkSize, nPoints-pointsOff)

      ! Compute the points lower and upper bound for the current chunk
      buf_start = pointsOff + 1
      buf_end = pointsOff + nChunkPoints

      if (me%subTree%useGlobalMesh) then
        counter = 0
        do iPoint = buf_start, buf_end
          counter = counter + 1
          points(counter, :) = tem_BaryOfId( tree,               &
            &                                tree%treeID(iPoint) )
        end do
      else
        points(1:nChunkPoints,:) = me%subTree%points(buf_start:buf_end,:)
      end if

      ! evaluate all variables on current chunk
      call tem_get_point_chunk(varSys  = varSys,                     &
        &                      varPos  = me%varMap%varPos          &
        &                                         %val(:nVars),    &
        &                      point   = points(1:nChunkPoints,:), &
        &                      time    = time,                     &
        &                      tree    = tree,                     &
        &                      nPnts   = nChunkPoints,             &
        &                      res     = res                       )

      ! preform spatial reduction
      call tem_reduction_spatial_append( me     = me%redSpatial,        &
        &                                chunk  = res,                  &
        &                                nElems = nChunkPoints,         &
        &                                tree   = tree,                 &
        &                                varSys = varSys,               &
        &                                varPos = me%varMap%varPos      &
        &                                                  %val(:nVars) )

    end do ! iChunk

    deallocate(points)

    ! Close reduction for current convergence
    call tem_reduction_spatial_close( me   = me%redSpatial, &
      &                               proc = me%proc        )

    ! evaluate residual and check convergence for each scalar
    ! and set steady state flag only in root process of convergence
    ! communicator
    if (me%proc%rank == 0) then
      ! Now check for convergence for each reduced scalars
      call tem_convergence_evaluate( me       = me,         &
        &                            achieved = isConverged )

      ! set steady state if any convergence object is converged
      status%bits(tem_stat_steady_state) = status%bits(tem_stat_steady_state)&
        &                                  .or. isConverged
      if (isConverged) then
        write(logUnit(10),*) 'Reached steady state ', time%iter, isConverged
      end if
    end if

  end subroutine tem_convergence_check_point
! **************************************************************************** !


! **************************************************************************** !
  !> Evaluate if the convergence was achieved
  !!
  subroutine tem_convergence_evaluate( me, achieved )
    ! ---------------------------------------------------------------------------
    !> Convergence description contians
    type(tem_convergence_type), intent(inout) :: me
    !> is all Scalars in current convergence_type are converged
    logical, intent(out)                  :: achieved
    ! ---------------------------------------------------------------------------
    integer :: iVar, iComp, iScalar
    real(kind=rk) :: residual, threshold_fac
    logical :: isConverged(me%varMap%nScalars)
    ! ---------------------------------------------------------------------------
    ! Increase the counter for the checks
    me%nChecks = me%nChecks + 1
    iScalar = 0
    do iVar = 1, me%varMap%varPos%nVals
      do iComp = 1, me%redSpatial(iVar)%nComponents
        iScalar = iScalar + 1
        ! Compare the results against threshold
        !norm = abs(state(1) - me%result_prev)
        residual = evaluate_residual( me      = me,                            &
          &                           state   = me%redspatial(ivar)%val(icomp),&
          &                           iScalar = iScalar                        )

        if( me%header%absoluteError ) then
          ! For absolute error, just use threshold
          threshold_fac = me%header%cond(iVar)%threshold
        else
          ! For relative errors, multiply threshold with current state value
          threshold_fac = me%redSpatial(iVar)%val(iComp) &
            &           * me%header%cond(iVar)%threshold
        end if

        isConverged(iScalar) = tem_comparator(                                 &
          &                        val       = residual,                       &
          &                        operation = me%header%cond(iVar)%operation, &
          &                        threshold = threshold_fac                   )
      end do
    end do
    ! update the convergence achieved if the last compare was succesful
    achieved = all(isConverged)

  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, iScalar ) result( res )
    ! ---------------------------------------------------------------------------
    !> Convergence description
    type(tem_convergence_type), intent(inout) :: me
    !> Spatially reduced variable value
    real(kind=rk), intent(in)              :: state
    !> Current scalar
    integer, intent(in) :: iScalar
    !> residual to check for convergence
    real(kind=rk) :: res
    ! ---------------------------------------------------------------------------
    integer :: pos_lastState
    real(kind=rk) :: average
    ! ---------------------------------------------------------------------------
    ! Reset the result
    res = huge( res )

    select case (me%norm_kind )
    case( norm_simple )
      if (me%nChecks > 1) then
        res = abs((state - me%lastState(1, iScalar)))
      end if
      ! update the result at t-1 to t as when we arrive at t+1, it will
      ! be required
      me%lastState(1, iScalar) = state
    case( norm_average )
      pos_lastState = mod( me%nChecks - 1, me%header%nLastVals ) + 1
      if ( me%nChecks <= me%header%nLastVals ) then
        average = sum( me%lastState(1:pos_lastState, iScalar) ) &
          &     / real( pos_lastState, kind=rk )
      else
        average = sum( me%lastState(:, iScalar) )    &
          &     / real( me%header%nLastVals, kind=rk )
      end if

      res = abs( (state - average ) )
      me%lastState( pos_lastState, iScalar ) = state
      !write(*,*) 'nCheck ', me%nChecks, iScalar
      !write(*,*) 'lastState', me%lastState
      !write(*,*) 'state:', state,'average', average, 'res', res, &
      !  &        'pos last', pos_lastState
    end select

  end function evaluate_residual
! **************************************************************************** !


! **************************************************************************** !
  !> Dumps array of convergence to given unit
  subroutine tem_convergence_dump_vector(me, outUnit)
    ! ---------------------------------------------------------------------------
    !> convergence to write into the lua file
    type(tem_convergence_type), intent(in) :: me(:)
    !> unit to write to
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------
    ! aotus type handling the output to the file in lua format
    type(aot_out_type) :: conf
    ! ---------------------------------------------------------------------------
    call aot_out_open( put_conf = conf, outUnit = outUnit )
    call tem_convergence_out_vector( me, conf )
    call aot_out_close( put_conf = conf )

  end subroutine tem_convergence_dump_vector
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Dump single convergence to given unit
  subroutine tem_convergence_dump_single(me, outUnit)
    ! ---------------------------------------------------------------------------
    !> convergence to write into the lua file
    type(tem_convergence_type), intent(in) :: me
    !> unit to write to
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------
    ! aotus type handling the output to the file in lua format
    type(aot_out_type) :: conf
    ! ---------------------------------------------------------------------------
    call aot_out_open( put_conf = conf, outUnit = outUnit )
    call tem_convergence_out_single( me, conf )
    call aot_out_close( put_conf = conf )

  end subroutine tem_convergence_dump_single
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Allows the output of array of convergence to lua out
  subroutine tem_convergence_out_vector(me, conf)
    ! ---------------------------------------------------------------------------
    !> convergence to write into the lua file
    type(tem_convergence_type), intent(in) :: me(:)
    !> aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf
    ! ---------------------------------------------------------------------------
    integer :: iConv
    ! ---------------------------------------------------------------------------
    call aot_out_open_table( put_conf = conf, tname='convergence' )
    do iConv = 1,size(me)
      call tem_convergence_out_single( me(iConv), conf, level=1 )
    end do
    call aot_out_close_table( put_conf = conf )

  end subroutine tem_convergence_out_vector
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Allows the output of the single convergence to lua out.
  !!
  !! The data is written into the file, the lunit is connected to.
  !! It is formatted as a Lua table.
  !!
  subroutine tem_convergence_out_single(me, conf, level)
    ! ---------------------------------------------------------------------------
    !> convergence to write into the lua file
    type(tem_convergence_type), intent(in) :: me
    !> aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf
    !> to dump variable with key or without key
    integer, optional, intent(in) :: level
    ! ---------------------------------------------------------------------------
    integer :: level_loc, iVar
    ! ---------------------------------------------------------------------------
    if (present(level)) then
      level_loc = level
    else
      level_loc = 0
    end if

    if( level_loc == 0) then
      call aot_out_open_table( put_conf = conf, tname = 'convergence' )
    else
      call aot_out_open_table( put_conf = conf )
    end if

    ! write variable names
    call aot_out_open_table( put_conf = conf, tname = 'variable' )
    do iVar = 1, me%header%nRequestedVars
      call aot_out_val( put_conf = conf,                         &
        &               val      = trim(me%header%varName(iVar)) )
    end do
    call aot_out_close_table( put_conf = conf )

    ! convergence norm
    call aot_out_val( put_conf = conf,           &
      &               val      = me%header%norm, &
      &               vname    = 'norm'          )
    call aot_out_val( put_conf = conf,                &
      &               val      = me%header%nLastVals, &
      &               vname    = 'nvals'              )

    call aot_out_val( put_conf = conf,                    &
      &               val      = me%header%absoluteError, &
      &               vname    = 'absolute'               )

    ! write reductions
    call tem_reduction_spatial_out( me%redSpatial, conf )

    ! write conditions
    call tem_condition_out( me%header%cond, conf )

    ! write timeControl info
    call tem_timeControl_out( me%header%timeControl, conf )
    ! write shapes
    call tem_shape_out( me%header%geometry, conf )

    call aot_out_close_table( put_conf = conf )

  end subroutine tem_convergence_out_single
  ! ************************************************************************** !


! **************************************************************************** !
  !> This function provides the assignment operator of two convergence.
  !! Temporary Solution as CRAY compiler dont have = Operator
  !! Copying a convegence object (right) into other convergence (left)
  !!
  subroutine Copy_convergence(left,right)
    ! ---------------------------------------------------------------------------
    !> convegence to copy to
    type(tem_convergence_type), intent(out) :: left
    !> convegence to copy from
    type(tem_convergence_type), intent(in) :: right
    ! ---------------------------------------------------------------------------

    left%header = right%header
    left%norm_kind = right%norm_kind
    left%varMap = right%varMap
    left%subTree = right%subTree
    left%proc = right%proc
    if (allocated(right%redSpatial)) then
      allocate(left%redSpatial(size(right%redSpatial)))
      left%redSpatial = right%redSpatial
    end if
    left%nDofs = right%nDofs

  end subroutine Copy_convergence
! **************************************************************************** !


end module tem_convergence_module
! **************************************************************************** !
