tem_init_elemLevels Subroutine

public subroutine tem_init_elemLevels(me, boundary, tree, stencils)

subroutine to find neighbours of cells

Typically every element requires information from its neighbors to perform an update of its state. All such required neighbors constitute a so called tem_stencil_module. tem_build_horizontalDependencies, or the connectivity of elements on the same refinement level in the tree, are therefore essential for the stencils in typical mesh-based numerical schemes. We now sketch the method to find the connectivity of elements in the mesh with the help of the described mesh layout. In a distributed mesh, the first distinction has to be made with respect to processor ownership.

Todo

tem_find_depNeigh does not exist any more Is the element in question tem_find_depNeigh a local or remote element?

For all identify_local_element, their actual position in the sparse mesh has to be identified for a given treeID.

A tem_path_type comparison of two nodes to decide their relative position in the ordered list of elements has to take into account the hierarchy of the mesh.

As in the actual simulation only leaves will be present, and no overlap of refinement levels is allowed, this is a sufficient determination. However, such overlapping regions might be created to enable the interpolation between different levels with virtual elements. These special elements are distinguished anyway, so this does not pose any problem in the search for neighborhood relations.

An element can be either identified by its treeID or by a tuple with four integer entries for the spatial coordinates and the refinement level: . This coordinate fully describes the spatial shape and position of an element and is important to determine spatial relations between elements. Conversion between treeIDs and coordinates can be achieved by the routines - tem_CoordOfId Coordinate from TreeID - tem_IdOfCoord TreeID from Coordinate

The spatial indices are limited by the refinement level: . To avoid undefined coordinates, movements through the mesh by additions to indices are done in a modulo() safeguard resulting in an periodic universe cube. This enables the movement in the mesh in horizontal direction by index alteration and translation into a treeID again. The described procedure is completely reversible, enabling the construction of the treeID for any given coordinate. Thus, the conversion between this coordinate and the serialized treeID encoding is fully described by the topology of the octree and the chosen space-filling curve ordering of elements, as explained above.

With this method we can describe the neighborhood of any given element independently of the actual solver by a simple list of relative offsets. In contrast to the generic horizontal relation, the vertical relation between child and parent nodes in the tree requires an interpolation operator between different refinement levels. This interpolation usually has to take into account solver specific requirements, but is otherwise quite isolated from the numerical operation on each refinement level. The TreElM library offers the solver a level-wise view, as suggested by the properties described above. To find all required neighbors in the distributed octree, the solver merely has to provide its horizontal dependencies. These are described with the help of an element specific tem_stencil_module. A stencil is basically a set of element-offsets , describing the relative positions of all required elements for a given element.

In this routine, level descriptor is allocated, all elements in tree are added into me%elem as fluid type, including their property, pntTID, stencil, neighID, sourceProc

Arguments

Type IntentOptional Attributes Name
type(tem_levelDesc_type), intent(out), allocatable :: me(:)

neighbor list containing all the neighbours for the cells given in treeidsubset. Result of this routine
type(tem_BC_prop_type), intent(in) :: boundary

boundaries for the elements with bnd property set
type(treelmesh_type), intent(in) :: tree

subset of tree ids for which the neighbours will be specified
type(tem_stencilHeader_type), intent(in) :: stencils(:)

the given stencil

Calls

