1. Harald Klimach
  2. TreElM

Commits

Harald Klimach  committed e841bbc

Merged in further Treelm development
* Improved Multilevel handling
* Unit test for the copy operator
* More comments
* Convergence test

  • Participants
  • Parent commits 0366cc9
  • Branches default

Comments (0)

Files changed (27)

File source/arrayMacros.inc

View file
  • Ignore whitespace
     module procedure truncate_GA_?tname?
   end interface
 
+  !> Empty the entries  without changing arrays
+  interface empty
+    module procedure empty_GA_?tname?
+  end interface
+
   !> Destroy the dynamic array
   interface destroy
     module procedure destroy_GA_?tname?
 
   end subroutine truncate_GA_?tname?
 
+  subroutine empty_GA_?tname?(me)
+    !------------------------------------------------------------------------
+    type(grw_?tname?Array_type) :: me !< Array to sortTruncate
+    !------------------------------------------------------------------------
+    ?tstring?, allocatable :: tArray(:)
+    !------------------------------------------------------------------------
+
+    me%nVals = 0
+
+  end subroutine empty_GA_?tname?
+
 
   subroutine append_GA_?tname?(me, val, pos)
     type(grw_?tname?Array_type) :: me !< Array to append the value to
     module procedure truncate_DA_?tname?
   end interface
 
+  !> Truncate the array, meaning
+  !! cut off the trailing empty entries
+  interface empty
+    module procedure empty_DA_?tname?
+  end interface
+
+
   !> Fix the dynamic array, meaning:
   !! store the array in the sorted order and cut
   !! off the trailing empty entries
     integer :: i
     !------------------------------------------------------------------------
 
+    ! Do a binary search on existing entries (returns closest entry next to
+    ! it if not found).
+    FoundPos = SortedPosOfVal(me, val, .true.)
+    if( present( wasAdded ) ) wasAdded = .false.
 
-    FoundPos = SortedPosOfVal(me, val, .true.)
-    ! If val was not found in 'me', the FoundPos must point to the
-    ! new index at the end
-    if( present( wasAdded ) ) wasAdded = .false.
     ! If it found the value, the position is smaller than nVals
     if (FoundPos <= me%nVals) then
+
+      ! The returned position might actually be the right entry already or
+      ! not, check for this here.
       if ((me%val(me%sorted(FoundPos)) == val) &
         & .and. me%Unique) then
 
         ! Found the value in a list of unique values,
         ! nothing to do, just return its position.
         pos = me%sorted(FoundPos)
+
       else
+
         ! Need to append a new value!
+
         if( present( wasAdded ) ) wasAdded = .true.
         if (me%nVals == me%ContainerSize) then
 
         ! sorted list at the now freed position.
         me%sorted(FoundPos) = me%nVals
         pos = me%nVals
+
       end if
+
     else
+
       ! Value to append is larger than all existing ones,
       ! just put it to the end of the list, this captures
       ! also the case of empty lists.
       me%val(FoundPos) = val
       me%sorted(FoundPos) = FoundPos
       pos = FoundPos
+
     end if
+
   end subroutine append_DA_?tname?
 !******************************************************************************!
 
 
 
 !******************************************************************************!
+  !> Empty all contents of the array without changing the size or status of any
+  !! array
+  !!
+  subroutine empty_DA_?tname?(me)
+    !------------------------------------------------------------------------
+    type(dyn_?tname?Array_type) :: me !< Array to sortTruncate
+    !------------------------------------------------------------------------
+    type(dyn_?tname?Array_type) :: tArray !< temporary Array
+    integer :: iVal
+    integer :: dPos
+    !------------------------------------------------------------------------
+    ! Reset the number of entries 
+    me%nVals = 0
+  end subroutine empty_DA_?tname?
+!******************************************************************************!
+
+
+!******************************************************************************!
   !> Fixing the dynamic array
   !!
   !! Truncate the array after the last valid entry and hence cut off the empty
     !------------------------------------------------------------------------
     allocate(swpval(me%nVals))
     swpval = me%Val(1:me%nVals)
-    !if( allocated( me%val)) 
+    if( allocated( me%val)) &
     deallocate(me%Val)
 
     me%ContainerSize = max(minLength, me%ContainerSize*2)

File source/env_module.f90

View file
  • Ignore whitespace
   integer, parameter, public :: tem_err_warning = 1
 
   integer, parameter, public :: chatter_none = -1
+  integer, parameter, public :: chatter_normal = 1
   integer, parameter, public :: chatter_little = 3
   integer, parameter, public :: chatter_aLot = 10
   integer, parameter, public :: chatter_all = 99
     character(len=128) :: line
     integer :: stat
     integer :: inUnit
-    integer :: iTask
     integer :: outUnit
     !-------------------------------------------------------------------------!
     if( present( info )) then

File source/mus_derived/mus_derivedQuantitiesIncomp_module.f90

View file
  • Ignore whitespace
 
   private
 
+  public :: derivePressureIncomp_forElemFromState
   public :: deriveVelocityIncomp_forElemFromState
   public :: deriveVelMagIncomp_forElemFromState
   public :: deriveEquilibriumIncomp_forElemFromState
 !******************************************************************************!
 
 !******************************************************************************!
+  !> Calculate the pressure variation of a given set of elements for the 
+  !! incompressible Navier-Stokes equations
+  !!
+  !! Pressure variation calculation according to the isentropic equation of state for 
+  !! the LBM \f$ p = (\rho - \rho0) c_s^2 \f$
+  !! with the calculation of density as in deriveDensity_forElemFromState
+  !!
+  !! The interface has to comply to the abstract interface
+  !! tem_var_system_mdoule#derive "derive".
+  subroutine derivePressureIncomp_forElemFromState( input, nElems, tree, requVar,    &
+    &                      globSys, stencil, res, conf, bary, minIDbound, time )
+    !---------------------------------------------------------------------------
+    !> The results are stored in the same manner in the same array.
+    real(kind=rk), target, intent(inout)    :: input(:)
+    !> number of elements in the passed input array (io-buffer)
+    integer, intent(in) :: nElems
+    !> tree to work on (e.g. global or local from tracking)
+    type(treelmesh_type), intent(in) :: tree
+    !> required variable
+    type( tem_variable_type ), intent(in) :: requVar
+    !> global variable system
+    type( tem_var_system_type ), intent(in) :: globSys
+    !> stencil information
+    type( tem_stencilHeader_type ) :: stencil
+    !> optional output array (mainly used in tracking)
+    real(kind=rk), optional, target, intent(out)    :: res(:)
+    !> aotus type for solver specific string or lua file (for harvester)
+    type( flu_state ) :: conf
+    !> array of barycenters of the elements
+    real(kind=rk), optional, intent(in) :: bary(:,:)
+    !> starting bound for the barycenters and the treeID list in the tree in
+    !! case of a chunk usage (needed for spacetime functions and meshInfo)
+    integer, optional, intent(in) :: minIDbound
+    type(tem_time_type), optional, intent(in) :: time
+    !---------------------------------------------------------------------------
+    !local variable
+    integer :: iDep, dep, iElem, buf_start, buf_end
+    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
+    real(kind=rk) :: dens, rho0
+    !---------------------------------------------------------------------------
+
+    if( present(res) )then
+      local_res => res
+    else
+      local_res => input
+    end if
+    !get reference density
+    rho0 = getFieldVariable( conf = conf,                                &
+      &                      varLabel = trim(requVar%label),             &
+      &                      varName = 'pressure',                       &
+      &                      fieldVar = 'rho0',                          &
+      &                      fieldProp = 'fluid' )
+
+    do iElem = 1, nElems
+      dens = 0.0_rk
+      do iDep = 1, size( globSys%variable%val( requVar%mapping )%varDep )
+        ! get the position of the depending variable in the global varSys
+        dep = globSys%variable%val( requVar%mapping )%varDep( iDep )
+      
+        ! set the start and end position in the state array
+        buf_start = (iElem - 1)*globSys%nScalars                               &
+          &   + globSys%variable%val( dep )%varPos( 1 )
+      
+        buf_end   = (iElem - 1)*globSys%nScalars                               &
+          &   + globSys%variable%val( dep )%varPos(                            &
+          &                            globSys%variable%val( dep )%nComponents )
+      
+        ! derive the density and store it in local_res
+        ! derive the density and store it in the input (for more than
+        ! 1 species (iDep) this will be the mixture density)
+        dens = dens + sum( input( buf_start : buf_end ))
+      
+      end do  ! iDep
+      local_res( iElem ) = (dens - rho0)*cs2
+    end do  ! iElem
+
+  end subroutine derivePressureIncomp_forElemFromState
+!******************************************************************************!
+
+
+!******************************************************************************!
   !> Calculate the velocity magnitude of a given element number with the given
   !! input vector (sum up all values)
   !!
 
 !******************************************************************************!
   !> Calculate the kinetic energy as
-  !!  E = 1/2 rho sqrt( ux^2 + uy^2 +uz^2 )
+  !!  E = 1/2 rho ( ux^2 + uy^2 +uz^2 )
   !!
   !! The interface has to comply to the abstract interface
   !! tem_var_system_mdoule#derive "derive".
     real(kind=rk) :: vel(3) ! temporary array for velocity
     integer :: iDir, iElem, dep, nComp
     integer :: buf_start, buf_end
-    real(kind=rk) :: ucx, usq
     real(kind=rk) :: weights(stencil%QQ) ! weights read from spec_string
     !---------------------------------------------------------------------------
   
       end do
       vel = vel / dens
       
-      local_res( iElem ) = 0.5_rk*dens*sqrt(                                   &
+      local_res( iElem ) = 0.5_rk*dens*(                                       &
         &       vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3) )
     end do ! iElem
   

File source/mus_derived/mus_derivedQuantitiesPhysics_module.f90

View file
  • Ignore whitespace
     &                                     deriveVelMag_forElemFromState
   use mus_derivedQuantitiesIncomp_module, only: deriveVelocityIncomp_forElemFromState,    &
     &                                     deriveVelMagIncomp_forElemFromState,      &
