! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> A tracking output (e.g. for macroscopic values) for tracking certain flow
!! Quantities over a period of time at certain locations in the grid.
!!
!! See the \ref tracking_mod "detailed tracking module description".
!!
!! \author Manuel Hasert, Simon Zimny, Nikhil Anand
!!
module tem_tracking_module

  ! incude treelm modules
  use env_module,                     only: rk, labelLen, pathLen, pathSep, &
    &                                       io_buffer_size
  use treelmesh_module,               only: treelmesh_type
  use tem_bc_prop_module,             only: tem_bc_prop_type
  use tem_aux_module,                 only: tem_abort
  use tem_comm_env_module,            only: tem_comm_env_type
  use tem_reduction_spatial_module,   only: tem_load_reduction_spatial,        &
    &                                       tem_reduction_spatial_config_type, &
    &                                       tem_reduction_spatial_init
  use tem_shape_module,               only: tem_shape_type, tem_load_shape
  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_time_module,                only: tem_time_type
  use tem_timeControl_module,         only: tem_timeControl_type,              &
    &                                       tem_timeControl_load,              &
    &                                       tem_timeControl_dump,              &
    &                                       tem_timeControl_check
  use tem_varSys_module,              only: tem_varSys_type
  use tem_varMap_module,              only: tem_varMap_type, tem_create_varMap
  use tem_tools_module,               only: tem_horizontalSpacer
  use tem_logging_module,             only: logUnit
  use tem_stencil_module,             only: tem_stencilHeader_type
  use tem_simControl_module,          only: tem_simControl_type
  use tem_status_module,              only: tem_stat_steady_state,             &
    &                                       tem_stat_stop_file,                &
    &                                       tem_status_run_terminate
  use tem_debug_module,               only: dbgUnit

  ! 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

  ! include libharvesting modules
  use hvs_output_module,   only: hvs_output_config_type, hvs_output_file_type, &
    &                            hvs_output_load, hvs_output_open,             &
    &                            hvs_output_write, hvs_output_close,           &
    &                            hvs_output_finalize, hvs_AsciiTransient,      &
    &                            hvs_output_init, hvs_VTK

  implicit none

  private

  public :: tem_tracking_type
  public :: tem_tracking_instance_type
  public :: tem_tracking_config_type
  public :: tem_trackingControl_type
  public :: tem_load_tracking
  public :: tem_init_tracker
  public :: tem_init_tracker_subTree
  public :: tem_tracking_finalize
  public :: tem_tracking_getData
  public :: tem_tracking_has_triggered
  public :: tem_tracker
  public :: tem_tracking_print_last_VTK_files

  !> General information about the tracking entities
  !! This data type is set in tem_load_tracking,
  !! then updated in tem_init_tracker_subTree
  !! After load balancing, it is reset in tem_reload_tracking
  type tem_trackingControl_type
    !> output status, activated?
    logical :: active = .false.
    !> number of tracking entities active on this process
    integer :: nActive = 0
    !> Total number of tracking entities defined in config file
    integer :: nDefined = 0
  end type tem_trackingControl_type

  !> Contains all config information about tracking
  !! Content in tracking config must NOT change!
  type tem_tracking_config_type

    !> log object labels
    character(len=labelLen) :: label

    !> folder to store files to
    character(len=pathLen) :: prefix

    !> array of requested variable labels
    character(len=labelLen), allocatable :: varName(:)

    !> stores time control parameters
    type(tem_timeControl_type) :: timeControl

    !> tracking shapes
    type(tem_shape_type), allocatable  :: geometry(:)

    !> originally set to true. But if false the exact polynomial is evaluated at
    ! the point tracked
    logical :: track_complete_element

    !> Data loaded from output table
    type(hvs_output_config_type) :: output_config

    !> Spatial reduction config which is loaded from disk
    type( tem_reduction_spatial_config_type ) :: redSpatial_config

  end type tem_tracking_config_type

  !> Tracking entity definition
  type tem_tracking_instance_type
    !> Contains name and position of variables to track in global varSys
    !! number of found variables can be accessed by me%varMap%varPos%nVals
    type(tem_varMap_type) :: varMap

    !> sub-tree resulting from the elements within the tracking shape
    !! The sub-tree also holds the sub-communicator
    !! This data needs to be UPDATED after balance
    type(tem_subTree_type) :: subTree

    !> Description for output file formats
    type(hvs_output_file_type) :: output_file

    !> Pointer to config array in tem_tracking_type
    integer :: pntConfig
  end type tem_tracking_instance_type

  type tem_tracking_type
    !> General information about the tracking entities
    type(tem_trackingControl_type) :: control
    !> tracking header for collecting the properties from the lua file
    type(tem_tracking_config_type), allocatable :: config(:)
    !> Instances of tracking type active on this process
    type(tem_tracking_instance_type), allocatable :: instance(:)
  end type tem_tracking_type

!TG  interface assignment(=)
!TG    module procedure Copy_tracking
!TG  end interface

  contains