proc~~tem_init_elemlevels~~CallsGraph proc~tem_init_elemlevels tem_init_elemLevels interface~append~11 append proc~tem_init_elemlevels->interface~append~11 interface~init~24 init proc~tem_init_elemlevels->interface~init~24 interface~tem_logging_isactive tem_logging_isActive proc~tem_init_elemlevels->interface~tem_logging_isactive proc~tem_abort tem_abort proc~tem_init_elemlevels->proc~tem_abort proc~tem_alloc_leveldesc tem_alloc_levelDesc proc~tem_init_elemlevels->proc~tem_alloc_leveldesc proc~tem_calc_neighbors tem_calc_neighbors proc~tem_init_elemlevels->proc~tem_calc_neighbors proc~tem_coordofid tem_CoordOfId proc~tem_init_elemlevels->proc~tem_coordofid proc~tem_elemlist_dump tem_elemList_dump proc~tem_init_elemlevels->proc~tem_elemlist_dump proc~tem_horizontalspacer tem_horizontalSpacer proc~tem_init_elemlevels->proc~tem_horizontalspacer proc~tem_require_dump tem_require_dump proc~tem_init_elemlevels->proc~tem_require_dump proc~tem_stencil_map_totreelmdef tem_stencil_map_toTreelmDef proc~tem_init_elemlevels->proc~tem_stencil_map_totreelmdef proc~append_ga_dynlong append_ga_dynlong interface~append~11->proc~append_ga_dynlong proc~append_ga_dynlong_vec append_ga_dynlong_vec interface~append~11->proc~append_ga_dynlong_vec proc~init_ga2d_real init_ga2d_real interface~init~24->proc~init_ga2d_real proc~tem_logging_isactive_for tem_logging_isActive_for interface~tem_logging_isactive->proc~tem_logging_isactive_for proc~tem_logging_isactive_primary tem_logging_isActive_primary interface~tem_logging_isactive->proc~tem_logging_isactive_primary mpi_abort mpi_abort proc~tem_abort->mpi_abort proc~tem_alloc_leveldesc->interface~init~24 proc~tem_firstidatlevel tem_FirstIdAtLevel proc~tem_calc_neighbors->proc~tem_firstidatlevel proc~tem_idofcoord tem_IdOfCoord proc~tem_calc_neighbors->proc~tem_idofcoord proc~tem_levelof tem_LevelOf proc~tem_coordofid->proc~tem_levelof proc~tem_element_dump tem_element_dump proc~tem_elemlist_dump->proc~tem_element_dump interface~expand~9 expand proc~append_ga_dynlong->interface~expand~9 proc~append_ga_dynlong_vec->interface~expand~9 interface~tem_stencil_dump tem_stencil_dump proc~tem_element_dump->interface~tem_stencil_dump proc~expand_ga_dynlong expand_ga_dynlong interface~expand~9->proc~expand_ga_dynlong proc~tem_stencilelement_dump tem_stencilElement_dump interface~tem_stencil_dump->proc~tem_stencilelement_dump proc~tem_stencilheader_dump tem_stencilHeader_dump interface~tem_stencil_dump->proc~tem_stencilheader_dump

Called by

proc~~tem_init_elemlevels~~CalledByGraph proc~tem_init_elemlevels tem_init_elemLevels proc~tem_create_leveldesc tem_create_levelDesc proc~tem_create_leveldesc->proc~tem_init_elemlevels proc~tem_dimbydim_construction tem_dimByDim_construction proc~tem_dimbydim_construction->proc~tem_create_leveldesc proc~tem_build_face_info tem_build_face_info proc~tem_build_face_info->proc~tem_dimbydim_construction