+    &                                     derivePressureIncomp_forElemFromState,      &
     &                                     deriveShearIncomp_forElemFromState,       &
     &                                     deriveKineticEnergyIncomp_forElemFromState,       &
     &                                     deriveWSSIncomp_forElemFromState
   public :: deriveTemperaturePhy_forElemFromState
 
   !incompressible model
+  public :: derivePressureIncompPhy_forElemFromState
   public :: deriveVelocityIncompPhy_forElemFromState
   public :: deriveWSSIncompPhy_forElemFromState
   public :: deriveShearIncompPhy_forElemFromState
 
 
 !******************************************************************************!
+  !> Calculate the pressure of a given set of elements (sum up all links).
+  !!
+  !! Pressure calculation according to the isentropic equation of state for 
+  !! the LBM \f$ p = \rho c_s^2 \f$
+  !! with the calculation of density as in deriveDensityPhy_forElemFromState
+  !!
+  !! The interface has to comply to the abstract interface
+  !! tem_var_system_mdoule#derive "derive".
+  subroutine derivePressureIncompPhy_forElemFromState( input, nElems, tree, requVar, &
+    &                      globSys, stencil, res, conf, bary, minIDbound, time )
+    !---------------------------------------------------------------------------
+    !> The results are stored in the same manner in the same array.
+    real(kind=rk), target, intent(inout)    :: input(:)
+    !> number of elements in the passed input array (io-buffer)
+    integer, intent(in) :: nElems
+    !> tree to work on (e.g. global or local from tracking)
+    type(treelmesh_type), intent(in) :: tree
+    !> required variable
+    type( tem_variable_type ), intent(in) :: requVar
+    !> global variable system
+    type( tem_var_system_type ), intent(in) :: globSys
+    !> stencil information
+    type( tem_stencilHeader_type ) :: stencil
+    !> optional output array (mainly used in tracking)
+    real(kind=rk), optional, target, intent(out)    :: res(:)
+    !> aotus type for solver specific string or lua file (for harvester)
+    type( flu_state ) :: conf
+    !> array of barycenters of the elements
+    real(kind=rk), optional, intent(in) :: bary(:,:)
+    !> starting bound for the barycenters and the treeID list in the tree in
+    !! case of a chunk usage (needed for spacetime functions and meshInfo)
+    integer, optional, intent(in) :: minIDbound
+    type(tem_time_type), optional, intent(in) :: time
+    !---------------------------------------------------------------------------
+    !local variable
+    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
+    real(kind=rk) :: pressFac
+    !---------------------------------------------------------------------------
+
+    if( present(res) )then
+      local_res => res
+    else
+      local_res => input
+    end if
+
+    pressFac = getConversionFac( conf = conf, facName = 'press' )
+
+    call derivePressureIncomp_forElemFromState( input = input,                       &
+      &                                   nElems = nElems,                     &
+      &                                   tree = tree,                         &
+      &                                   requVar = requVar,                   &
+      &                                   globSys = globSys,                   &
+      &                                   stencil = stencil,                   &
+      &                                   res = local_res,                     &
+      &                                   conf = conf,                         &
+      &                                   bary = bary,                         &
+      &                                   minIDbound = minIDbound,             &
+      &                                   time = time )
+  
+   local_res = local_res*pressFac
+
+  end subroutine derivePressureIncompPhy_forElemFromState
+!******************************************************************************!
+
+
+!******************************************************************************!
   !> Calculate the velocity of a given element number with the given input
   !!        vector (sum up all values)
   !!

File source/mus_derived/mus_derivedQuantities_module.f90

View file
  • Ignore whitespace
   use tem_variable_module, only: tem_variable_type
   use tem_stencil_module, only: tem_stencilHeader_type
   use tem_tools_module, only: tem_write
+  use tem_topology_module, only: tem_levelOf
   use tem_time_module, only: tem_time_type
   use treelmesh_module, only: treelmesh_type
   use mus_solSpecHelpers_module, only: getWeights, getFieldVariable
     type(tem_time_type), optional, intent(in) :: time
     !---------------------------------------------------------------------------
     !local variable
-    integer :: iDep, dep, iElem, buf_start, buf_end
+    integer :: iElem
     real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    real(kind=rk) :: dens
     !---------------------------------------------------------------------------
 
     if( present(res) )then
         &       + globSys%variable%val( pdf_pos )%varPos( pdf_nComps )
       
       ! first compute shear stress, save it to array tau
+      ! here we pass elements one by one, so minIDbound need to adapt to it.
       tau = 0._rk
       call globSys%derive( shear_pos )%do(                                     &
         &                          input = input( buf_start:buf_end ),         &
         &                          res = tau,                                  &
         &                          conf = conf,                                &
         &                          bary = bary,                                &
-        &                          minIDbound = minIDbound,                    &
+        &                          minIDbound = (minIDbound + iElem - 1),      &
         &                          time = time )
       
       a2 = -1._rk * ( tau(1) + tau(4) + tau(6) )
   !! \f$ f^{neq}_i = f_i - f^{eq}_i\f$ is the non-equilibirium density.\n
   !! For more information, please refer to:\n
   !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann
-  !! simulations. Physical Review E. 2009;79(4):1-14.
-  !! \todo: update formula
+  !! simulations. Physical Review E. 2009;79(4):1-14.\n
   !!
+  !! For multi-level mesh, Omega on finer level needs to be adjusted in order to
+  !! get the correct shearstress calculation.\n
+  !! First, we defines c as the dx ratio between finer and coarse level.\n
+  !! \f$ c={ \Delta dx }_{ c }/{ \Delta dx }_{ f } \f$
+  !! Then the viscosity on the different levels must satisfy:\n
+  !! \f$ \frac { { \nu  }_{ f } }{ { \nu  }_{ c } } =c \f$
+  !! This constrain leads to a relationship of omega on different levels:\n
+  !! \f$ {\omega}_f = \frac {1}{ {\lambda}(\frac{1}{{\omega}_c}-0.5)+0.5 } \f$
+  !! For more information, please refer to:\n
+  !! Manuel H, Harald K, Joerg B, Sabine R. Aeroacoustic validation of the
+  !! lattice boltzmann method on non-uniform grids. ECCOMAS 2012
+  !! 
   subroutine deriveShear_forElemFromState( input, nElems, tree, requVar,       &
     &                     globSys,  stencil, res, conf, bary, minIDbound, time )
     !---------------------------------------------------------------------------
     real(kind=rk) :: stenCxDir(6, stencil%QQ) ! stencil cxDir
     real(kind=rk) :: tau(requVar%nComponents), NonEq(stencil%QQ)
     integer       :: iDir, iElem, buf_start, buf_end
-    real(kind=rk) :: omega, conv, cxsq( stencil%QQ )
+    integer       :: my_level, dx_ratio
+    real(kind=rk) :: ref_omega, my_omega, conv, cxsq( stencil%QQ )
     !---------------------------------------------------------------------------
 
     if( present(res) )then
       local_res => input
     end if
 
-    ! Calculate Omega
-    omega = getFieldVariable( conf = conf,                                     &
-      &                       varLabel = trim(requVar%label) ,                 &
-      &                       varName = 'shearstress',                         &
-      &                       fieldVar = 'omega',                              &
+    ! Get the omega on the coarst level as reference omega
+    ref_omega = getFieldVariable( conf = conf,                     &
+      &                       varLabel = trim(requVar%label) ,     &
+      &                       varName = 'shearstress',             &
+      &                       fieldVar = 'omega',                  &
       &                       fieldProp = 'fluid' )
 
-    if( omega /= 1._rk ) then
-      !JQ: conversion factor. Use it to convert post-coll to pre-coll
-      conv = 1._rk / ( 1._rk - omega )
-    else
-      conv = 0._rk
-      call tem_write( "Caution! Shear stress can NOT be calculated "//         &
-        &             "when omega equals to 1 !" )
-      call tem_write( "Setting it to zero!" )
-    endif
-
     ! convert stencil cxDir into real type
     stenCxDir(1, :) =   real(stencil%cxDir(1, :), kind=rk)  &
       &               * real(stencil%cxDir(1, :), kind=rk)
 
 
     do iElem = 1, nElems
+      ! get my current level
+      my_level = tem_levelOf( tree%treeID( minIDbound + iElem ) )
+      ! dx ratio between my level and minimum level
+      dx_ratio = 2 ** ( my_level - tree%global%minLevel )
+      ! get the correct omega value
+      my_omega = 1._rk / &
+        &   ( real(dx_ratio, kind=rk) * ( 1._rk / ref_omega - 0.5_rk ) + 0.5_rk )
+
+      if( my_omega /= 1._rk ) then
+        !JQ: conversion factor. Use it to convert post-coll to pre-coll
+        conv = 1._rk / ( 1._rk - my_omega )
+      else
+        conv = 0._rk
+        call tem_write( "Caution! Shear stress can NOT be calculated "//         &
+          &             "when omega equals to 1 !" )
+        call tem_write( "Setting it to zero!" )
+      endif
+
       tau = 0._rk
       nonEq = 0._rk
-      
+
       ! set the begin and end position in the input state array
       buf_start = (iElem-1) * globSys%nScalars + 1
       buf_end =   (iElem-1) * globSys%nScalars + stencil%QQ
-      
+
       ! Calculate NonEquilibrium density
       ! for PULL, this is post-collision Non-Eq
       ! for PUSH, this is  pre-collision Non-Eq
         &               bary = bary,                                           &
         &               minIDbound = minIDbound,                               &
         &               time = time )
-      
+
       ! Shear stress requires pre-collision NonEq
       ! Convert post-collision to pre-collision by a factor of 1/(1-omega)
       nonEq = nonEq * conv
-      
+
       do iDir = 1, stencil%QQ
         tau(1) = tau(1) + nonEq(iDir) * ( stenCxDir(1, iDir) - cxsq(iDir) * div1_3 ) ! tau_x
         tau(4) = tau(4) + nonEq(iDir) * ( stenCxDir(4, iDir) - cxsq(iDir) * div1_3 ) ! tau_y
         tau(3) = tau(3) + nonEq(iDir) * stenCxDir(3, iDir) ! tau_xz
         tau(5) = tau(5) + nonEq(iDir) * stenCxDir(5, iDir) ! tau_yz
       enddo ! iDir