! ****************************************************************************** !
  !> Read the tracker configuration from the main lua file
  !!
  !! Setup the values for the tracking entities
  !!
  subroutine tem_load_tracking(me, conf, parent)
    !--------------------------------------------------------------------------*
    !> list of the trackingeentities to create
    type( tem_tracking_type ), intent(out) :: me
    !> handle of the lua config file
    type( flu_state ) :: conf
    !> if the tracking table is a child-table of some other table,
    !! use the parent as a reference
    integer, optional :: parent
    !---------------------------------------------------------------------------
    integer :: tc_handle, sub_handle
    integer :: iTrack, nTracks
    !---------------------------------------------------------------------------

    ! Read the number of trackings in the lua file
    call aot_table_open( L       = conf,       &
      &                  thandle = tc_handle,  &
      &                  key     = 'tracking', &
      &                  parent  = parent      )

    if (tc_handle == 0) then
      write(logUnit(1),*) 'No Tracking entities found!'
      call aot_table_close(L=conf, thandle=tc_handle)
      call tem_horizontalSpacer(fUnit=logUnit(1))
      me%control%nActive = 0
      me%control%nDefined = 0
      allocate( me%config(0) )
      allocate( me%instance(0) )
      return
    else ! track entity exists.
      me%control%active = .true.
    end if

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

    ! Only single table
    if (sub_handle == 0) then
      nTracks = 1
      write(logUnit(1),*) 'Tracking is a single table'
      allocate( me%config(1) )
      call tem_load_trackingConfig( conf       = conf,         &
        &                           sub_handle = tc_handle,    &
        &                           config     = me%config(1)  )
      call aot_table_close(L=conf, thandle=sub_handle)
    else ! Multiple table
      call aot_table_close(L=conf, thandle=sub_handle)
      nTracks = aot_table_length(L=conf, thandle=tc_handle)
      ! Allocate the defined number of tracking entities
      allocate( me%config( nTracks ))
      write(logUnit(1),"(A,I0)") 'Number of Tracking entities: ', nTracks

      ! Loop over all the definitions and assign the variables from the lua
      ! file on the tem_tracking_type.
      ! Inside this routine it will open tracking subtable. Each subtable
      ! contains one or more tracking variables the stuff is done in the
      ! routine tem_load_trackingConfig
      do iTrack = 1, nTracks
        write(logUnit(1),"(A,I0)") 'Loading tracker: ', iTrack
        call aot_table_open( L       = conf,       &
          &                  parent  = tc_handle,  &
          &                  thandle = sub_handle, &
          &                  pos     = iTrack      )
        call tem_load_trackingConfig( conf       = conf,             &
          &                           sub_handle = sub_handle,       &
          &                           config     = me%config(iTrack) )
        call aot_table_close(L=conf, thandle=sub_handle)
        write(logUnit(1),"(A,I0)") 'Done tracker ', iTrack
      end do
    end if ! sub_handle

    ! me%control%nActive = nTracks
    me%control%nDefined = nTracks
    allocate(me%instance(nTracks))

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

  end subroutine tem_load_tracking
! ****************************************************************************** !