Source Code

  subroutine tem_init_elemLevels( me, boundary, tree, stencils )
    ! -------------------------------------------------------------------- !
    !> neighbor list containing all the neighbours for the
    !! cells given in treeidsubset. Result of this routine
    type(tem_levelDesc_type), allocatable, intent(out)  :: me(:)
    !> boundaries for the elements with bnd property set
    type(tem_bc_prop_type), intent(in)                  :: boundary
    !> subset of tree ids for which the neighbours will be specified
    type(treelmesh_type), intent(in)                    :: tree
    !> the given stencil
    type(tem_stencilHeader_type), intent(in)            :: stencils(:)
    ! -------------------------------------------------------------------- !
    type(tem_stencilElement_type) :: tStencil
    integer :: posInTree, nElemsBnd, iQQN, iLevel, nProcs, hashpos
    integer :: x(4), nStencils, iStencil, elemPos, nStencilElems
    integer :: indElem, minLevel, maxLevel, QQN
    integer :: initlen
    integer(kind=long_k) :: treeID
    integer(kind=long_k), allocatable :: neighIDs(:)
    integer :: addedPos
    logical :: wasAdded
    integer :: posInBCID
    ! -------------------------------------------------------------------- !

    call tem_horizontalSpacer( fUnit = dbgUnit(1), before = 1 )
    write(dbgUnit(1),*) 'Inside routine: tem_init_elemLevels '

    ! precondition for this routine: the first stencil must be compute stencil
    if ( .not. stencils(1)%useAll ) then
      write(logUnit(1),*)'tem_init_elemLevels:'
      write(logUnit(1),*)'Error: The first stencil does not have the'//    &
        &                'useAll flag set'
      write(logUnit(1),*)'       The first one is considered the compute ' &
        &            //'stencil and'
      write(logUnit(1),*)'       needs to be set for all elements.'
      call tem_abort()
    end if

    nProcs   = tree%global%nParts
    minLevel = tree%global%minLevel
    maxLevel = tree%global%maxLevel
    nElemsBnd = 0
    nStencils = size( stencils )

    ! -------------   Prepare a hash table   -----------------------------------
    ! Find a suitable hash size.
    ! The hash lookup is O(1), while the identification grows with
    ! log(nElems) and log(nProcs), yet allocation, deallocation and
    ! filling will cause some overhead of O(nHashes), this is an attempt to
    ! find a compromise of these parameters.
    nHashes = int(max(tree%nElems, nProcs, 64), kind=long_k)

    ! @todo Find the closest smaller prime for the hash size, to minimize the
    !       number of collisions?
    ! Limit the size of the hash:
    nHashes = min(nHashes, int(io_buffer_size, kind=long_k))

    allocate(hash(0:int(nHashes-1)))
    allocate(hash_elemPos(0:int(nHashes-1)))
    ! Reset the hash
    hash = -1_long_k
    hash_elemPos = -1
    ! --------------------------------------------------------------------------


    initlen = ceiling(1.1 * real(tree%nElems)/real(maxlevel-minlevel+1))
    call tem_alloc_levelDesc( me, minLevel, maxLevel, initlen, nStencils )

    write(logUnit(5),*) 'Searching the neighbors for each element ...'

    ! Loop over all given stencils
    stencilLoop: do iStencil = 1, nStencils
      write(logUnit(6), "('   on stencil: ', I0)") iStencil
      QQN = stencils(iStencil)%QQN
      allocate( neighIDs( QQN ) )

      ! map user stencil to treelm definition
      call tem_stencil_map_toTreelmDef( me = stencils( iStencil ))

      if ( stencils(iStencil)%useAll ) then
        ! This stencil belongs to all elements in the tree (fluid stencil)
        nStencilElems = tree%nElems
      else
        ! This stencil belongs to only a subset of the tree
        ! (boundary and IBM stencils)
        nStencilElems = stencils( iStencil )%nElems
      end if

      write(logUnit(6),"(A,I0)") '   nElems: ', nStencilElems
      write(logUnit(6),"(A,I0)") '   QQN: ', QQN

      ! Loop over all the elements connected to the current stencil
      ! stencils( iStencil )
      elemLoop: do indElem = 1, nStencilElems
        ! Reset the empty stencil
        ! Create the temporary stencil for assigning the positions of the
        ! neighbor treeIDs which are added in elem%neighID
        ! also store the original stencil it depends on (depends = iStencil )
        call init( me        = tStencil, &
          &        QQN       = QQN,      &
          &        headerPos = iStencil )

        if( stencils( iStencil )%useAll ) then
          ! if the stencil is defined to be connected to all elements,
          ! then simply loop over all the fluid elements
          posInTree = indElem
        else
          ! if the stencil is defined for a subset, loop over the
          ! defined subset in %elem
          posInTree = stencils( iStencil )%elem%val( indElem )
        end if
        ! assign fluid element with info from the tree list and properties
        ! from disk
        treeID = tree%treeID( posInTree )
        ! Get topological info about coordinates+level
        x = tem_CoordOfId( treeID )

        ! Find position of this element in the elemList as it might be
        ! there already from previous stencil iterations
        hashpos = int(mod( treeID, nHashes))
        if (hash(hashpos) /= treeID) then ! cache miss
          call append( me         = me( x(4) )%elem,                  &
            &          tID        = treeID,                           &
            &          pntTID     = posInTree,                        &
            &          eType      = eT_fluid,                         &
            &          nNeighIDs  = QQN,                              &
            &          sourceProc = tree%global%myPart+1,             &
            &          property   = tree%ElemPropertyBits(posInTree), &
            &          pos        = elemPos                           )
          ! for differentiating between fluid and ghost nodes
          me( x(4) )%elem%property%val(elemPos) & 
            & = ibset( me( x(4) )%elem%property%val(elemPos), prp_fluid)
          !write(dbgunit(1),*) "prp elems = ", me( x(4) )%elem%property%val(elemPos)
          !flush(dbgUnit(1))
          !stop
          hash(hashpos) = treeID
          hash_elemPos( hashpos ) = elemPos
        else ! cache hit
          elemPos = hash_elemPos( hashpos )
        end if

        posInBCID = -1
        ! Only perform the boundary check based on the treeID boundary
        ! information for stencils which work on all elements (so we can
        ! actually read out the boundary information from the tree )
        if ( iStencil == 1 ) then
          ! First check, if I have prp_hasBnd
          ! If I don't take all the elements, then I can't analyze boundaries
          if( btest( me(x(4))%elem%property%val( elemPos ), prp_hasBnd )) then
            ! If yes, there might be no neighbors
            nElemsBnd = nElemsBnd +1
            posInBCID = nElemsBnd
          end if
        end if ! useAll elements? (boundary information)

        ! calculate possible neighbors or get BCID
        call tem_calc_neighbors( posInBCID   = posInBCID,            &
          &                      boundary_ID = boundary%boundary_ID, &
          &                      stencil     = stencils(iStencil),   &
          &                      x           = x,                    &
          &                      neighIDs    = neighIDs              )

        ! append neighbors into elem%neighID or require list
        do iQQN = 1, QQN

          call append( me       = me(x(4))%elem%neighID%val( elemPos ), &
            &          val      = neighIDs(iQQN),                       &
            &          pos      = addedPos,                             &
            &          wasAdded = wasAdded )
          tStencil%tIDpos(iQQN) = addedPos

          if ( stencils(iStencil)%requireNeighNeigh .and. neighIDs(iQQN)>0 ) then
            call append( me       = me(x(4))%require, &
              &          val      = neighIDs(iQQN),   &
              &          pos      = addedPos,         &
              &          wasAdded = wasAdded )
          end if

        end do

        ! Append the temporary stencil to the element
        call append( me  = me(x(4))%elem%stencil%val( elemPos ),             &
          &          val = tStencil )
      end do elemLoop ! indElem = 1, nStencilElems

      deallocate( neighIDs )
    end do stencilLoop ! iStencil = 1, nStencils
    write(logUnit(5),*) 'Finished searching the neighbors ID for each fluid.'

    if (tem_logging_isActive(main_debug%logger, 7)) then
      do iLevel = minLevel, maxLevel
        write(dbgUnit(1),"(A,I0)") 'Dump levelDesc%elem on level: ', iLevel
        call tem_elemList_dump( me      = me( iLevel )%elem,                  &
          &                     compact = .true.,                             &
          &                     nUnit   = dbgUnit(5),                         &
          &                     stencil = .true.,                             &
          &                     string  = 'after initializing level elements' &
          &                               //' i.e. only fluids')
        if( me( iLevel )%require%nVals > 0 ) then
          write(dbgUnit(1),"(A,I0)") 'Dump levelDesc%require on level: ', iLevel
          call tem_require_dump( me      = me( iLevel )%require,               &
            &                    nUnit   = dbgUnit(5),                         &
            &                    string  = 'after initializing level elements' )
        end if
      end do
    end if

    write(dbgUnit(1),*) 'Leave  routine: tem_init_elemLevels'
    call tem_horizontalSpacer( fUnit = dbgUnit(1), after = 1 )

  end subroutine tem_init_elemLevels