-      
-      local_res( (iElem-1)*6+1 : iElem*6 ) = ( omega * 0.5_rk - 1._rk ) * tau(:)
+
+      local_res( (iElem-1)*6+1 : iElem*6 ) = ( my_omega * 0.5_rk - 1._rk ) * tau(:)
     enddo ! iElem
 
   end subroutine deriveShear_forElemFromState
 
 !******************************************************************************!
   !> Calculate the kinetic energy as
-  !!  E = 1/2 rho sqrt( ux^2 + uy^2 +uz^2 )
+  !!  E = 1/2 rho ( ux^2 + uy^2 +uz^2 )
   !!
   !! The interface has to comply to the abstract interface
   !! tem_var_system_mdoule#derive "derive".
     real(kind=rk) :: vel(3) ! temporary array for velocity
     integer :: iDir, iElem, dep, nComp
     integer :: buf_start, buf_end
-    real(kind=rk) :: ucx, usq
     real(kind=rk) :: weights(stencil%QQ) ! weights read from spec_string
     !---------------------------------------------------------------------------
   
       end do
       vel = vel / dens
       
-      local_res( iElem ) = 0.5_rk*dens*sqrt(                                   &
-        &       vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3) )
+      local_res( iElem ) = 0.5_rk*dens*(                                       &
+        &       vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3))
     end do ! iElem
   
   end subroutine deriveKineticEnergy_forElemFromState

File source/mus_derived/mus_physics_type.f90

View file
  • Ignore whitespace
     real(kind=rk) :: force 
     !> Pressure (kg/m/s^2) = rho0*dx^2/dt^2
     real(kind=rk) :: press
-    !> Energy (kg*m^2/s^2) = rho0*dx^5/dt^2
+    !> specific Energy (rho*m^2/s^2) = rho0*dx^2/dt^2
     real(kind=rk) :: energy
   end type mus_convertFac_type
 
     !Pressure
     me%fac%press = me%rho0*me%dx*me%dx/(me%dt*me%dt)
     !Energy
-    me%fac%energy = me%rho0*me%dx*me%dx*me%dx*me%dx*me%dx/(me%dt*me%dt)
+    me%fac%energy = me%rho0*me%dx*me%dx/(me%dt*me%dt)
 
   end subroutine computeConvertFac
 !******************************************************************************!
     character(len=*) :: fun_str
     !---------------------------------------------------------------------------
     character(len=256) :: dtVel, dtOmega, dtVisc, dxLevel, dxnL, omega_dt, &
-      & vel_dt, tempT0
+      & vel_dt
     !---------------------------------------------------------------------------
     write(dxLevel,*) 'function getdxFromLevel(t) &
       & return t.len_bnd/math.pow(2,t.level) end'

File source/mus_derived/mus_variable_module.f90

View file
  • Ignore whitespace
     &                                     deriveTemperature_forElemFromState,            &
     &                                     deriveVelMag_forElemFromState
   use mus_derivedQuantitiesIncomp_module, only: deriveVelocityIncomp_forElemFromState,    &
-    &                                     deriveVelMagIncomp_forElemFromState,      &
-    &                                     deriveEquilibriumIncomp_forElemFromState, &
-    &                                     deriveEquilibriumIncomp_forElemFromMacro, &
-    &                                     deriveShearIncomp_forElemFromState,       &
+    &                                     derivePressureIncomp_forElemFromState,          &
+    &                                     deriveVelMagIncomp_forElemFromState,            &
+    &                                     deriveEquilibriumIncomp_forElemFromState,       &
+    &                                     deriveEquilibriumIncomp_forElemFromMacro,       &
+    &                                     deriveShearIncomp_forElemFromState,             &
     &                                     deriveKineticEnergyIncomp_forElemFromState,     &
     &                                     deriveWSSIncomp_forElemFromState
   use mus_derivedQuantitiesMultispc_module, only: deriveMolFracMultispc_forElemFromState,&
     &                                     deriveEquilibriumMultispc_2d_forElemFromState, &
     &                                     deriveEquilibriumMultispc_2d_forElemFromMacro
   use mus_derivedQuantitiesPhysics_module, only: derivePressurePhy_forElemFromState, &
+    &                                     derivePressureIncompPhy_forElemFromState,        &  
     &                                     deriveVelocityPhy_forElemFromState,        &  
     &                                     deriveWSSPhy_forElemFromState,             &  
     &                                     deriveKineticEnergyPhy_forElemFromState,       &
     pressure_var%variable%nComponents = 1
     allocate(pressure_var%variable%depVars(1))
     pressure_var%variable%depVars(1) = 'pdf'
-    pressure_var%do => derivePressure_forElemfromState
 
     pressurePhy_var%variable%mapping = 1
     pressurePhy_var%variable%label = 'pressure_phy'
     pressurePhy_var%variable%nComponents = 1
     allocate(pressurePhy_var%variable%depVars(1))
     pressurePhy_var%variable%depVars(1) = 'pdf'
-    pressurePhy_var%do => derivePressurePhy_forElemfromState
 
     ! velocity
     velocity_var%variable%label = 'velocity'
         wss_var%do            => deriveWSS_forElemfromState
         kineticEnergy_var%do      => deriveKineticEnergy_forElemfromState
         temperature_var%do      => deriveTemperature_forElemfromState
+        pressure_var%do => derivePressure_forElemfromState
 
         velocityPhy_var%do       => deriveVelocityPhy_forElemfromState
         velMagPhy_var%do         => deriveVelMagPhy_forElemfromState
         wssPhy_var%do            => deriveWSSPhy_forElemfromState
         kineticEnergyPhy_var%do  => deriveKineticEnergyPhy_forElemfromState
         temperaturePhy_var%do    => deriveTemperaturePhy_forElemfromState
+        pressurePhy_var%do => derivePressurePhy_forElemfromState
 
         ! create derVarList
         allocate( derVarList( 17 ))
         shearStress_var%do    => deriveShearIncomp_forElemfromState
         wss_var%do            => deriveWSSIncomp_forElemfromState
         kineticEnergy_var%do      => deriveKineticEnergyIncomp_forElemfromState
+        pressure_var%do       => derivePressureIncomp_forElemfromState
 
         velocityPhy_var%do       => deriveVelocityIncompPhy_forElemfromState
         velMagPhy_var%do         => deriveVelMagIncompPhy_forElemfromState
         shearStressPhy_var%do    => deriveShearIncompPhy_forElemfromState
         wssPhy_var%do            => deriveWSSIncompPhy_forElemfromState
-        kineticEnergyPhy_var%do      => deriveKineticEnergyIncompPhy_forElemfromState
+        kineticEnergyPhy_var%do  => deriveKineticEnergyIncompPhy_forElemfromState
+        pressurePhy_var%do       => derivePressureIncompPhy_forElemfromState
         ! create derVarList
         allocate( derVarList( 15 ))
         derVarList = (/ density_var, velocity_var, pressure_var,               &
           &                  nComponents = derVarList(iVar)%variable%nComponents,  &
           &                  pos = pos,                                            &
           &                  depVar = ['all'] )