! ****************************************************************************** !
  !> Read the tracker variables from tracking subtables defined in
  !! configuration from the main lua file
  !!
  !! If tracking is just a single table with single tracking entry
  !! then load only one tracking log exists with
  !! one or more variables using tem_load_trackingHeader_single.
  !! Else if tracking is table of many log then allocate log and load each
  !! log type using tem_load_trackingHeader_single
  !! Setup the values for the tracking entities
  !!
  subroutine tem_load_trackingConfig(config, conf, sub_handle)
    !--------------------------------------------------------------------------
    !> list of the tracking entities to create
    type( tem_tracking_config_type ), intent(out) :: config
    !> handle of the lua config file
    type( flu_state ) :: conf
    !> table sub-handle for the tracking table
    integer, intent(in) :: sub_handle
    !--------------------------------------------------------------------------
    integer :: iError            ! error flag handle
    integer, allocatable :: vError(:)
    !> number of requested variables
    integer :: nRequestedVars
    ! --------------------------------------------------------------------------

    call aot_get_val( L       = conf,                                          &
      &               thandle = sub_handle,                                    &
      &               val     = config%label,                                  &
      &               ErrCode = iError,                                        &
      &               key     = 'label',                                       &
      &               default = 'unnamed_track')

    write(logUnit(1),*) 'Tracking label: '//trim( config%label )

    call aot_get_val( val       = config%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 track in '//trim(config%label)
      call tem_abort()
    end if

    nRequestedVars = size(config%varName)

    ! load time control to output tracking
    call tem_timeControl_load( conf   = conf,              &
      &                        parent = sub_handle,        &
      &                        me     = config%timeControl )
    call tem_timeControl_dump(config%timeControl, logUnit(2))

    ! Where to store the tracking file?
    call aot_get_val( L       = conf,          &
      &               thandle = sub_handle,    &
      &               val     = config%prefix, &
      &               ErrCode = iError,        &
      &               key     = 'folder',      &
      &               default = '.'//pathSep   )

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

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

    ! Load SPATIAL reductions
    call tem_load_reduction_spatial(                                   &
      &                   conf              = conf,                    &
      &                   parent            = sub_handle,              &
      &                   redSpatial_config = config%redSpatial_config )

    if( config%redSpatial_config%active ) then
      ! Check if the number of reductions correspond to the number of variables
      ! in the system
      if( size( config%redSpatial_config%reduceType ) /= nRequestedVars ) then
        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
    end if

    ! Load output table for vis_kind
    call hvs_output_load( me       = config%output_config,          &
      &                   conf     = conf,                          &
      &                   parent   = sub_handle,                    &
      &                   isReduce = config%redSpatial_config%active )

  end subroutine tem_load_trackingConfig
! ****************************************************************************** !


! ****************************************************************************** !
  !> Routine creates subTree for each tracking object and removes tracking
  !! objects on process which do not include any elements to track
  !!
  !! Identify, how many and which elements exist on my local process and are
  !! requested from the trackers
  !! Empty tracking entities are removed, so the track(:) might be re-allocated
  subroutine tem_init_tracker_subTree( me, tree, solver, bc_prop, stencil, &
    &                                  prefix )
    !---------------------------------------------------------------------------
    !> tracking entities
    type(tem_tracking_type), intent(inout)             :: me
    !> Global mesh from which the elements are identified and then stored to
    !! sub-meshes inside the trackers
    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
    !> Global solver information
    type(tem_solveHead_type), intent(in)               :: solver
    !> stencil used to create subTree of boundary type
    type(tem_stencilHeader_type), optional, intent(in) :: stencil
    !> Prefix for output filename
    !! Usually: solver%simName
    character(len=labelLen), optional, intent(in)      :: prefix
    !---------------------------------------------------------------------------
    integer :: iLog, nActive
    ! temporary tracker array
    type( tem_tracking_instance_type ), allocatable :: tempTrack(:)
    ! prefix for tracking label
    character(len=pathLen) :: prefix_loc
    ! tracking%config%prefix//tracking%config%label
    character(len=pathLen) :: basename
    ! ---------------------------------------------------------------------------
    call tem_horizontalSpacer(fUnit=logUnit(1))
    write(logUnit(3),*) 'Initialize tracking subTree to remove empty objects'
    call tem_horizontalSpacer(fUnit=logUnit(1))

    nActive = 0

    if (present(prefix)) then
      prefix_loc = trim(prefix)
    else
      ! prefix for tracking label
      prefix_loc = trim(solver%simName)//'_'
    end if

    if( me%control%active ) then
      ! Allocate the temporary track
      allocate(tempTrack( me%control%nDefined ) )

      do iLog = 1, me%control%nDefined

        basename = trim(me%config(iLog)%prefix) // trim(prefix_loc) // &
          &        trim(me%config(iLog)%label)

        write(logUnit(3),*) 'Creating subTree for tracking object ' &
          &                 // trim( me%config(iLog)%label )

        !-----------------------------------------------------------------------
        ! identify tracker elements
        !-----------------------------------------------------------------------
        call tem_create_subTree_of( inTree    = tree,                       &
          &                         bc_prop   = bc_prop,                    &
          &                         stencil   = stencil,                    &
          &                         subTree   = me%instance(iLog)%subTree,  &
          &                         inShape   = me%config(iLog)%geometry,   &
          &                         storePnts = me%config(iLog)             &
          &                                     %output_config%useGetPoint, &
          &                         prefix    = trim(basename)              )

        ! get rid of the empty track in order to avoid empty writes to disk
        if ( me%instance(iLog)%subTree%useGlobalMesh .or. &
          &  ( me%instance(iLog)%subTree%nElems > 0 ) ) then
          nActive = nActive + 1
          tempTrack( nActive ) = me%instance(iLog)
          ! Pointer to array of tracking headers loaded from config file
          tempTrack( nActive )%pntConfig = iLog
        end if

      end do  ! nActive

      deallocate(me%instance)
      allocate( me%instance(nActive) )
      me%control%nActive = nActive

      do iLog = 1, nActive
        ! Copy the stuff from the temporary track
        me%instance(iLog) = temptrack(iLog)
      end do

      deallocate(temptrack)
    end if ! if tracking active

  end subroutine tem_init_tracker_subTree
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initialize the tracker entities:
  !! * create varMap, i.e. map requested variables to global variable system
  !! * initialize spatial reduction
  !! * initialize hvs output
  !!
  subroutine tem_init_tracker( me, tree, solver, varSys, nDofs, globProc, &
    &                          solSpec_unit                               )
    !---------------------------------------------------------------------------
    !> tracking entities
    type(tem_tracking_type),intent(inout) :: me
    !> Global mesh from which the elements are identified and then stored to
    !! sub-meshes inside the trackers
    type(treelmesh_type), intent(in)                  :: tree
    !> Global solver information
    type(tem_solveHead_type),intent(in)               :: solver
    !> solver-provided variable systems
    type(tem_varSys_type), intent(in)                 :: varSys
    !> The number of dofs for each scalar variable of the equation system
    integer, intent(in), optional                     :: nDofs
    !> Process description to use.
    type(tem_comm_env_type), intent(in)               :: globProc
    !> Solver specific unit for restart header
    integer, optional, intent(in)                  :: solSpec_unit
    !---------------------------------------------------------------------------
    integer :: iLog, nVars, iVar, iConfig
    ! prefix for tracking label to differiate tracking for different scheme
    ! with same tracking label
    character(len=pathLen) :: prefix
    ! tracking%config%prefix//tracking%config%label
    character(len=pathLen) :: basename
    !---------------------------------------------------------------------------

    call tem_horizontalSpacer(fUnit=logUnit(1))
    write(logUnit(1),*) 'Initialize tracking objects'
    call tem_horizontalSpacer(fUnit=logUnit(1))

    ! prefix for tracking label to differiate tracking for different scheme
    ! with same tracking label
    prefix = trim(solver%simName)//'_'

    if( me%control%active ) then

      do iLog = 1, me%control%nActive
        iConfig = me%instance(iLog)%pntConfig

        write(logUnit(3),"(A,I0,A)") 'Track object: ', iLog, ', label: ' &
          &                          // trim( me%config(iConfig)%label )

        ! map variables
        ! create tracking variable position in the global varSys
        call tem_create_varMap( varname = me%config(iConfig)%varname, &
          &                     varSys  = varSys,                     &
          &                     varMap  = me%instance(iLog)%varMap    )

        nVars = me%instance(iLog)%varMap%varPos%nVals
        ! Abort if none variables of the variable defined in current
        ! tracking object are found in varSys
        if ( nVars /= size( me%config(iConfig)%varname ) ) then
          write(logUnit(1),*) ' Some of the following variables are not found:'
          do iVar = 1, size(me%config(iConfig)%varName)
            write(logUnit(1),*) trim(me%config(iConfig)%varName(iVar))
          end do
          call tem_abort()
        end if

        basename = trim(me%config(iConfig)%prefix) // trim(prefix) &
          &        // trim(me%config(iConfig)%label)

        ! Init spatial reduction
        me%instance(iLog)%output_file%ascii%isReduce = me%config(iConfig)     &
          &                                            %redSpatial_config%active
        if ( me%config(iConfig)%redSpatial_config%active ) then
          ! Initialize reduction
          call tem_reduction_spatial_init(                                     &
            &                 me = me%instance(iLog)%output_file%ascii         &
            &                        %redSpatial,                              &
            &  redSpatial_config = me%config(iConfig)%redSpatial_config,       &
            &             varSys = varSys,                                     &
            &             varPos = me%instance(iLog)%varMap%varPos%val(:nVars) )
        end if

        ! Initialize output
        if ( me%instance(iLog)%subTree%useGlobalMesh ) then
          call hvs_output_init(out_file    = me%instance(iLog)%output_file,    &
            &                  out_config  = me%config(iConfig)%output_config, &
            &                  tree        = tree,                             &
            &                  varSys      = varSys,                           &
            &                  varPos      = me%instance(iLog)%varMap%varPos   &
            &                                                 %val(:nVars),    &
            &                  basename    = trim(basename),                   &
            &                  globProc    = globProc,                         &
            &                  timeControl = me%config(iConfig)%timeControl,   &
            &                  solver      = solver,                           &
            &                  geometry    = me%config(iConfig)%geometry,      &
            &                  nDofs       = nDofs,                            &
            &                  solSpec_unit = solSpec_unit                     )
        else
          call hvs_output_init(out_file    = me%instance(iLog)%output_file,    &
            &                  out_config  = me%config(iConfig)%output_config, &
            &                  tree        = tree,                             &
            &                  varSys      = varSys,                           &
            &                  varPos      = me%instance(iLog)%varMap%varPos   &
            &                                                 %val(:nVars),    &
            &                  subTree     = me%instance(iLog)%subTree,        &
            &                  basename    = trim(basename),                   &
            &                  globProc    = globProc,                         &
            &                  timeControl = me%config(iConfig)%timeControl,   &
            &                  solver      = solver,                           &
            &                  geometry    = me%config(iConfig)%geometry,      &
            &                  nDofs       = nDofs,                            &
            &                  solspec_unit = solSpec_unit                     )
        end if
      end do

    end if ! if tracking active

  end subroutine tem_init_tracker
! ****************************************************************************** !


! ****************************************************************************** !
  !> A routine to obtain tracked data.
  !!
  !! This routine will return all requested variables in the tracking object
  !! me and return it for all elements of the subtree in the res field.
  subroutine tem_tracking_getData(varMap, subTree, varSys, mesh, time, nDofs, res)
    !> varMap from tem_tracking_instance_type
    type(tem_varMap_type) :: varMap

    !> subTree from tem_tracking_instance_type
    type(tem_subTree_type) :: subTree

    !> Variable system describing available data.
    type(tem_varsys_type), intent(in) :: varsys

    !> Mesh definition of the input data.
    type(treelmesh_type), intent(in) :: mesh

    !> Time information for the current data.
    type(tem_time_type), intent(in) :: time

    !> Number of degrees of freedom.
    integer, intent(in) :: nDofs

    !> Tracked data, has to match the subtree definition.
    !!
    !! The memory layout is like this:
    !!  1. All variable components
    !!  2. nDofs
    !!  3. nElems (subtree%nElems)
    real(kind=rk), intent(out) :: res(:)
    !------------------------------------------------------------------------!
    integer :: maxComponents
    integer :: nComponents
    integer :: compOff
    integer :: elemOff
    integer :: nElems
    integer :: nScalars, nVars
    integer :: elemSize
    integer :: nChunks
    integer :: chunksize
    integer :: nChunkElems
    integer :: res_size
    integer :: buf_start, buf_end
    integer :: e_start, d_start, t_start
    integer :: iElem, iChunk, iDoF, iVar
    integer :: varpos
    real(kind=rk), allocatable :: tmpdat(:)
    !------------------------------------------------------------------------!

    nElems = subTree%nElems
    nScalars = varMap%nScalars
    nVars = varMap%varpos%nVals

    ! Need to obtain the data variable for variable, and store it in an
    ! intermediate array, because all components should be put together in the
    ! res array.
    ! The temporary array therefore needs to be sufficiently large to store the
    ! maximal number of components.
    maxComponents = maxval(varSys%method%val(varMap%varPos &
      &                          %val(:nVars))%nComponents )

    ! Number of elements to fit into a single chunk.
    chunkSize = min( io_buffer_size / (maxComponents*nDofs), nElems )

    ! Size of a single element
    elemsize = nScalars*nDofs

    if ( (nElems > 0) .and. (chunkSize == 0) ) then
      write(logUnit(0),*) 'The chosen io_buffer_size of ', io_buffer_size
      write(logUnit(0),*) 'is too small for outputting ', maxComponents
      write(logUnit(0),*) 'scalar values with ', nDofs
      write(logUnit(0),*) 'degrees of freedom!'
      write(logUnit(0),*) 'Please increase the io_buffer_size to at least ', &
        &                 real(maxComponents*nDofs) / real(131072), ' MB!'
      call tem_abort()
    end if

    ! Using a temporary array to store the variables and transfer them to res
    ! in the correct ordering afterwards.
    allocate(tmpdat(chunkSize*maxComponents*nDofs))

    nChunks = 0
    if (chunkSize > 0) then
      nChunks = ceiling( real(nElems, kind=rk)      &
        &                / real(chunkSize, kind=rk) )
    end if

    chunks: do iChunk=1,nChunks
      elemOff = ((iChunk-1)*chunkSize)
      nChunkElems = min(chunkSize, nElems - elemOff)
      buf_start = elemOff + 1
      buf_end = elemOff + nChunkElems

      compOff = 0
      vars: do iVar=1,varMap%varPos%nVals
        varpos = varMap%varPos%val(iVar)
        nComponents = varSys%method%val(varPos)%nComponents
        res_size = nChunkElems * nDofs * nComponents
        ! derive the quantities for all the elements in the current chunk
        call varSys%method%val(varpos)%get_element(                          &
          &                                varSys  = varSys,                 &
          &                                elemPos = subtree%map2global(     &
          &                                              buf_start:buf_end), &
          &                                time    = time,                   &
          &                                tree    = mesh,                   &
          &                                nElems  = nChunkElems,            &
          &                                nDofs   = nDofs,                  &
          &                                res     = tmpdat(:res_size)       )
        do iElem=1,nChunkElems
          e_start = (elemOff+iElem-1)*elemsize
          t_start = (iElem-1)*nComponents*nDofs
          do iDof=1,nDofs
            d_start = (iDof-1)*nScalars + compOff
            res( (e_start+d_start+1) : (e_start+d_start+nComponents) ) &
              &  = tmpdat( t_start + (iDof-1)*nComponents + 1 &
              &            :t_start + iDof*nComponents        )
          end do
        end do
        ! Increase the component offset for the next variables.
        compOff = compOff + nComponents
      end do vars
    end do chunks

  end subroutine tem_tracking_getData
! ***************************************************************************** !


! ***************************************************************************** !
  !> Decision on whether the giving tracker should be written in the current
  !! iteration.
  function tem_tracking_has_triggered( timeControl, simControl, proc ) &
    &      result(triggered)
    type(tem_timeControl_type), intent(inout) :: timeControl
    type(tem_comm_env_type),    intent(inout) :: proc
    type(tem_simControl_type),  intent(in) :: simControl
    logical :: triggered

    logical :: tc_triggered

    call tem_timeControl_check( me        = timeControl,    &
      &                         now       = simControl%now, &
      &                         comm      = proc%comm,      &
      &                         triggered = tc_triggered    )

    triggered = tc_triggered                                       &
      &         .or. simControl%status%bits(tem_stat_steady_state) &
      &         .or. simControl%status%bits(tem_stat_stop_file)    &
      &         .or. tem_status_run_terminate(simControl%status)
  end function tem_tracking_has_triggered

! ***************************************************************************** !


! ****************************************************************************** !
  !> This routine runs over each tracking object and dump requested quantities
  !! if timeControl is active on current time
  subroutine tem_tracker( track, simControl, varSys, tree )
    ! ---------------------------------------------------------------------------
    !> tracking object containing all tracking relevant information
    type(tem_tracking_type ), intent(inout)   :: track
    !> Simulation control contains current simulation time and
    !! status bits
    type(tem_simControl_type), intent(in)     :: simControl
    !> global variable system
    type(tem_varSys_type), intent(in)         :: varSys
    !> global tree
    type(treelmesh_type ), intent(in)         :: tree
    ! ---------------------------------------------------------------------------
    integer :: iLog, iConfig
    ! ---------------------------------------------------------------------------

    ! Run over all tracking objects
    do iLog = 1, track%control%nActive
      iConfig = track%instance(iLog)%pntConfig

      ! Skip this tracking object, if there are no entries in the
      ! variable system
      if (track%instance( iLog )%varMap%nScalars < 1) cycle

      ! dump tracking when at least one of following conditions is triggered:
      !   tracking timeControl is triggered
      !   simulation reached steady state
      !   stop file is defined
      !   simulation is terminated abruptly
      if ( tem_tracking_has_triggered(                           &
        &    timeControl = track%config(iConfig)%timeControl,    &
        &    simControl  = simControl,                           &
        &    proc        = track%instance(iLog)%output_file%proc ) ) then

        if( track%instance( iLog )%subTree%useGlobalMesh ) then
          ! Open the output files, this also generates the vertices for the mesh,
          ! and writes the mesh data to disk. Also writes header file depends
          ! on output vis_kind
          call hvs_output_open(                                 &
            &    out_file = track%instance(iLog)%output_file,   &
            &    use_iter = track%config(iConfig)%output_config &
            &                    %vtk%iter_filename,            &
            &    mesh     = tree,                               &
            &    varSys   = varSys,                             &
            &    time     = simControl%now                      )

          ! Evaluate and write results to disk
          call hvs_output_write( out_file = track%instance(iLog)%output_file, &
            &                    varSys   = varSys,                           &
            &                    mesh     = tree                              )

          ! Close opened files
          call hvs_output_close( out_file = track%instance(iLog)%output_file, &
            &                    varSys   = varSys,                           &
            &                    mesh     = tree                              )
        else
          ! Open the output files, this also generates the vertices for the mesh,
          ! and writes the mesh data to disk. Also writes header file depends
          ! on output vis_kind
          call hvs_output_open(                                 &
            &    out_file = track%instance(iLog)%output_file,   &
            &    use_iter = track%config(iConfig)%output_config &
            &                    %vtk%iter_filename,            &
            &    mesh     = tree,                               &
            &    varSys   = varSys,                             &
            &    subTree  = track%instance(iLog)%subTree,       &
            &    time     = simControl%now                      )

          ! Evaluate and write results to disk
          call hvs_output_write( out_file = track%instance(iLog)%output_file, &
            &                    varSys   = varSys,                           &
            &                    mesh     = tree,                             &
            &                    subTree  = track%instance(iLog)%subTree      )

          ! Close opened files
          call hvs_output_close( out_file = track%instance(iLog)%output_file, &
            &                    varSys   = varSys,                           &
            &                    mesh     = tree,                             &
            &                    subTree  = track%instance(iLog)%subTree      )
        end if !Global mesh

      end if  ! do tracking? interval, tmin, tmax check
    end do ! iLog

  end subroutine tem_tracker
! ****************************************************************************** !


! ****************************************************************************** !
  !> Close all the units of the tracking objects
  !!
  subroutine tem_tracking_finalize( me )
    ! ---------------------------------------------------------------------------
    !> tracker object to close
    type(tem_tracking_type ), intent(inout) :: me
    ! ---------------------------------------------------------------------------
    integer :: iTrack !, iError
    ! ---------------------------------------------------------------------------

    do iTrack = 1, me%control%nActive
      call hvs_output_finalize( out_file = me%instance(iTrack)%output_file )
    end do

  end subroutine tem_tracking_finalize
! ****************************************************************************** !


!TG ! ****************************************************************************** !
!TG   !> This function provides the assignment operator of two tracking.
!TG   !! Temporary Solution as CRAY compiler dont have = Operator
!TG   !! Copying a tracking object (right) into other tracking (left)
!TG   !!
!TG   subroutine Copy_tracking(left,right)
!TG     ! ---------------------------------------------------------------------------
!TG     !> tracking to copy to
!TG     type(tem_tracking_type), intent(out) :: left
!TG     !> tracking to copy from
!TG     type(tem_tracking_type), intent(in) :: right
!TG     ! ---------------------------------------------------------------------------
!TG
!TG     left%config = right%config
!TG     left%instance%varMap = right%instance%varMap
!TG     left%instance%subTree = right%instance%subTree
!TG     left%instance%output_file = right%instance%output_file
!TG
!TG   end subroutine Copy_tracking
!TG ! ****************************************************************************** !

  ! ************************************************************************** !
  !> Print the filenames of the last written VTK files on the screen.
  !!
  !! Mainly useful for the harvesting tools to get information on the written
  !! visualization files.
  subroutine tem_tracking_print_last_VTK_files(track)
    !> List of tracking entities
    type(tem_tracking_type), intent(in) :: track

    integer :: iTrack, iConfig
    integer :: nTracks

    nTracks = track%control%nActive

    do iTrack = 1, nTracks
      iConfig = track%instance(iTrack)%pntConfig
      ! Only the root process for this output file needs to print the name
      if (track%instance(iTrack)%output_file%proc%rank == 0) then
        ! Only the VTK outputs are considered here
        if (track%config(iConfig)%output_config%vis_kind == hvs_VTK) then
          write(*,*) 'Wrote VTK file: ', &
            &        trim(track%instance(iTrack)%output_file%vtk%last_opened_file)
        end if
      end if
    end do

  end subroutine tem_tracking_print_last_VTK_files

end module tem_tracking_module
! ****************************************************************************** !
!> \page tracking_mod Tracking Module Description
!! \section tracking Tracking
!!
!! During the runtime of simulations, the user needs to keep track of the
!! simultation state without being obliged to write out the complete domain
!! to disk.
!! This is achieved by tracking objects, at which the state is evaluated
!! at certain defined intervals and then written to disk.
!!
!! \subsection tracking_format Tracking format
!!
!! We support three formats for tracking.
!! <em>Note: We strongly encourage you to use the Harvester format. </em>
!! Herefore you need to run Harvester on the simulation folder to create
!! files for the Visualization tools.
!! - \ref tem_tracking_dumpasciitransient "ascii" -
!! Transient data: simple gnuplot-style
!! A simple ascii file which creates a small header and then writes at each
!! defined time step a column wise representation of all variables and points.
!! This is not efficient and should only be used for simple point-tracking
!! objects.
!! Each process writes its own file.
!! - \ref tem_tracking_dumpasciispatial "asciiSpatial" -
!! Spatial Data: This format stores solution data along with coordinates in
!! single data file for every timestep. Therefore, there would be multiple
!! files depending on the total number of time-steps.
!! It is ideal for tracking *line* and *plane*.
!!
!! - \ref tem_restart_module::restart_format "Harvester Format"
!! Writes all participating elements to a treelm mesh and creates for each
!! required time step a restart file.
!! This is a scalable tracking object which works regardless of number of
!! processes.
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! format='harvester', folder = './tracking/'
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! <em>Note: We will demonstrate each format in the examples of tracking point,
!! line and plane in the following sections</em>
!!
!! \subsection tracking_variabletypes Tracking variable types
!!
!! - \ref tem_variable_module "variables":
!!  all variables defined in the \ref tem_varSys_module "global variable
!!                   system" can be tracked.
!!                   See \ref tem_variable_module::tem_variable_from
!!                   "how to define variables"
!! - derived variables: all variables defined in the global variable system
!!    can be tracked
!! - spacetime: as analytical result, can be defined in lua file
!! - \ref tem_reduction_spatial_module "reductions":
!!     e.g. difference, sum or average of defined variables
!!
!! \subsection tracking_shapes Tracking shapes
!!
!! We show the tracking with simple shapes here.
!! Get an overview over all existing shapes \ref tem_shape_module "here".
!!
!! - zero-dimensional: point
!! - one-dimensional: line
!! - two-dimensional: plane
!! - three-dimensional: cubes, cylinders
!!
!! So far point, line, plane are implemented in \ref tem_shape_module "shapes".
!!
!! \subsection tracking_usage Usage
!!
!! The user specifies the desired parameters in the configuration file
!! (in the scheme table) as
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! -- if you want to have a analytical solution as reference, you can add
!! --  a spacetime variable first
!! -- define single new (spacetime) variable
!! add_variable = { name='spacetime', ncomponents=1, spacetime=luaFun }
!! -- define multiple new (spacetime) variables
!! add_variable = {
!! {name='spacetime', ncomponents=1, spacetime=luaFun },
!! {name='spacetime2', ncomponents=1, spacetime=luaFun }
!! }
!!
!! -- Tracking
!! tracking = {
!! { label = 'probe_error',
!! variable = {
!!   { 'velocity' },
!!   { name='difference',
!!     ncomponents=1,
!!     dep={'density','spacetime'}},
!!   { name='difference',
!!     ncomponents=1,
!!     dep={'wss','spacetime2'}}},
!! folder = './tracking/',
!! format = 'ascii',
!! time = {interval = 1, min = 1000, max = 4000 },
!! reduction = {'l2norm'},
!! shape = { kind = 'canoND',
!!           object = { origin ={0.,0.,0.},
!!                      vec = {1.,0.,0.},
!!                      segments={100}} }
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!! Here, a point at the location {0.,0.,0.} is defined as a pressure probe.
!! The interval is set to 1, meaning at each time step, the pressure has to
!! be written.
!! The probing starts at simulation time 1000 and lasts until 4000.
!! The tracking write both pressure and velocityX in a single file.
!! Use ref_value to normalize the output.
!! \verbatim
!! # Name:       GaussPulse_TestSimulation     Date 2011-060-27  Time 11:43 CEST
!! # Probe type: Point at {0.0 0.0 0.0}
!! # Quantities: pressure, velocityX
!! # tMin:       1000,     tMax:    4000,     Interval: 1
!! # ref_value: 1.0, 0.01
!! #   iTime     rTime      Pressure           velocityX
!! 1000          1000.0     1.0241818          0.00158287
!! 1001          1001.0     1.0251427          0.00152617
!! ...
!! 4000          4000.0     1.7162511          0.01273673
!! \endverbatim
!!
!! The shape is defined as a canonical object.
!! The different shapes are basically derivatives of a simple cube.
!! A cube with no length in any direction is a point, a cube with length along
!! only one axis is a line whereas a cube with length along two axes is simply
!! treated as plane by the underlying code.
!! One example each for Point, Line and Plane is shown below.
!!
!! \subsection tracking_point Point Tracking
!!
!! <em>Example: tracking a point in ASCII format</em>
!!
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! tracking = { label = 'Velocity',
!! variable = {'velocity'},   -- options: density, velocity
!!
!! --   Tracking over a point
!! shape={kind = 'canoND', origin = {0.5,0.5,0.5} },
!!
!! time = {min = 0, interval = 1},     -- if no max given: active for all times
!! format = 'ascii',
!! folder = 'tracking/'}
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! Once, the above section is executed within the musubi.lua file, an output in
!! the asciiTransient format is generated in the output folder specified
!! (tracking in this particular example).
!! Now, we open the result file inside tracking folder, the structure of which
!! looks like below:
!!
!! \verbatim
!! # Rank of the process:       0
!! #               time     velocity_01     velocity_02     velocity_03
!! 1.0000  0.000000000000  0.000000000000  0.000000000000
!! 2.0000  0.000000000000  0.000000000000  0.000000000000
!! 3.0000  0.000000000000  0.000000000000  0.000000000000
!! 4.0000  0.000000000000  0.000000000000  0.000000000000
!! ...
!! ...
!! 1997.0000 -0.006065558329  0.000443400622 -0.000003481789
!! 1998.0000 -0.006063901533  0.000443795703 -0.000003286949
!! 1999.0000 -0.006061596302  0.000443327890 -0.000003534752
!! 2000.0000 -0.006059665212  0.000443720604 -0.000003339771
!! \endverbatim
!! Here, it can be noticed that 3 velocity components for all time steps at the
!! point (0.5, 0.5, 0.5) are stored. The interval was chosen as 1, so they are
!! tracked for each time step. If the interval was chosen as say 10, then the
!! tracking would have been done after 10 units of simulated time.
!!
!!
!! \subsection tracking_line Line Tracking
!!
!! We need to set the origin and the direction vector,
!! as e.g. {origin={0.,0.,0.}, vec={1.0,1.0,1.0}}.
!! The line is represented as a list of points in the tracking facility.
!! The number of subdivisions and its distribution of the points on the line
!! needs to be set, as e.g. @segments=500, distribution='equal'@.
!! Make sure to specify enough segments to take into account all the required
!! elements on the line.
!! Elements are added only once, so it might be beneficial to specify much more
!! segments than actually needed to consider elements on the highest levels as
!! well.
!!
!! <em>Example: Tracking a line in ASCII Spatial format</em>
!!
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! tracking = {
!!   label = 'line_probe',
!!   variable = {'state'}, -- state (=pdfs), pressure, density etc
!!   shape = { kind = 'canoND',
!!             object={
!!               origin = {0.5,0.5,0.5},
!!               vec = { {0.0,0.2,0.0}},
!!               segments = {300},
!!               distribution='equal'}
!!           },
!!   format='asciispatial',
!!   folder='tracking/',
!!   time = {interval = 5, min = 0, max = 800 } }
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! After, the above tracking is executed within musubi.lua file multiple res
!! files are obtained inside the tracking folder each specific to particular
!! time step. If you open any one of the 'res' files, there will be six data
!! columns: three coordinate and three velocity components.
!! An example is shown below:
!! \verbatim
!! # Rank of the process:       0
!! #         coordX          coordY          coordZ     velocity_01     velocity_02     velocity_03
!! 0.51562500      0.51562500      0.51562500 -0.000998448727 -0.000118182367  0.000007560959
!! 0.51562500      0.54687500      0.51562500 -0.000922746216  0.000011023243 -0.000000130462
!! 0.51562500      0.57812500      0.51562500 -0.001083956104  0.000038534363 -0.000003774649
!! ...
!! 0.51562500      0.64062500      0.51562500 -0.001343981219 -0.000131107478 -0.000004381864
!! 0.51562500      0.67187500      0.51562500 -0.001261956611  0.000081027056  0.000000215364
!! 0.51562500      0.70312500      0.51562500 -0.001733303068  0.000060154818  0.000005793207
!! \endverbatim
!!
!! Out of three coordinate components, only 'coordY' is varying and 'coordX',
!! 'coordZ' remain constant. This is according to our tracking input with
!! @origin = {0.5, 0.5, 0.5}, vec = {0.0, 0.2, 0.0}@.
!! You can now make a line plot with any one of velocity components with
!! respect to 'coordY' in "gnuplot":http://www.gnuplot.info/.
!!
!! \subsection tracking_plane Plane Tracking
!!
!! We need to set the origin and the direction vectors, as e.g.
!! object = {origin={0.,0.,0.}, vec={{1.0,1.0,1.0},{0.0,-1.0,1.0}}}.
!! The plane is, just as the line, represented as a list of points in the
!! tracking facility.
!! The number of subdivisions in each direction and its distribution of the
!! points on the plane needs to be set,
!! as e.g. @segments1=500, segments2=300, distribution='equal'@.
!! Make sure to specify enough segments to take into account all the required
!! elements on the line.
!! Elements are added only once, so it might be beneficial to specify much more
!! segments than actually needed to consider elements on the highest levels
!! as well.
!!
!! <em> Example: tracking a plane in harvester format</em>
!!
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! tracking = {
!!   label = 'plane_probe',
!!   variable = {'state'}, -- state (=pdfs), pressure, density etc
!!   shape = { kind = 'canoND',
!!             object={
!!               origin = {-75.,-15.,0.},
!!               vec = { {150.,0,0.}, {0.,30,0.}},
!!               segments = {300, 300},
!!               distribution='equal'}
!!           },
!!   format='harvester',
!!   folder='tracking/',
!!   time = {interval = 5, min = 0, max = 800 } }
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! \subsection tracking_global_mesh Global mesh tracking
!!
!! <em> Example: tracking the whole mesh</em>
!!
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! tracking = {
!!   ...
!!   shape = { kind = 'global' },
!!   ...
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! \subsection tracking_boundary Boundary elements tracking
!!
!! <em> Example: tracking the boundary elements</em>
!!
!! ~~~~~~~~~~~~~~~~~~~~~{.lua}
!! tracking = {
!!   ...
!!   shape = { kind = 'property', property = {'boundary'} },
!!   ...
!! }
!! ~~~~~~~~~~~~~~~~~~~~~
!!
!! \subsection tracking_design Design and Implementation
!!
!! The user-specified tracking shapes have to be translated to a set of
!! elements, on which the defined quantities are calculated and written to disk
!! in a suitable format.
!!
!! All elements are converted to a set of points.
!! For each point, the correspoinding treeID is added to the element list of
!! the tracker as \verbatim track%elems( iLevel )%elem(1:nElems)\endverbatim
!! Each element will only be added once.
!! All these elements are then treated in the main loop within the tracker
!! routine.