-        var_sys%derive( pos )%do => derVarList( iVar )%do
+        ! check if the variable has been added (only appended if requirements are
+        ! fullfilled)
+        if(pos .eq. 0 )then
+          write( buffer_depVars,*) trim(derVarList(iVar)%variable%depVars(1))
+          call tem_write('------------------------------------------------')
+          call tem_write('The following variable has not been appended to the list')
+          call tem_write('label: '//trim(buffer)//trim(derVarList(iVar)%variable%label))
+          call tem_write('dep: '//trim(buffer)//trim(buffer_depVars))
+          call tem_write('------------------------------------------------')
+        ! check if the variable to append is already a state variable
+        else if( pos .gt. var_sys%nVars )then
+          var_sys%derive( pos )%do => derVarList( iVar )%do
+          var_sys%variable%val( pos )%mapping = pos
+        end if
       endif
     end do ! loop derVarList
 
     do iVar = 1, size(tVar)
       if ( trim(tVar(iVar)%st_fun(1)%fun_kind) .ne. 'none') then
         call tem_var_append( sys = var_sys, var = tVar( iVar ), pos = pos)
-        var_sys%derive(pos)%do => evalSpacetime
-        var_sys%variable%val( pos )%mapping = pos
+        ! check if the variable has been added (only appended if requirements are
+        ! fullfilled)
+        if(pos .eq. 0 )then
+          write( buffer_depVars,*) trim(tVar(iVar)%depVars(1))
+          call tem_write('------------------------------------------------')
+          call tem_write('The following variable has not been appended to the list')
+          call tem_write('label: '//trim(buffer)//trim(tVar(iVar)%label))
+          call tem_write('dep: '//trim(buffer)//trim(buffer_depVars))
+          call tem_write('------------------------------------------------')
+        ! check if the variable to append is already a state variable
+        else if( pos .gt. var_sys%nVars )then
+          var_sys%derive(pos)%do => evalSpacetime
+          var_sys%variable%val( pos )%mapping = pos
+        end if
       endif
     enddo
 
     do iVar = 1, size(meshInfoList)
       call tem_var_append( sys = var_sys, var = meshInfoList( iVar )%variable, &
         &                  pos = pos )
-      var_sys%derive(pos)%do => meshInfoList( iVar )%do
-      var_sys%variable%val( pos )%mapping = pos
+      ! check if the variable has been added (only appended if requirements are
+      ! fullfilled)
+      if(pos .eq. 0 )then
+        write( buffer_depVars,*) trim(meshInfoList(iVar)%variable%depVars(1))
+        call tem_write('------------------------------------------------')
+        call tem_write('The following variable has not been appended to the list')
+        call tem_write('label: '//trim(buffer)//trim(meshInfoList(iVar)%variable%label))
+        call tem_write('dep: '//trim(buffer)//trim(buffer_depVars))
+        call tem_write('------------------------------------------------')
+      ! check if the variable to append is already a state variable
+      else if( pos .gt. var_sys%nVars )then
+        var_sys%derive( pos )%do => meshInfoList( iVar )%do
+        var_sys%variable%val( pos )%mapping = pos
+      end if
     enddo
 
 

File source/tem_balance_module.f90

View file
  • Ignore whitespace
 
   !> Load balancing information and control 
   type tem_balance_type
+    !> load balancing type
+    character(len=labelLen) :: balanceKind
     !> is load balancing activated?
     logical :: active
     !> Imbalance: max( timeMainLoop( iProc )) - average( timeMainLoop )
     real( kind=rk ) :: imbalance
+    !> base weight for each element which is set to be
+    !! weight( element ) = baseWeight^( level(element) -refLevel )
+    real( kind=rk ) :: baseWeight
   end type tem_balance_type
 
   contains
       &                    val = me%active, &
       &                    ErrCode = iError, default = .false. )
 
+    call aot_get_val(L = conf, thandle = thandle, &
+      &                    key = 'kind',          &
+      &                    val = me%balanceKind,  &
+      &                    ErrCode = iError, default = 'simple' )
+
+    call aot_get_val(L = conf, thandle = thandle, &
+      &                    key = 'baseweight',          &
+      &                    val = me%baseWeight,  &
+      &                    ErrCode = iError, default = 2._rk )
+
     if( present( verbose )) then
       if( verbose ) then
-        if( me%active ) &
-        call tem_write( 'Load Balancing is active' )
+        if( me%active ) then
+          call tem_write( 'Load Balancing is active' )
+          call tem_write( '  kind:       '//trim(me%balanceKind))
+          call tem_write( '  baseweight: ', me%baseWeight )
+        endif
       endif
     endif
     call aot_table_close(L=conf, thandle=thandle)

File source/tem_canonical_module.f90

View file
  • Ignore whitespace
     !> kind of canonicalND (line, plane, point, box)
     character( len=labellen ) :: kind
     logical :: active(3)     !< identify which vectors are active (not equal 0)
+    !> dimension of canonical object 
+    !! nDim=0 - point
+    !! nDim=1 - line
+    !! nDim=2 - plane
+    !! nDim=3 - box
+    integer :: nDim
     !> To choose what to do with intersection of box
     !! if only_surface = true than the only the surface of the object
     !! is intersected
     integer :: vError(3), errfatal(3)
     ! counter and total number of vectors
     integer :: iVec, nVecs   
-    integer :: nActive
     ! lua handles
     integer :: vec_thandle, halfvec_thandle, single_thandle, halfsingle_thandle
     real(kind=rk), allocatable :: local_vecs(:,:)
       &               key = 'distribution', default = 'equal' )
 
     me%active = .false.
-    nActive = 0
+    me%nDim = 0
 
     ! check the kind and the active vectors
     do iVec = 1, 3
       vecAbsSum = me%vec(1,iVec) + me%vec(2,iVec) + me%vec(3,iVec)
       if(vecAbsSum .ne. 0._rk) then
         me%active(iVec) = .true.
-        nActive = nActive + 1
+        me%nDim = me%nDim + 1
       end if
     end do
 
-    if(nActive .eq. 0) me%kind = 'point'
-    if(nActive .eq. 1) me%kind = 'line'
-    if(nActive .eq. 2) me%kind = 'plane'
-    if(nActive .eq. 3) me%kind = 'box'
+    if(me%nDim .eq. 0) me%kind = 'point'
+    if(me%nDim .eq. 1) me%kind = 'line'
+    if(me%nDim .eq. 2) me%kind = 'plane'
+    if(me%nDim .eq. 3) me%kind = 'box'
 
     !if canonical object is box then check for only_surface
     if(me%kind=='box') then
 
     call tem_write( 'Initializing canonical shapes')
 
+    !Initialize countElems
+    countElems = 0
+
     ! loop over all canoND objects
     do iCano = 1, size(me)
       ! \todo: SZ: set an upper limit in the load routine of the shapes,

File source/tem_comm_module.fpp

View file
  • Ignore whitespace
     !! \todo request handle array could exist during complete code runtime
     integer :: rq_handle( recv%nProcs + send%nProcs )
     integer :: status( mpi_status_size, recv%nProcs + send%nProcs )
-    integer :: iErr             !< error flag
+    integer :: iErr             ! error flag
     integer :: iProc
 
     !> Values for send me must have been copied from the actual state array
     do iProc = 1, recv%nProcs
       !> Start receive communications
-      call mpi_irecv(                                       &
-       &      state, & !< me
-       &      1, & ! number of datatype elements
-       &      recv%buf_?tname?(iProc)%memIndexed, & !< data type
-       &      recv%proc(iProc), & !< sending process
-       &      message_flag, & 
-       &      process%comm, & !< communicator
-       &      rq_handle(iProc), & !< handle
-       &      iErr ) !< error status
+      call mpi_irecv(                                                          &
+                      ! me
+       &              state,                                                   &
+                      ! number of datatype elements
+       &              1,                                                       &
+                      ! data type
+       &              recv%buf_?tname?(iProc)%memIndexed,                      &
+                      ! sending process
+       &              recv%proc(iProc),                                        &
+                      ! tag
+       &              message_flag,                                            & 
+                      ! communicator
+       &              process%comm,                                            &
+                      ! handle
+       &              rq_handle(iProc),                                        &
+                      ! error status
+       &              iErr )
 
     enddo
 
     !>  Start the sending communications
     do iProc = 1, send%nProcs
-      call mpi_isend(                                       &
-       &      state, & !< me
-       &      1, & !< number of datatype elements
-       &      send%buf_?tname?(iProc)%memIndexed, & !< data type
-       &      send%proc(iProc), & !< target process
-       &      message_flag, & !< tag
-       &      process%comm, & !< communicator
-       &      rq_handle( iProc + recv%nProcs ), & !< handle
-       &      iErr ) !< error status
+      call mpi_isend(                                                          &
+                      ! me
+       &              state,                                                   &
+                      ! number of datatype elements
+       &              1,                                                       & 
+                      ! data type
+       &              send%buf_?tname?(iProc)%memIndexed,                      &
+                      ! target process
+       &              send%proc(iProc),                                        &
+                      ! tag
+       &              message_flag,                                            &
+                      ! communicator
+       &              process%comm,                                            &
+                      ! handle
+       &              rq_handle( iProc + recv%nProcs ),                        &
+                      ! error status
+       &              iErr )
 
     enddo !< iProc
 
     !> Wait for above communications to complete
-    call mpi_waitall( &
-      & send%nProcs+recv%nProcs, & !< number of reqs to wait for
-      & rq_handle, & !< request handles
-      & status, & !< mpi status
-      & iErr ) !< error status
+    call mpi_waitall(                                                          &
+                      ! number of reqs to wait for
+      &               send%nProcs+recv%nProcs,                                 &
+                      ! request handles
+      &               rq_handle,                                               & 
+                      ! mpi status
+      &               status,                                                  & 
+                      ! error status
+      &               iErr )
 
   end subroutine comm_typed_isend_irecv_?tname?
   !****************************************************************************!

File source/tem_construction_module.f90

View file
  • Ignore whitespace
 module tem_construction_module
   use mpi
   use env_module, only: long_k, rk, globalMaxLevels, dbgUnit, io_buffer_size, &
-                        chatter_none, chatter_little, chatter_aLot, chatter_all
+                        chatter_none, chatter_little, chatter_aLot, chatter_all, pathLen
   use tem_aux_module, only: tem_abort
   use tem_dyn_array_module, only: dyn_longArray_type, dyn_intArray_type, &
     &                             dyn_realArray_type, &
   public :: tem_cleanupDependencyArrays
   public :: tem_init_elemLevels
   public :: tem_eTypeOfId
+  public :: tem_eTypeOfTotalPos
   public :: tem_treeIDinTotal
   public :: tem_killArrays
   public :: tem_find_depProc
     !! gets the source elements for the interpolation
     integer                  ::  dependencyLevel = -1
     !> position of the source elements in the totalList
-    type( dyn_intArray_type ) :: elem
+    type( grw_intArray_type ) :: elem
     !> position of the source elements in the subbuffer
-    type( dyn_intArray_type ) :: elemBuffer
+    type( grw_intArray_type ) :: elemBuffer
     !> Interpolation weight for each source element specified above
     real(kind=rk), allocatable :: weight(:)
     real(kind=rk) :: coord(3)
 
     !---------new types for element description---------------------
     type(tem_element_type) :: elem
+    !> This list includes treeIDs for which additionally
+    !! neighbors have to be identified
     type(dyn_longArray_type) :: require
     type(dyn_haloArray_type) :: halosPrc
     !> points to the position of the element in the elem array
     call check_additionalComm( levelDesc, process, doAdditional, tree%global%minlevel ) 
    
     if( doAdditional ) then
+      ! passing only first stencil as this is the required compute stencil
       call identify_additionalNeigh( tree, process, levelDesc,               &
         &                        pathFirst, pathLast, computeStencil(1), verbosity)
       call communicate_elements( tree, process, levelDesc, commPattern,          &
           &                     nUnit = dbgUnit,  stencil = .true.,         &
           &                     verbose = verbose,                          &
           &                     string = 'after initializing level elements')
+        if( me( iLevel )%require%nVals .gt. 0 ) then
+        call tem_require_dump( me = me( iLevel )%require, compact = .true., &
+          &                    nUnit = dbgUnit,                             &
+          &                    verbose = verbose,                           &
+          &                    string = 'after initializing level elements')
+        end if
       end do
     end if
 
   end subroutine tem_killArrays
   !*****************************************************************************
 
-!  !*****************************************************************************
-!  !>  Invoke the identify_elements for each neighbor of the stencil 
-!  !! and store the positions of the encountered elements
-!  subroutine identify_stencilNeigh( iElem, elemStencil, tree, pathFirst, pathLast, &
-!    &                                     levelDesc, elemPos, prc, &
-!    &                                     computeStencil, nesting, verbose )
-!    !--------------------------------------------------------------------------
-!    integer, intent(in) :: iElem     !< element position in levelDesc to identify
-!    type(treelmesh_type), intent(in) :: tree !< the tree
-!    type(tem_path_type), intent(in) :: pathFirst(:), pathLast(:)
-!    type(tem_levelDesc_type), intent(inout) :: levelDesc(tree%global &
-!      &                                                  %minLevel:)
-!    type(tem_comm_env_type), intent(in) :: prc    !< process
-!    integer, intent(out) :: elemPos
-!    type(tem_stencilHeader_type), intent(in) :: computeStencil
-!    integer, intent(in) :: nesting !< nesting level
-!    !> position in levelDesc( iLevel )%Elem where the element was reconstructed
-!    integer, intent(in) :: verbose
-!    !--------------------------------------------------------------------------
-!    !--------------------------------------------------------------------------
-! 
-!  end subroutine identify_stencilNeigh
-!  !*****************************************************************************
+  !*****************************************************************************
+  !>  Invoke the identify_elements for each neighbor of the stencil 
+  !! and store the positions of the encountered elements
+  subroutine identify_stencilNeigh( iElem, iLevel, iStencil, tree, pathFirst, pathLast, &
+    &                               levelDesc, prc, &
+    &                               computeStencil, nesting, verbose )
+    !--------------------------------------------------------------------------
+    integer, intent(in) :: iElem     !< element position in levelDesc to identify
+    integer, intent(in) :: iLevel    !< element level
+    integer, intent(in) :: iStencil  !< stencil within the element to act on 
+    type(treelmesh_type), intent(in) :: tree !< the tree
+    type(tem_path_type), intent(in) :: pathFirst(:), pathLast(:)
+    type(tem_levelDesc_type), intent(inout) :: levelDesc(tree%global &
+      &                                                  %minLevel:)
+    type(tem_comm_env_type), intent(in) :: prc    !< process
+    type(tem_stencilHeader_type), intent(in) :: computeStencil
+    integer, intent(in) :: nesting !< nesting level
+    !> position in levelDesc( iLevel )%Elem where the element was reconstructed
+    integer, intent(in) :: verbose
+    !--------------------------------------------------------------------------
+    integer :: iStencilElem, elemPos
+    integer :: neighIDpos
+    integer(kind=long_k) :: nnTreeID
+    !--------------------------------------------------------------------------
+    ! identify all the compute neighbors of the current element
+    do iStencilElem = 1, computeStencil%QQN
+      neighIDpos = levelDesc( iLevel )%elem%stencil%val( iElem )           &
+        &                  %val( iStencil )%tIDpos( iStencilElem )
+      if( neighIDpos .gt. 0 ) then
+        nnTreeID = levelDesc( iLevel )%elem%neighID%val( iElem )             &
+          &                  %val( neighIDpos )
+        ! This call might add new halo elements
+        if ( nnTreeID .gt. 0_long_k ) then
+          call identify_elements( nTreeID = nnTreeID, tree = tree,               &
+            &                     pathFirst = pathFirst,                         &
+            &                     pathLast = pathLast,                           &
+            &                     levelDesc = levelDesc,                         &
+            &                     elemPos   = elemPos,                           &
+            &                     prc = prc,                                     &
+            &                     computeStencil = computeStencil,               &
+            &                     nesting = nesting,                             &
+            &                     verbose = verbose )
+        else
+          elemPos = 0
+        end if
+      else
+        elemPos = 0
+      end if
+      ! And add the encountered neighbor elements to the current element's
+      ! stencil neighbors
+      levelDesc( iLevel )%elem%stencil%val( iElem )       &
+        &           %val( iStencil )%totalPos( iStencilElem ) = elemPos
+    end do
+ 
+  end subroutine identify_stencilNeigh
+  !*****************************************************************************
 
   !*****************************************************************************
   !> Check, on which partition a given element is located
     integer, intent(in) :: verbose
     !--------------------------------------------------------------------------
     integer(kind=long_k) :: children(8) ! child elements
-    integer :: nDepProcs, depProc(8) !< processes that this neigh depend on
+    integer :: nDepProcs, depProc(64) !< processes that this neigh depend on
     integer :: iChild, neighLevel
     type(tem_path_type) :: neighPath
     type(tem_stencilElement_type) :: emptyStencil(1)
     integer :: p_lb, p_ub !< process lower and upper bound
     integer :: nNesting !< new nesting level
     !> Position of the current element in the hash
-    integer :: hashpos
+    integer :: hashpos, neighIDpos, neighPos, elemNesting
     logical :: cacheHit, hashPos_avail
+    logical :: updated !< was the element updated during identify_local_element
     !--------------------------------------------------------------------------
     ! Increase the nesting
     nNesting = nesting + 1
     cacheHit = .false.
 
     if( verbose .ge. chatter_all ) then
-       write(dbgUnit,*) 'identifying required TreeID for', nTreeID
+      write(dbgUnit,*) 'identifying required TreeID for', nTreeID, 'nesting', nesting
     end if
     if( nTreeID .gt. 0 ) then
       neighLevel = tem_LevelOf( nTreeID )
           &           hash_elemPos( hashPos )) ) then
             cacheHit = .true.
           end if
+          if( levelDesc( neighlevel )%elem%needsUpdate%val(                    &
+          &           hash_elemPos( hashPos )) ) then
+            cacheHit = .false.
+           ! Set the needs update to false as it was now updated
+           levelDesc( neighLevel )%elem%needsUpdate%val( hash_elemPos( hashPos )) = .false.
+          end if
         end if
       end if
+      !\todo: Reactivate the cache! The problem is that the nesting of an
+      !element is updated to 0 without actually looping over its neighbors!
+      !cacheHit = .false.
       cachemiss: if( .not. cacheHit ) then
         ! Probably did not hit this element yet, put it in the hash now,
         ! and identify it.
           call tem_find_depProc(depProc, nDepProcs, &
             &                    neighPath, p_lb, p_ub, PathFirst, PathLast)
         end if ! local or non-local
-
+!        write(dbgUnit,*) 'nDepProcs', nDepProcs
         ! How many processes possess a part of the requested treeID
         if ( nDepProcs .eq. 1 ) then
+          updated = .false.
           ! might be a local, halo of same level or a ghostFromCoarse
           call single_process_neighbor( nTreeID = nTreeID,                     &
             &                           levelDesc = levelDesc,                 &
             &                           minLevel = tree%global%minlevel,       &
             &                           foundLevel = neighLevel,               &
             &                           elemPos = elemPos,                     &
+            &                           updated = updated,                     &
             &                           computeStencil = computeStencil,       &
             &                           nesting = nesting )
-          levelDesc( neighLevel )%elem%haloNesting%val( elemPos ) =            &
-             min( nesting,                                                     &
+          elemNesting = min( nesting,                                          &
                   levelDesc( neighLevel )%elem%haloNesting%val( elemPos ) )
-
-          if ( levelDesc( neighLevel )%elem%eType%val( elemPos )               &
-            &   == eT_ghostFromCoarserElem                                     &
-            & .and. nesting .lt. 1 ) then
-            if( verbose .ge. chatter_all ) then
-              write(dbgUnit,*) 'creating all parent neighbors for', &
-              & levelDesc( neighLevel )%elem%tID%val(elemPos)
-            end if
-            ! Create all elements required up to the actual existing fluid
-            ! element these include the neighbors of the parents. In a level
-            ! jump >1, these intermediate levels have to provide valid
-            ! quantities over two of their computation updates to account for
-            ! the recursive algorithm.
-            call create_allParentNeighbors(                           &
-              &    targetElem = elemPos,                              &
-              &    nesting = levelDesc( neighLevel )%elem%haloNesting &
-              &                                     %val(elemPos),    &
-              &    level = neighLevel, tree = tree,                   &
-              &    computeStencil = computeStencil,                   &
-              &    levelDesc = levelDesc,                             &
-              &    pathFirst = pathFirst, pathLast = pathLast,        &
-              &    prc = prc, verbose = verbose )
-          endif
+          if( nesting .lt. levelDesc( neighLevel )%elem%haloNesting%val( elemPos )) then
+            ! element needs updating as the current nesting is smaller than the
+            ! one in the element property
+            levelDesc( neighLevel )%elem%needsUpdate%val( elemPos ) = .true.
+            levelDesc( neighLevel )%elem%haloNesting%val( elemPos ) = elemNesting
+            updated = .true.
+          end if
+
+!          write(dbgUnit,*) 'done single proc neigh', nTreeID, ' nst', elemNesting, elemPos, 'upd?', updated
+          !updated = .true.
+          if( updated ) then
+            if ( levelDesc( neighLevel )%elem%eType%val( elemPos )               &
+              &   == eT_ghostFromCoarserElem ) then  
+!                call tem_element_dump( me = levelDesc( neighLevel )%elem, elemPos = elemPos, &
+!                   &         nUnit = dbgUnit, stencil = .true. )
+              if( nesting .lt. 1 ) then
+                ! Create the direct neighbors of the ghostFromCoarser
+!                write(dbgUnit,*) 'creating neighbor for gfC ', &
+!                & levelDesc( neighLevel )%elem%tID%val(elemPos)
+
+                ! identify all the compute neighbors of the current element
+                call identify_stencilNeigh( iElem = elemPos, iLevel = neighLevel,   &
+                       &             tree = tree, iStencil = 1,                     &
+                       &             pathFirst = pathFirst,                         &
+                       &             pathLast = pathLast,                           &
+                       &             levelDesc = levelDesc,                         &
+                       &             prc = prc,                                     &
+                       &             computeStencil = computeStencil,               &
+                       &             nesting = elemNesting+1,                     &
+                       &             verbose = verbose )
+
+              end if ! nesting < 1?
+        
+              if( verbose .ge. chatter_all ) then
+                write(dbgUnit,*) 'creating all parent neighbors for', &
+                & levelDesc( neighLevel )%elem%tID%val(elemPos)
+              end if
+              ! Create all elements required up to the actual existing fluid
+              ! element these include the neighbors of the parents. In a level
+              ! jump >1, these intermediate levels have to provide valid
+              ! quantities over two of their computation updates to account for
+              ! the recursive algorithm.
+              call create_allParentNeighbors(                           &
+                &    targetElem = elemPos,                              &
+                &    nesting = nNesting,                                &
+                &    level = neighLevel, tree = tree,                   &
+                &    computeStencil = computeStencil,                   &
+                &    levelDesc = levelDesc,                             &
+                &    pathFirst = pathFirst, pathLast = pathLast,        &
+                &    prc = prc, verbose = verbose )
+        
+            endif ! ghostFromCoarser?
+          endif ! updated?
 
         elseif( nDepProcs .gt. 1 ) then
           ! more than one depending processes
           end do
         else
           ! cell not existing. stop.
-          !write(*,*) '---------- WARNING!!! -------------------------'
+          write(*,*) '---------- WARNING!!! -------------------------'
           write(*,*) 'cell', nTreeID ,' not existing. ' // &
             &        'Number of dependencies', nDepProcs
-          write(*,*) '        coord',tem_coordOfId( nTreeID )
-          write(*,*) '        bary ',tem_baryOfId(  tree, nTreeID )
+          call tem_tIDinfo( me = nTreeID, tree = tree )
           write(*,*) 'WARNING!!! This should never occur and points to a bug' &
                      // ' or a buggy mesh.'
         end if
 
 
   !*****************************************************************************
-  !>
+  !> create all the neighbors of an element's parent
+  !!
+  !! Create all elements required up to the actual existing fluid
+  !! element these include the neighbors of the parents. In a level
+  !! jump >1, these intermediate levels have to provide valid
+  !! quantities over two of their computation updates to account for
+  !! the recursive algorithm.                                         
   recursive subroutine create_allParentNeighbors( &
     &                    targetElem, level, computeStencil, tree, levelDesc, &
     &                    pathFirst, pathLast, prc, verbose, nesting )
     !> position in levelDesc( iLevel )%Elem where the element was reconstructed
     integer, intent(in) :: verbose
     !--------------------------------------------------------------------------
-    integer :: iStencilElem, nNesting
-    integer(kind=long_k) :: nnTreeID
-    integer :: coarserLevel, neighPos, elemPos, neighIDpos
+    integer :: nNesting
+    integer(kind=long_k) :: nnTreeID  !< neighbor tree ID
+    integer(kind=long_k) :: cTreeID  !< current tree ID
+    integer :: coarserLevel, cPos, elemPos, neighIDpos, parentNesting
     !--------------------------------------------------------------------------
-    nNesting = nesting + 1
+    !nNesting = nesting + 1
     coarserLevel = level
-    neighPos = targetElem
-    do
-!      call tem_element_dump( me = levelDesc( coarserLevel )%elem &
-!        &                 %val( neighPos ), nunit = dbgUnit, stencil =.true. )
-      ! find the neighbors of all parent ghosts up to the existing fluid element
-      do iStencilElem = 1, computeStencil%QQN
-        neighIDpos = levelDesc( coarserLevel )%elem%stencil%val( neighPos )            &
-          &                  %val(1)%tIDpos( iStencilElem )
-!        write(dbgUnit, *) iStencilElem, coarserLevel,  'treeID from stencil', neighIDpos
-        nnTreeID = levelDesc( coarserLevel )%elem%neighID%val( neighPos )              &
-          &                  %val( neighIDpos )
-        ! This call might add new halo elements
-        if ( nnTreeID .gt. 0_long_k ) then
-        call identify_elements( nTreeID = nnTreeID, tree = tree,               &
+    ! The parent again must be computed, so the nesting has to be 0
+    parentNesting = 0
+    cPos = targetElem
+    ! exit if we have reached the minimal level
+    if( coarserLevel == tree%global%minlevel ) return
+
+    ! current element treeID
+    cTreeID = levelDesc( coarserLevel )%elem%tID%val( cPos )
+
+    ! Get the parent of the current treeID
+    coarserLevel = coarserLevel - 1
+    cTreeID = tem_parentOf( cTreeID )
+
+    ! ... and identify the parent (new current tree ID)
+    call identify_elements( nTreeID = cTreeID, tree = tree,               &
+      &                     pathFirst = pathFirst,                         &
+      &                     pathLast = pathLast,                           &
+      &                     levelDesc = levelDesc,                         &
+      &                     elemPos   = cPos,                          &
+      &                     prc = prc,                                     &
+      &                     computeStencil = computeStencil,               &
+      &                     nesting = parentNesting,                       &
+      &                     verbose = verbose )
+!    call tem_element_dump( me = levelDesc( coarserLevel )%elem, elemPos = cPos, &
+!       &         nUnit = dbgUnit, stencil = .true. )
+
+    if( cPos .le. 0 ) then
+      call tem_write(' Element not found ', int( cTreeID ))
+      call tem_write('   cPos         ', cPos )
+    end if
+
+    call identify_stencilNeigh( iElem = cPos, iLevel = coarserLevel,           &
+          &                     tree = tree, iStencil = 1,                     &
           &                     pathFirst = pathFirst,                         &
           &                     pathLast = pathLast,                           &
           &                     levelDesc = levelDesc,                         &
-          &                     elemPos   = elemPos,                           &
           &                     prc = prc,                                     &
           &                     computeStencil = computeStencil,               &
-          &                     nesting = nNesting,                            &
+          &                     nesting = parentNesting,                     &
           &                     verbose = verbose )
-        levelDesc( coarserLevel )%elem%stencil%val( neighPos )       &
-          &           %val(1)%totalPos( iStencilElem ) = elemPos
-        end if
-      end do
-      if( coarserLevel == tree%global%minlevel ) exit
-
-      nnTreeID = maxval( levelDesc( coarserLevel )                     &
-        &                 %elem                                        &
-        &                 %neighID%val(neighPos)                       &
-        &                         %val( levelDesc( coarserLevel )      &
-        &                               %elem%stencil%val( neighPos )  &
-        &                                            %val(1)%tIDpos(:) &
-        &                             )                                &
-        &              )
-      coarserLevel = coarserLevel - 1
-      nnTreeID = tem_parentOf(nnTreeID)
-
-      neighPos = PositionOfVal(levelDesc( coarserLevel )%elem%tID, nnTreeID )
-      if( neighPos .le. 0 ) then
-        call tem_write(' Element not found ', int( nnTreeID ))
-        call tem_write('   neighPos         ', neighPos )
-      end if
-
-      if ( levelDesc( coarserLevel )%elem%eType%val( neighPos ) &
-        &  .ne. eT_ghostFromCoarserElem ) exit
-    end do
 
   end subroutine create_allParentNeighbors
   !*****************************************************************************
   !!                     local process: identify if ghost or fluid
   subroutine single_process_neighbor( nTreeID, levelDesc, tree, prc, iProc, &
     &                                 minLevel, elemPos, foundLevel,        &
-    &                                 computeStencil, nesting )
+    &                                 computeStencil, nesting, updated )
     !---------------------------------------------------------------------------
     integer(kind=long_k), intent(in) :: nTreeID   !< neighboring treeID
     type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:)
     integer, intent(out) :: elemPos
     integer, intent(out) :: foundLevel
     integer, intent(in) :: nesting
+    logical, intent(out) :: updated !< was the element updated in this call?
     !---------------------------------------------------------------------------
     type(tem_stencilElement_type) :: emptyStencil(1)
     integer :: neighLevel
     neighLevel = tem_LevelOf( nTreeID ) !< Has to be same as tLevel!?
     if ( (neighLevel .lt. minLevel) &
       &  .or. (neighLevel .gt. uBound(levelDesc,1)) ) then
-      call tem_write(' Warning: level lower than minLevel was demanded in' &
-        &            // ' singleProcNeigh')
+      call tem_write(' Warning: level which is not included in the fluid tree ' &
+        &            // 'was demanded in singleProcNeigh')
       call tem_write('          treeID', nTreeID )
+      call tem_abort
     end if
+    ! Set the element updated flag as a default to false
+    updated = .false.
 
     ! If it is a remote cell on only one process -> regular halo
     if( iProc .ne. prc%rank + 1 )  then
         &        eType = eT_haloElem, property = 0_long_k,                     &
         &        sourceProc  = iProc, haloNesting = nesting,                   &
         &        stencil = emptyStencil, pos = elemPos, wasAdded = wasAdded )
+        write(dbgUnit,*) 'adding remote halo', nTreeID, 'nesting', nesting
       if( .not. wasAdded ) then
-        levelDesc( neighLevel )%elem%haloNesting%val( elemPos ) &
-          &  = min( levelDesc( neighLevel )%elem%haloNesting%val( elemPos ), &
-          &         nesting )
+        ! update the nesting to the current one
+        if( nesting .lt.  levelDesc( neighLevel )%elem%haloNesting%val( elemPos )) then
+          levelDesc( neighLevel )%elem%needsUpdate%val( elemPos ) = .true.
+          levelDesc( neighLevel )%elem%haloNesting%val( elemPos ) &
+            &  = min( levelDesc( neighLevel )%elem%haloNesting%val( elemPos ), &
+            &         nesting )
+          updated = .true.
+        end if
+      else
+        updated = .true.
       end if
+!      call tem_element_dump( me = levelDesc( neighLevel )%elem, elemPos = elemPos, &
+!      &  stencil = .true., nUnit = dbgUnit )
     else
       ! LOCAL
       ! Either a ghost or fluid cell
         &                          tree = tree, elemPos = elemPos,             &
         &                          level = foundLevel,                         &
         &                          minLevel = minLevel,                        &
+        &                          nesting = nesting,                          &
+        &                          updated = updated,                          &
         &                          computeStencil = computeStencil )
     end if
 
   !!    not existing( localPos=0): add to halo (MH: why to %halo and not
   !!    halosPrc(iProc)%mine ?)
   subroutine identify_local_element( nTreeID, levelDesc, tree,                 &
-    &              minLevel, level, elemPos, computeStencil, verbose)
+    &        minLevel, level, elemPos, nesting, updated, computeStencil, verbose)
     !---------------------------------------------------------------------------
     integer(kind=long_k), intent(in) :: nTreeID   !< neighboring treeID
     type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:)
     type(treelmesh_type), intent(in) :: tree
     integer, intent(in) :: minLevel
+    integer, intent(in) :: nesting
     type(tem_stencilHeader_type), intent(in) :: computeStencil
     integer, intent(out) :: elemPos
     integer, intent(out) :: level
+    logical, intent(out) :: updated !< was the element updated in this call?
     integer, intent(in), optional :: verbose !< verbosity level
     !---------------------------------------------------------------------------
     integer :: verbosity  !< local verbosity
       verbosity = chatter_none
     end if
 !    write(dbgUnit,*) ' identify_local_element: ', nTreeID
+    ! Set the element updated flag as a default to false
+    updated = .false.
 
     ! Position of neighbor treeID
     neighLevel = tem_LevelOf( nTreeID )
       foundLevel = tem_LevelOf( foundTID )
       if( foundLevel .eq. neighLevel ) then
         ! Fluid is already part of the elem list. Adding not necessary.
+        updated = .false.
       elseif( foundLevel .lt. neighLevel ) then
+        ! ghostFromCoarser
         ! existing fluid is coarser than current element
         ! Add all the children of found treeID down to me (intermediate levels)
         ! LBM specific. \todo hand over eligible child
         ! is added ghostFromCoarser
         call add_all_virtual_children( sourceID = foundTID,                   &
-          &                            foundPos = localPos,                   &
+          &               foundPos = localPos,                                &
+          &               elemPath = neighPath,                               &
           &               sourceProperty = tree%ElemPropertyBits( localPos ), &
-          &                            sourceLevel = foundLevel,              &
-          &                            targetLevel = neighLevel,              &
-          &                            levelDesc = levelDesc,                 &
-          &                            minlevel = minLevel,                   &
-          &                            tree = tree,                           &
-          &                            computeStencil = computeStencil )
+          &               sourceLevel = foundLevel,                           &
+          &               targetLevel = neighLevel,                           &
+          &               levelDesc = levelDesc,                              &
+          &               minlevel = minLevel,                                &
+          &               nesting = nesting,                                  &
+          &               updated = updated,                                  &
+          &               tree = tree,                                        &
+          &               computeStencil = computeStencil              )
       end if ! on same level?
     elseif( localPos .lt. 0) then
 !      write(dbgUnit,*) ' filling to real children for', nTreeID
         &                         minLevel = minlevel,                         &
         &                         tree = tree,                                 &
         &                         foundPos = dPos,                             &
+        &                         updated = updated,                           &
         &                         computeStencil = computeStencil )
     else
-      write(*,*) 'error, element not existing',nTreeID,                        &
+      write(*,*) 'Warning: element not existing',nTreeID,                        &
         &        'adding to nonexisting ...'
       call tem_tIDinfo( me = nTreeID, tree = tree )
       ! This case occurs, when a remote halo was added, which was not existing
   !> Find all the virtual children of the sourceTreeID down to the targetLevel
   !! and add to the level-wise ghostFromCoarser list in the level descriptor
   recursive subroutine add_all_virtual_children( &
-    &                    sourceID, sourceProperty, foundPos, sourceLevel, &
-    &                    targetLevel, levelDesc, minLevel, tree, &
-    &                    computeStencil )
+    &       sourceID, sourceProperty, foundPos, elemPath, sourceLevel, &
+    &       targetLevel, levelDesc, minLevel, tree, computeStencil,    &
+    &       nesting, updated )
     !--------------------------------------------------------------------------
     integer(kind=long_k), intent(in) :: sourceID
     integer(kind=long_k), intent(in) :: sourceProperty
+    type(tem_path_type), intent(in)  :: elemPath
+    integer,intent(in) ::  nesting
     integer,intent(in) ::  foundPos
     integer,intent(in) ::  sourceLevel
     integer,intent(in) ::  targetLevel
     integer, intent(in) :: minLevel
     type(treelmesh_type), intent(in) :: tree
     type(tem_stencilHeader_type), intent(in) :: computeStencil
+    logical, intent(out) :: updated !< was the element updated in this call?
     !--------------------------------------------------------------------------
     integer :: targetPos, iChild, iDir, nVals
     !> position of the existing (source) tID in the elem list
-    integer :: sourcePos
+    integer :: sourcePos, levelPos 
     logical :: wasAdded
-    integer(kind=long_k) :: children(8) ! child elements
+    logical :: childUpdated
+    integer(kind=long_k) :: cTreeID     ! current treeID to identify
     integer(kind=long_k) :: childID     ! child treeID
     integer(kind=long_k) :: curNeighborID, coarserNeighborID
     integer(kind=long_k), allocatable :: tNeighID(:)
     type(tem_stencilElement_type) :: tStencil(1)
     integer :: x(4), curLevel, offset(4), xc(4), xElem(4), neighIDpos, addedPos
     !--------------------------------------------------------------------------
-!    write(dbgUnit,*) ' add_virtual children from', sourceElem%tID, 'targetLevel', targetLevel
+!    write(dbgUnit,*) ' add_virtual children from', sourceID, 'targetLevel', targetLevel
     !Position of the coarser source element
     sourcePos = PositionOfVal( me = levelDesc( sourceLevel )%elem%tID, &
       &                        val = sourceID )
     allocate( tNeighID( computeStencil%QQN ))
     offset(4)   = 0
+    ! By default, set that no element was updated
+    updated = .false.
     if( sourceLevel .lt. targetLevel ) then
-      children = tem_directChildren( sourceID )
+      !children = tem_directChildren( sourceID )
       curLevel = sourceLevel + 1
       call init( me = tStencil(1), QQN = computeStencil%QQN, headerPos = 1 )
-      do iChild = 1, 8
+      !do iChild = 1, 8
         ! Add to the level-wise ghost list
-        xElem = tem_coordOfId( children( iChild ))
+        levelPos = targetLevel - sourceLevel
+        cTreeID = elemPath%node( levelPos )
+        ! determine which child it is
+        iChild = tem_childNumber( cTreeID )
+        xElem = tem_coordOfId( cTreeID ) !children( iChild ))
+!        write(dbgUnit,*) '      virtual: adding child ', cTreeID
         call append( me = levelDesc( curLevel )%elem,                          &
-          &          tID = children( iChild ),                                 &
+          &          tID = cTreeID,                                 &
           &          property = sourceProperty,                                &
           &          eType = eT_ghostFromCoarserElem,                          &
           &          sourceProc = tree%global%myPart+1,                        &
           &          wasAdded = wasAdded )
 
         if( wasAdded ) then
+          updated = .true.
           ! inherit boundaries and neighbors from parent
           x(:) = tem_coordOfId( int(iChild, kind=long_k) )
           ! inherit the boundary infos from the parent
           do iDir = 1, 3
             call tem_find_BCs_fromCoarser( dir = iDir, childCoord = x,         &
               & sourceLevel = sourceLevel, sourcePos = sourcePos,              &
-!              & targetLevel = curLevel, targetPos = targetPos,                   &
-!              & sourceElem = levelDesc( sourceLevel )%elem%tID%val( sourcePos ), &
-!              & targetElem = levelDesc( curLevel )%elem%tID%val( targetPos ),    &
               & neighID    = tNeighID,                                         &
               & minLevel = minLevel,                                           &
               & levelDesc = levelDesc,                                         &
           end if
         else
           ! Existing element encountered.
+          updated = .false.
         endif
         ! Overwrite the eventually existing nesting with the smallest value.
         ! The smallest nesting determines if further neighbors have to be
         ! retrieved in communicate_elements
-
+        if( nesting .lt. levelDesc( curLevel )%elem%haloNesting%val( targetPos )) then
+          ! needs update
+          updated = .true.
+          levelDesc( curLevel )%elem%needsUpdate%val( targetPos ) = .true.
+          levelDesc( curLevel )%elem%haloNesting%val( targetPos )                &
+            &       = min( nesting, levelDesc( curLevel )%elem%haloNesting%val( targetPos ))
+        end if
         ! In any case we have to recurse down to the target level
         ! lower levels might not yet exist.
-        call add_all_virtual_children( sourceID = children( iChild ),          &
+!        call tem_element_dump( me = levelDesc( curLevel )%elem, elemPos = targetPos, &
+!           &         nUnit = dbgUnit, stencil = .true. )
+
+        !MH: What about childUpdated?
+        call add_all_virtual_children( sourceID = cTreeID,          &
           &                            foundPos = targetPos,                   &
+          &                            elemPath = elemPath,                    &
+          &                            nesting = nesting,                      &
           &                            sourceLevel = curLevel,                 &
           &                            targetLevel = targetLevel,              &
           &                            levelDesc = levelDesc,                  &
           &                            minLevel = minLevel,                    &
+          &                            updated = childUpdated,                 &
           &                            tree = tree,                            &
           &                            sourceProperty = sourceProperty,        &
           &                            computeStencil = computeStencil )
-      end do
+      !end do
     end if
+    updated = ( updated .or. childUpdated )
   end subroutine add_all_virtual_children
   !*****************************************************************************
 
     type(tem_stencilElement_type) :: tStencil
     integer(kind=long_k) :: tNeighID
     !--------------------------------------------------------------------------
-!    write(dbgUnit,*) 'findBCs from finer targetElem', targetElem%tID
+    !write(dbgUnit,*) 'findBCs from finer targetElem', &
+    !& levelDesc( targetLevel )%elem%tID%val( targetPos )
     if  ( .not. allocated( levelDesc( targetLevel )%elem%stencil &
       &                                                 %val(targetPos) &
       &                                                 %val ) &
       ! link direction in the stencil. (only direct neighbors)
       ! Distinguishing three different kinds of neighbors:
       ! SIDE, EDGE and CORNER
-!      write(dbgUnit,*) 'current dir', iStencilElem
       linktype = sum( computeStencil%cxDir( :, iStencilElem )**2)
       select case(linktype)
       case(1)
           &             levelDesc = levelDesc( sourceLevel )                   &
           &             )
       end select
-!      write(dbgUNit,*) 'appending tNeighID', tNeighID
       ! Append the required treeID
       call append( me= levelDesc( targetLevel )%elem%neighID%val(targetPos), &
         &          val = tNeighID,    &
   !> Find all the children of the requested reqTreeID
   !! add to the level-wise ghostFromFine list in the level descriptor
   recursive subroutine fill_to_real_children( reqTreeID, levelDesc,            &
-    &                                         minLevel, tree,                  &
-    &                                         foundPos, computeStencil )
+    &       minLevel, tree, updated, foundPos, computeStencil )
     !--------------------------------------------------------------------------
     integer(kind=long_k), intent(in) :: reqTreeID   !< requested treeID
     type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:)
     integer, intent(in) :: minLevel
     type(treelmesh_type), intent(in) :: tree
     integer, intent(out) :: foundPos
+    logical, intent(out) :: updated !< was the current element updated in this call?
     type( tem_stencilHeader_type ), intent(in) :: computeStencil
     !--------------------------------------------------------------------------
     integer :: iChild, addedPos, reqLevel
     integer(kind=long_k) :: children(8), property
-    logical :: wasAdded
-    integer :: childPos(8)
+    logical :: wasAdded, childUpdated
+    integer :: childPos(8), neighLevel
+    type(tem_path_type) :: neighPath
     !--------------------------------------------------------------------------
 !    write(dbgUnit,*) 'filling to read children', reqTreeID
+    ! Set as not updated by default
+    updated = .false.
 
     ! Create the ghostFromFiner
     reqLevel        = tem_LevelOf( reqTreeID )
       &          pos = addedPos, wasAdded = wasAdded )
 
     if( wasAdded ) then
+      updated = .true.
       children = tem_directChildren( reqTreeID )
       childPos = 0 ! reset child positions. non-existing children are 0
       ! reset property
       property = 0_long_k
       do iChild = 1, 8
         ! What's the configuration of the requested element
-        childPos( iChild ) = PositionOfVal( levelDesc( reqLevel + 1 )%elem%tID,&
-                                            children(iChild))
+      !  childPos( iChild ) = PositionOfVal( levelDesc( reqLevel + 1 )%elem%tID,&
+      !                                      children(iChild))
+        neighLevel = reqLevel + 1
+        neighPath  = tem_PathOf( children( iChild ))
+
+        ! Return position in the treeIDlist
+        childPos( iChild ) = tem_PosOfPath( neighPath, tree%pathList )
+
+!        write(dbgUnit,*) ' childPos ', childPos( iChild ), 'ID',children( iChild )
         if( childPos( iChild ) .lt. 0 ) then
            ! finer elements! add them.
            call fill_to_real_children( reqTreeID = children( iChild ),         &
              &                         levelDesc = levelDesc,                  &
              &                         minLevel = minlevel,                    &
              &                         tree = tree,                            &
+             &                         updated = childUpdated,                 &
              &                         computeStencil = computeStencil,        &
              &                         foundPos = childPos( iChild ))
            ! Unify all properties of the children
            property = ieor( property, &
            &   levelDesc( reqLevel + 1)%elem%property%val( childPos( iChild )) )
+           updated = ( updated .or. childUpdated )
         else
-          ! Fluid cells are already in the element list
-          ! Not existing cells do not have to be added
+          ! Fluid cells are already in the element list (Not existing cells do not
+          ! have to be added). If the element is a fluid, we have to get its
+          ! position in the levelwise element list of the level descriptor.
+          childPos( iChild ) = PositionOfVal( levelDesc( reqLevel + 1 )%elem%tID,&
+                                            children(iChild))
         end if
       end do
+
       ! Now reconstruct current element's neighborhood
       ! based on the children's information
 !      write(dbgUnit,*) 'find_BCs_fromFiner'
     else
       ! ghostFromFiner element was not added
       ! Hence, we found available information which now has to be returned
+      updated = .false.
     end if
+
+!    call tem_element_dump( me = levelDesc( reqLevel )%elem, elemPos = addedPos, &
+!       &         nUnit = dbgUnit, stencil = .true. )
+
     ! Return the position of where the element was encountered
     foundPos = addedPos
 
     do iLevel = tree%global%minLevel, tree%global%maxLevel
       do iElem = 1, levelDesc( iLevel )%nElems_fluid
         do iStencil = 1, levelDesc( iLevel )%elem%stencil%val(iElem )%nVals
-!            call identify_stencilNeigh( iElem = iElem,                             &
-!              &                     elemStencil =  levelDesc( iLevel )%elem             &
-!            &                   %stencil%val( iElem )%val( iStencil )                            &
-!              &                     tree = tree,                               &
-!              &                     pathFirst = pathFirst,                     &
-!              &                     pathLast = pathLast,                       &
-!              &                     levelDesc = levelDesc,                     &
-!              &                     prc = prc,                                 &
-!              &                     computeStencil = computeStencil(1),        &
-!              &                     verbose = verbosity)
           do iNeighElem = 1, levelDesc( iLevel )%elem             &
             &                   %stencil%val( iElem )%val( iStencil )%QQN
             !MH: Same as in Line 2092 (identify_additionalNeigh)
             end if
             levelDesc( iLevel )%elem%stencil%val( iElem )%val( iStencil )      &
               &                %totalPos( iNeighElem ) = elemPos
-          end do ! iNeighElem
+          end do ! iNeighElem!
         end do ! iStencil
       end do ! iElem
     end do ! iLevel
 
     ! Find neighbors of neighbors from the require list
+    ! the require list
     do iLevel = tree%global%minLevel, tree%global%maxLevel
       do iElem = 1, levelDesc( iLevel )%require%nVals
         iStencil = 1
         elemPos = PositionOfVal( me = levelDesc( iLevel )%elem%tID,            &
           &                      val = levelDesc( iLevel )%require%val( iElem ))
-        if( elemPos .gt. 0 ) then
+        if( elemPos .gt. 0 .and. levelDesc(iLevel )%elem%eType%val( elemPos ) .gt. 0) then
          do iNeighElem = 1, levelDesc( iLevel )%elem             &
            &                   %stencil%val( elemPos )%val( iStencil )%QQN
             neighPos =levelDesc( iLevel )%elem%stencil%val( iElem )            &
   !*****************************************************************************
 
 
+ !*****************************************************************************
+  !> Return the element type of a treeID by defining the position in the totalList
+  function tem_eTypeOfTotalPos( pos, levelDesc ) result( eType )
+    !--------------------------------------------------------------------------
+    integer, intent(in) :: pos !< the element you are looking for
+    !> the descriptor you use for searching
+    type(tem_levelDesc_type), intent(in) :: levelDesc
+    integer :: eType
+    !--------------------------------------------------------------------------
+    integer :: levelPos
+    !--------------------------------------------------------------------------
+    eType   = eT_undefinedElem
+    levelPos = PositionOfVal( me = levelDesc%elem%tID,           &
+      &                       val = levelDesc%total( pos ))
+    if( levelPos .gt. 0 ) then
+      eType = levelDesc%elem%eType%val( levelPos )
+    end if
+
+  end function tem_eTypeOfTotalPos
   !*****************************************************************************
-  !> Returns the absolute position in the total list of a given treeID
-  !! opposed to PosOfId, where the relative position in one of the separate
-  !! lists is returned. Herefore, total list has to be created beforehand.
+
+  !*****************************************************************************
+  !> Return the element type of a treeID 
   function tem_eTypeOfId( tID, levelDesc ) result( eType )
     !--------------------------------------------------------------------------
     integer(kind=long_k), intent(in) :: tID !< the element you are looking for
 
 
   !*****************************************************************************
-  !> 
-  !!
+  !> identify additionally required neighbor elements
+  !! run over the 'require' list of elements, which was accumulated before in
+  !! init_elemLevels. The list includes neighbor elements of stencil neighbors,
+  !! for stencils with the requireNeighNeigh attribute set. 
+  !! This is needed for example for LBM boundary stencil elements, which in turn
+  !! require their compute stencil neighborhood to allow PULL operations from
+  !! there
   subroutine identify_additionalNeigh( tree, prc, levelDesc,                       &
     &                  pathFirst, pathLast, computeStencil, verbose )
     !--------------------------------------------------------------------------
     type(tem_levelDesc_type), intent(inout) :: levelDesc(tree%global           &
       &                                                      %minlevel:)
     type(tem_path_type), intent(in) :: pathFirst(:), pathLast(:)
+    !> the compute stencil, for which the additional neighbors are reconstructed
     type(tem_stencilHeader_type), intent(in) :: computeStencil
     integer, intent(in) :: verbose
     !--------------------------------------------------------------------------
     integer :: nesting, verbosity
     integer(kind=long_k) :: nTreeID
     !--------------------------------------------------------------------------
+    ! The position of the compute stencil
     iStencil = 1
     do iLevel = tree%global%minlevel, tree%global%maxLevel
+      ! Run over the additionally required element list
       do indElem = 1, levelDesc( iLevel )%require%nVals 
+        ! get the position of the treeID in the element list
         iElem = PositionOfVal( me = levelDesc( iLevel )%elem%tID,            &
           &                      val = levelDesc( iLevel )%require%val( indElem ))
-        if( iElem .gt. 0 ) then
+        ! get the element position 
+        if( iElem .gt. 0 .and. levelDesc( iLevel )%elem%eType%val( iElem ) .gt. 0) then
+          ! Run over all the neighbors of the compute stencil
           do iNeighElem = 1, levelDesc( iLevel )%elem%stencil%val( iElem )   &
               &                                   %val( iStencil )%QQN
+            ! position
             neighPos = levelDesc( iLevel )%elem%stencil%val( iElem )           &
               &                  %val( iStencil )%tIDpos( iNeighElem ) 
             if( neighPos .gt. 0 ) then
               &                %totalPos( iNeighElem ) = neighPos
           end do
         else
-          call tem_write('not found required element in mesh ')
+          call tem_write('not found required element in mesh ', nTreeID)
         end if
       end do
     end do
   !! return type of added element in levelPos(2)
   !! Also, non-existing elements are reported as such (levelPos(2))
   subroutine identify_halo( haloTreeID, elemPos, haloLevel, levelDesc, tree, &
-    &                       minLevel, computeStencil )
+    &            updated,   nesting, minLevel, computeStencil )
     !--------------------------------------------------------------------------
     type(treelmesh_type), intent(in) :: tree
     integer, intent(in) :: minLevel
+    integer, intent(in) :: nesting
     integer(kind=long_k), intent(in) :: haloTreeID   !< neighboring treeID
     integer, intent(out) :: elemPos !< type and position in list of found treeID
     integer, intent(out) :: haloLevel
+    logical, intent(out) :: updated  
     type(tem_levelDesc_type), intent(inout) :: levelDesc(minlevel:)
     type(tem_stencilHeader_type), intent(in) :: computeStencil
     !--------------------------------------------------------------------------
         &                          levelDesc = levelDesc,                      &
         &                          tree = tree, elemPos = elemPos,             &
         &                          level = haloLevel,                          &