! Copyright (c) 2011-2018,2021 Harald Klimach <harald.klimach@dlr.de> ! Copyright (c) 2011-2012 Manuel Hasert <m.hasert@grs-sim.de> ! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de> ! Copyright (c) 2011 Metin Cakircali <m.cakircali@grs-sim.de> ! Copyright (c) 2011 Aravindh Krishnamoorthy <aravindh28.4@gmail.com> ! Copyright (c) 2011-2012 Khaled Ibrahim <k.ibrahim@grs-sim.de> ! Copyright (c) 2011, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de> ! Copyright (c) 2011-2016, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de> ! Copyright (c) 2012 Jens Zudrop <j.zudrop@grs-sim.de> ! Copyright (c) 2013 Melven Zoellner <yameta@freenet.de> ! Copyright (c) 2013-2015 Kartik Jain <kartik.jain@uni-siegen.de> ! Copyright (c) 2014, 2017 Verena Krupp <verena.krupp@uni-siegen.de> ! Copyright (c) 2014 Monika Harlacher <monika.harlacher@uni-siegen.de> ! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de> ! Copyright (c) 2016, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de> ! Copyright (c) 2016 Neda Ebrahimi Pour <neda.epour@uni-siegen.de> ! Copyright (c) 2017 Daniel Petró <daniel.petro@student.uni-siegen.de> ! Copyright (c) 2017 Raphael Haupt <Raphael.Haupt@student.uni-siegen.de> ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this ! list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, ! this list of conditions and the following disclaimer in the documentation ! and/or other materials provided with the distribution. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE ! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE ! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR ! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, ! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! **************************************************************************** ! ! !> TREe based ELemental Mesh module, providing the definition and methods to !! act on a octree based mesh. !! !! The mesh may either be externally generated by seeder or some simple mesh !! configurations can also be created internally. !! In case of an external mesh from seeder, the directory containing the mesh !! needs to be provided in the configuration: !! !!```lua !! mesh = '/path/to/mesh/' !!``` !! !! Note that the string provided here is used as a prefix. Thus, if the mesh !! is in a directory it needs to end in a path seperator. !! The default is 'mesh/', which results in the application to look for the !! mesh in the `mesh` subdirectory in the working directory of the application. !! !! Internal meshes are defined in a table, where the first option describes !! which of the predefined meshes is to be generated, then there needs to be !! a definition of the universe (or root) cube within which the octree lives. !! This is given by an `origin` and then the `length` by which it extends along !! each axis. Further options may be required by the respective predefined mesh !! types. !! An example for a fully periodic uniform cubical mesh: !! !!```lua !! mesh = { !! predefined = 'cube', !! origin = {0., 0., 0.}, !! length = 10., !! refinementLevel = 4 !! } !!``` !! !! The following predefined meshes are available: !! !! * 'cube' create a fully periodic cube without boundary conditions and uniform !! element sizes by refining the root element (universe). For this kind of !! mesh you need to provide an `refinementLevel` which defines how often the !! encompassing cube should be bisected in each dimension leading to meshes !! with `8^refinementLevel` elements. !! * 'slice' create a slice of the fully periodic cube with uniform elements. !! For the slice you need to provide the `refinementLevel`, as for the cube. !! The resulting mesh will consist of `4^refinementLevel` elements. !! * 'line' creates a periodic line. For this kind of mesh it is possible to !! specify either an `refinementLevel` to create meshes with !! `2^refinementLevel` elements, or the number of elements in the mesh by !! setting `element_count`. !! * 'line_bounded' create a line with boundary conditions in west (on the left) !! and east (on the right). For this mesh the length needs to be defined by !! `element_count`. !! !! The mesh is stored in [[treelmesh_type]]. !! The definition is designed to be used for parallel executions. !! It includes a reference to [[tem_global_module:tem_global_type]] !! and the number of elements, which were loaded on the local partition, as well !! as the first and last tree ID of each partition in the current communicator. !! This enables the determination of any tree ID's position among the processes, !! independent of their actual existence. !! This is important to identify the location of neighbor elements when creating !! connectivity lists for the linearized element lists locally on each process !! within the solvers. !! !!  !! !! See also !! --- !! !! - [[tem_balance_module]] "Load Balancing" !! - [[tem_dyn_array_module]] "Dynamic Data Structures" !! module treelmesh_module ! include treelm modules use mpi use env_module, only: rk, & & long_k, & & long_k_mpi, & & rk_mpi, & & PathLen, & & tem_create_EndianSuffix, & & globalMaxLevels, & & io_buffer_size use tem_aux_module, only: tem_open, tem_abort, check_mpi_error use tem_global_module, only: tem_global_type, & & load_tem_global, & & dump_tem_global, & & tem_global_mesh_free use tem_property_module, only: tem_property_type, & & gather_property, & & prp_hasbnd, & & prp_sendHalo, & & prp_hasNormal, & & prp_hasQVal use tem_topology_module, only: tem_path_type, & & tem_firstIdAtLevel, & & tem_lastIdAtLevel, & & tem_levelOf, & & tem_IDofCoord use tem_logging_module, only: logUnit use tem_tools_module, only: tem_horizontalSpacer use tem_sparta_module, only: tem_balance_sparta, & & tem_sparta_type, & & tem_init_sparta, & & tem_destroy_sparta, & & tem_exchange_sparta ! include aotus modules use aotus_module, only: aot_get_val, & & flu_state, & & aoterr_NonExistent, & & aoterr_WrongType, & & aoterr_Fatal use aot_table_module, only: aot_table_open, & & aot_table_close, & & aot_table_length implicit none private public :: load_tem public :: unload_treelmesh public :: treelmesh_type, load_treelmesh, dump_treelmesh public :: free_treelmesh public :: generate_treelm_cube public :: treelmesh_fromList public :: tem_load_internal public :: dump_MeshHeader public :: exchange_elements public :: tem_dump_weights public :: tem_load_weights !> declaration of the overall mesh on the local partition type treelmesh_type !> This entry provides some global informations on the mesh, which is needed !! across all partitions. type(tem_global_type) :: global ! @todo: need update after mesh changes integer :: nElems !< Total number of Elements on this partition integer(kind=long_k) :: elemOffset !< Offset of this partition !> List of treeIDs of the first element on each partition, has a length of !! global%nParts integer(kind=long_k), allocatable :: Part_First(:) !> List of treeIDs of the last element on each partition, has a length of !! global%nParts integer(kind=long_k), allocatable :: Part_Last(:) !> The treeIDs of each element (the array has a length of nElems) !! these IDs identify the element in the complete mesh, with a breadth-first !! numbering scheme. Allowing easy computation of parents and children. !! [treeid] \ref treeID !! !! # TreeID {#treeID} !! The coordinate description of an element has three spatial entries and a !! fourth entry indicating the refinement level in the octree, resulting in !! a tuple of four integers: \f$(x, y, z, L)\f$. !! !! This coordinate fully describes the spatial shape and position of an !! element and is important to determine spatial relations between elements. !! It can be obtained from a given treeID by the following procedure. !! 1. First, the \ref tem_topology_module::tem_levelof "refinement level" is !! found. !! 2. Then a \ref tem_topology_module::tem_firstidatlevel "level-offset" is !! calculated. !! 3. After the offset is known, the Morton index of the node in question is !! obtained by subtracting the offset \f$s\f$ from the treeID of the !! node. !! 4. Finally a simple bit interleaving rule allows the !! \ref tem_topology_module::tem_coordofid "transformation from the !! Morton index to the three dimensional coordinates". !! !! The spatial indices are limited by the refinement level: !! \f$x, y, z \in \mathbf{Z}_{\ge 0}: x, y, z < 2^L\f$. !! To avoid undefined coordinates, movements through the mesh by additions !! to indices are done in a modulo(\f$2^L\f$) safeguard resulting in a !! periodic universe cube. !! This enables the movement in the mesh in horizontal direction (e.g. !! for stencil neighbor identification) by index alteration and translation !! into a treeID again. !! @todo: change to dynamic array integer(kind=long_k), allocatable :: treeID(:) !> Bit field of element properties for each element (the array has a length !! of nElems), each bit in the 8 byte integer indicates the presence or !! absence of a given property, which are further described in the !! Property component. !! @todo: change to dynamic array integer(kind=long_k), allocatable :: ElemPropertyBits(:) !> declaration of additional elemental properties, the array has a length of !! global%nProperties, each property provides a list of local element !! indices of elements with this property. type(tem_property_type), pointer :: Property(:) => null() !> Path for all elements in the treeID array type(tem_path_type), allocatable :: pathList(:) !> Name of the file that contains weight information for each element to !! allow load balancing. !! If empty string, balancing will be deactivated. character(len=pathLen) :: weights_file = '' !> File to write weights to after simulation. !! Defaults to the same name as weights_file, if you do not want to have !! the input file overwritten, you can provide a different file to write !! to by setting write_weights accordingly in the configuration. character(len=pathLen) :: write_weights !> Level weights as provided by the solver. real(kind=rk) :: levelWeight(globalMaxLevels) end type contains ! ************************************************************************ ! !> Load the treelmesh !! !! Depending on the configuration given in conf, this will either load a !! a \ref tem_distributed_octree "tree mesh" from disk with the given prefix !! or \ref tem_load_internal "create the mesh internally". !! @todo: to calculate levelWeight, global has to be loaded fisrt, !! maybe load_global shoule be moved out of this routine subroutine load_tem( me, conf, myPart, nParts, comm, levelWeight, meshDir, & & parent) ! -------------------------------------------------------------------- ! !> Structure to load the mesh to type(treelmesh_type), intent(out) :: me !> Directory containing the mesh informations type(flu_State) :: conf !> Partition to use on the calling process (= MPI Rank in comm) integer, intent(in) :: myPart !> Number of partitions, the mesh is partitioned into (= Number of MPI !! processes in comm). integer, intent(in) :: nParts !> MPI Communicator to use !! now defaults to the one given in the tree%global%comm integer, intent(in) :: comm !> Balancing weight for elements on different levels. !! If these weights are present, the mesh is partitioned initially according !! to the given weights. !! This can not be combined with given offsets. If offset is also provided !! they will overwrite the levelWeight distribution. real(kind=rk), intent(in), optional :: levelWeight(globalMaxLevels) !> output Mesh directory name character(len=*), intent(out), optional :: meshDir !> optional parent handle integer, intent(in), optional :: parent ! -------------------------------------------------------------------- ! integer :: tem_handle character(len=pathLen) :: dirname character(len=4) :: EndianSuffix integer :: mesh_error, weights_error ! The communicator which is used integer :: commLocal real(kind=rk), allocatable :: weights(:) integer(kind=long_k) :: chunksize, remainder logical :: found_weights type( tem_sparta_type ) :: sparta ! -------------------------------------------------------------------- ! ! Use the incoming communicator commLocal = comm call tem_horizontalSpacer(fUnit = logUnit(1)) EndianSuffix = tem_create_EndianSuffix() ! Make sure the dirname has some sane value. dirname = '' call aot_table_open(L = conf, & & parent = parent, & & thandle = tem_handle, & & key = 'mesh' ) if (tem_handle /= 0) then ! key mesh is a table ! The mesh is actually given as a table, parse it and ! generate a mesh internally. write(logUnit(1),*) 'Generating an internally defined mesh' call tem_load_internal(me, conf, tem_handle, myPart, nParts, commLocal) else ! tem_handle == 0, mesh is not a table ! The mesh is not a table, try to interpret it as a string. call aot_get_val( L = conf, & & thandle = parent, & & key = 'mesh', & & val = dirname, & & ErrCode = mesh_error, & & default = 'mesh/' ) if (btest(mesh_error, aoterr_Fatal)) then write(logUnit(0),*) 'FATAL Error occured, while retrieving mesh folder!' if (btest(mesh_error, aoterr_NonExistent))then write(logUnit(0),*) 'Variable not existent!' write(logUnit(0),*) 'Using the default value mesh folder: mesh/' else if (btest(mesh_error, aoterr_WrongType))then write(logUnit(0),*) 'Variable has wrong type!' write(logUnit(0),*) 'STOPPING' call tem_abort() end if end if ! Load global information ! @todo: this routine may already be called outside if (associated(me%global%property)) deallocate(me%global%property) call load_tem_global(me%global, dirname, myPart, nParts, commLocal) ! Homogeneous distribution chunksize = me%global%nElems / nParts remainder = mod(me%global%nElems, int(nParts, kind=long_k)) me%nElems = int(chunksize) if (myPart < int(remainder)) then me%nElems = me%nElems + 1 end if me%elemOffset = chunksize * int(myPart, kind=long_k) & & + min(int(myPart, kind=long_k), remainder) call load_treelmesh( me = me, & & nParts = nParts ) end if ! If requested, provide the directory name to the caller. if ( present(meshDir) ) then meshDir = trim( dirname ) end if ! Make sure each process has at least one element! if (me%global%nElems/me%global%nParts == 0) then write(logunit(1),*) '' write(logunit(1),*) 'FATAL Error occured:' write(logunit(1),*) 'Less elements then processes, this is not' write(logunit(1),*) 'supported! And you should reconfigure your' write(logunit(1),*) 'setup.' write(logunit(1),*) '' write(logunit(1),*) 'Restart this simulation with at most ', & & me%global%nElems, ' processes!' call tem_abort() end if me%levelWeight = 0.0_rk if (present(levelWeight)) then me%levelWeight = levelWeight end if ! Get the name for the weights file. ! This defaults to the mesh prefix with weights appended. call aot_get_val( L = conf, & & key = 'weights', & & thandle = parent, & & val = me%weights_file, & & ErrCode = weights_error, & & default = '' ) if (.not. btest(weights_error, aoterr_NonExistent)) then if (myPart == 0 .and. trim(me%weights_file) /= '') then inquire(file=trim(me%weights_file)//EndianSuffix, exist=found_weights) if (.not. found_weights) then write(logUnit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(logUnit(1),*) 'Weight file ' & & // trim(me%weights_file)//EndianSuffix & & // ' can NOT be found!' write(logUnit(1),*) 'May proceed WITHOUT BALANCING!' write(logUnit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' end if end if end if call aot_get_val( L = conf, & & key = 'write_weights', & & thandle = parent, & & val = me%write_weights, & & ErrCode = weights_error, & & default = me%weights_file ) if (trim(me%write_weights) /= '') then write(logunit(1),*) 'The computational weight of each element will be' write(logunit(1),*) 'written at the end of the simulation into file' write(logunit(1),*) trim(me%write_weights)//EndianSuffix write(logunit(1),*) 'If this is not desired, set write_weights explicitly' write(logunit(1),*) 'to an empty string.' end if ! Do sparta balancing, but only if there are actually multiple processes. if (me%global%nParts > 1) then allocate( weights(me%nElems) ) call tem_load_weights( me = me, & & weights = weights, & & success = found_weights ) if ( found_weights ) then call tem_init_sparta( sparta, nParts ) call tem_balance_Sparta(weight = weights, & & myPart = myPart, & & nParts = nParts, & & comm = commLocal, & & myElems = me%nElems, & & sparta = sparta, & & offset = me%elemOffset ) call exchange_elements( me, sparta ) call tem_destroy_sparta( sparta ) end if deallocate( weights ) end if ! Output basic mesh infos to the screen call dump_meshHeader( me ) call aot_table_close(L=conf, thandle=tem_handle) call tem_horizontalSpacer(fUnit = logUnit(1)) end subroutine load_tem ! ************************************************************************ ! ! ************************************************************************ ! !> Unload the treelmesh !! !! Deallocate the data within tree, while keep the global information !! subroutine unload_treelmesh( me ) ! -------------------------------------------------------------------- ! !> Structure to load the mesh to type(treelmesh_type), intent(inout) :: me ! -------------------------------------------------------------------- ! write(logUnit(6),*) 'Deallocate the tree data structure' if (allocated(me%Part_first)) deallocate( me%Part_First ) if (allocated(me%Part_last)) deallocate( me%Part_Last ) if (allocated(me%treeID)) deallocate( me%treeID ) if (allocated(me%ElemPropertyBits)) deallocate( me%ElemPropertyBits ) if (associated(me%property)) deallocate(me%property) if (allocated(me%pathList)) deallocate( me%pathList ) end subroutine unload_treelmesh ! ************************************************************************ ! ! ************************************************************************ ! !> Completely free the memory of the treelmesh again. subroutine free_treelmesh(me) !> Structure to free type(treelmesh_type), intent(inout) :: me call unload_treelmesh(me) call tem_global_mesh_free(me%global) end subroutine free_treelmesh ! ************************************************************************ ! ! ************************************************************************ ! !> Load an internally generated mesh !! !! The predefined meshes are illustrated in the following configuration !! snippet: !! !!```lua !! -- use the predefined full cube !! mesh = { predefined = 'cube', !! -- for a single element in Z-direction use slice: !! -- predefined = 'slice', !! !! -- for just a line (single element in Y and Z) use line: !! -- predefined = 'line', !! -- If you set the refinementLevel to 0, a single element will !! -- be created with treeID=1 on level 1, with the bounding cube !! -- redefined such, that the element fills the original !! -- specification for the bounding cube. !! -- The line can also be generated with an arbitrary number of !! -- elements instead of a level. In this case the length refers !! -- to the overall length of the given elements not the extent of !! -- the cube. Periodicity is assumed in X direction as well. !! -- To create the line with a certain number of elements set the !! -- element_count instead of the refinementLevel. If both are !! -- given, the element_count overrules the refinementLevel. !! !! -- for a line with boundary conditions line_bounded: !! -- predefined = 'line_bounded' !! -- Like line, but needs to have `element_count` given, not !! -- `refinementLevel`. This mesh defines boundary conditions in the !! -- X direction as west (on the left) and east (on the right), that !! -- then need to be defined in the solver configuration. !! -- origin of the cube !! origin = {0.,0.,0.}, !! -- length of the cube !! length = 10., !! -- refinement level to resolve the cube !! refinementLevel = 4 } !!``` !! subroutine tem_load_internal( me, conf, thandle, myPart, nParts, comm ) ! -------------------------------------------------------------------- ! !> Structure to load the mesh to type(treelmesh_type), intent(out) :: me !> Directory containing the mesh informations type(flu_State) :: conf !> Handle for the table to read the description !! of the mesh from. integer, intent(in) :: thandle !> Partition to use on the calling process (= MPI Rank in comm) integer, intent(in) :: myPart !> Number of partitions, the mesh is partitioned into (= Number of MPI !! processes in comm). integer, intent(in) :: nParts !> MPI Communicator to use integer, intent(in) :: comm ! -------------------------------------------------------------------- ! character(len=40) :: meshtype real(kind=rk) :: origin(3) real(kind=rk) :: length integer :: elementcount integer :: level integer :: sub_handle integer :: iError, i integer :: orig_err(3) ! -------------------------------------------------------------------- ! call aot_get_val( L = conf, & & thandle = thandle, & & val = meshtype, & & ErrCode = iError, & & key = 'predefined', & & default = 'cube' ) select case(trim(meshtype)) case('cube') write(logUnit(1),*)'Creating a full cubical mesh without boundaries' ! Get the origin of the cube: call aot_table_open( L = conf, & & parent = thandle, & & thandle = sub_handle, & & key = 'origin' ) if (aot_table_length(L = conf, thandle = sub_handle) == 3) then do i=1,3 call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(i), & & ErrCode = iError, & & pos = i ) end do else write(logUnit(1),*)'Origin specified for the mesh to generate has to ' & & //'have 3 entries!' write(logUnit(1),*)'Default to {0., 0., 0.} ' origin = 0.0_rk end if call aot_table_close( L=conf, thandle=sub_handle ) ! Get the length of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = length, & & ErrCode = iError, & & key = 'length', & & default = 1.0_rk ) ! Get the refinement level: call aot_get_val( L = conf, & & thandle = thandle, & & val = level, & & ErrCode = iError, & & key = 'refinementLevel' ) ! Now generate a mesh with the gathered configuration parameters if (level > 0) then ! Only if the level is greater than 0, actually create mesh. call generate_treelm_cube( me = me, & & origin = origin, & & length = length, & & level = level, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) else ! Meshes on level 0 have only one element, but we need treeids > 0, ! thus we create single elements on level 1 with treeid=1, and ! redefine the bounding cube accordingly. call generate_treelm_single( me = me, & & origin = origin, & & length = 2*length, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) end if case('slice') write(logUnit(1),*)'Creating a slice (single element in Z) mesh' & & //' without boundaries' ! Get the origin of the cube: origin = 0.0_rk call aot_table_open( L = conf, & & parent = thandle, & & thandle = sub_handle, & & key = 'origin' ) if (aot_table_length(L = conf, thandle = sub_handle) >= 2) then do i=1,2 call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(i), & & ErrCode = iError, & & pos = i ) end do else write(logUnit(1),*)'Origin specified for the mesh to generate has to ' & & // 'have 2 entries!' write(logUnit(1),*)'Default to {0., 0.} ' end if call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(3), & & ErrCode = iError, & & Default = 0.0_rk, & & pos = 3 ) call aot_table_close( L=conf, thandle=sub_handle ) ! Get the length of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = length, & & ErrCode = iError, & & key = 'length', & & default = 1.0_rk ) ! Get the refinement level: call aot_get_val( L = conf, & & thandle = thandle, & & val = level, & & ErrCode = iError, & & key = 'refinementLevel' ) ! Now generate a mesh with the gathered configuration parameters if (level > 0) then ! Only if the level is greater than 0, actually create mesh. call generate_treelm_slice( me = me, & & origin = origin, & & length = length, & & level = level, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) else ! Meshes on level 0 have only one element, but we need treeids > 0, ! thus we create single elements on level 1 with treeid=1, and ! redefine the bounding cube accordingly. call generate_treelm_single( me = me, & & origin = origin, & & length = 2*length, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) end if case('line') write(logUnit(1),*) 'Creating a line (single element in Y and Z) mesh' & & // ' without boundaries' ! Get the origin of the cube: origin = 0.0_rk call aot_table_open( L = conf, & & parent = thandle, & & thandle = sub_handle, & & key = 'origin' ) if (aot_table_length(L = conf, thandle = sub_handle) >= 1) then call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(1), & & ErrCode = iError, & & pos = 1 ) else write(logUnit(1),*) 'Origin specified for the mesh to generate has' & & // ' to have 1 entry!' write(logUnit(1),*)'Default to {0.} ' end if call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(2), & & ErrCode = iError, & & Default = 0.0_rk, & & pos = 2 ) call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(3), & & ErrCode = iError, & & Default = 0.0_rk, & & pos = 3 ) call aot_table_close( L=conf, thandle=sub_handle ) ! Get the length of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = length, & & ErrCode = iError, & & key = 'length', & & default = 1.0_rk ) ! Get the refinement level: call aot_get_val( L = conf, & & thandle = thandle, & & val = level, & & ErrCode = iError, & & key = 'refinementLevel', & & default = -1 ) ! Get the element count: call aot_get_val( L = conf, & & thandle = thandle, & & val = elementcount, & & ErrCode = iError, & & key = 'element_count', & & default = -1 ) if ( (elementcount < 1) .and. (level < 0) ) then write(logUnit(1),*) 'For the line mesh you need to specify either' write(logUnit(1),*) 'the element_count or the refinementLevel!' write(logUnit(1),*) 'STOPPING...' call tem_abort() end if ! Now generate a mesh with the gathered configuration parameters if (elementcount > 1) then call generate_treelm_elements( me = me, & & origin = origin, & & length = length, & & elementcount = elementcount, & & myPart = myPart, & & nParts = nParts, & & comm = comm, & & predefined = trim(meshtype), & & bclabel = 'internal 1D BC' ) else if ((level > 0) .and. (elementcount < 1)) then ! Only if the level is greater than 0, actually create mesh. call generate_treelm_line( me = me, & & origin = origin, & & length = length, & & level = level, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) else ! Meshes on level 0 have only one element, but we need treeids > 0, ! thus we create single elements on level 1 with treeid=1, and ! redefine the bounding cube accordingly. call generate_treelm_single( me = me, & & origin = origin, & & length = 2*length, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) end if end if case('line_bounded') write(logUnit(1),*) 'Creating a line (single element in Y and Z) mesh' & & // ' WITH boundaries' write(logUnit(1),*) 'Boundaries west and east have to be provided!' ! Get the origin of the cube: origin = 0.0_rk call aot_table_open( L = conf, & & parent = thandle, & & thandle = sub_handle, & & key = 'origin' ) if (aot_table_length(L = conf, thandle = sub_handle) >= 1) then call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(1), & & ErrCode = iError, & & pos = 1 ) else write(logUnit(1),*) 'Origin specified for the mesh to generate has' & & // ' to have 1 entry!' write(logUnit(1),*)'Default to {0.} ' end if call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(2), & & ErrCode = iError, & & Default = 0.0_rk, & & pos = 2 ) call aot_get_val( L = conf, & & thandle = sub_handle, & & val = origin(3), & & ErrCode = iError, & & Default = 0.0_rk, & & pos = 3 ) call aot_table_close( L=conf, thandle=sub_handle ) ! Get the length of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = length, & & ErrCode = iError, & & key = 'length', & & default = 1.0_rk ) ! Get the element count: call aot_get_val( L = conf, & & thandle = thandle, & & val = elementcount, & & ErrCode = iError, & & key = 'element_count', & & default = -1 ) if ( (elementcount < 1) ) then write(logUnit(1),*) 'For the bounded line mesh you need to specify' write(logUnit(1),*) 'the element_count!' write(logUnit(1),*) 'STOPPING...' call tem_abort() end if ! Now generate a mesh with the gathered configuration parameters call generate_treelm_elements( me = me, & & origin = origin, & & length = length, & & elementcount = elementcount, & & myPart = myPart, & & nParts = nParts, & & comm = comm, & & predefined = trim(meshtype), & & bclabel = 'bounded internal 1D BC' ) case('single') write(logUnit(1),*) 'Creating a single element mesh without boundaries' ! Get the origin of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = origin, & & ErrCode = orig_err, & & key = 'origin', & & default = [0.0_rk, 0.0_rk, 0.0_rk] ) ! Get the length of the cube: call aot_get_val( L = conf, & & thandle = thandle, & & val = length, & & ErrCode = iError, & & key = 'length', & & default = 1.0_rk ) ! Now generate a mesh with the gathered configuration parameters call generate_treelm_single( me = me, & & origin = origin, & & length = length, & & myPart = myPart, & & nParts = nParts, & & comm = comm ) case('default') write(logUnit(1),*) 'Do not know how to generate mesh ' & & // trim(meshtype) end select end subroutine tem_load_internal ! ************************************************************************ ! ! ************************************************************************ ! !> This subroutine reads a mesh in treelm format from disk at the !! specified directory name. !! subroutine load_treelmesh(me, nParts ) ! -------------------------------------------------------------------- ! !> Structure to load the mesh to type(treelmesh_type), intent(inout) :: me !> Directory containing the mesh informations ! character(len=*), intent(in) :: dirname !> Partition to use on the calling process (= MPI Rank in comm) ! integer, intent(in) :: myPart !> Number of partitions, the mesh is partitioned into (= Number of MPI !! processes in comm). integer, intent(in) :: nParts ! -------------------------------------------------------------------- ! logical :: ex integer :: iElem integer :: iProp integer :: fh, etype, ftype integer :: iostatus( MPI_STATUS_SIZE ) integer :: iError integer :: file_rec_len integer :: typesize character(len=300) :: ElemFileName character(len=4) :: EndianSuffix integer(kind=long_k), allocatable :: buffer(:) integer(kind=long_k), allocatable :: globbuffer(:) integer(kind=MPI_OFFSET_KIND) :: displacement ! -------------------------------------------------------------------- ! EndianSuffix = tem_create_EndianSuffix() ElemFileName = trim(me%global%dirname)//'elemlist'//EndianSuffix write(logUnit(1), *) 'Load mesh from file: '//trim(ElemFileName) if (associated(me%property)) deallocate(me%property) allocate(me%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) allocate(me%treeID(me%nElems)) if (me%global%myPart == 0) then ! Only one process should check the file existence inquire(file=trim(ElemFileName), exist=ex) if (.not. ex) then write(logUnit(1),*) 'ERROR: File ' // trim(ElemFileName) & & // ' not found!' write(logUnit(1),*) 'Aborting.' call tem_abort() end if end if allocate(me%ElemPropertyBits(me%nElems)) allocate(buffer(me%nElems*2)) if (me%global%nElems*2 > io_buffer_size) then write(logUnit(1),*) 'using MPI to read' ! Create a contiguous type to describe the vector per element call MPI_Type_contiguous( 2, long_k_mpi, etype, iError ) call check_mpi_error(iError,'type etype in load_treelmesh') call MPI_Type_commit( etype, iError ) call check_mpi_error(iError,'commit etype in load_treelmesh') call MPI_Type_size(etype, typesize, iError ) call check_mpi_error(iError,'typesize in load_treelmesh') ! Create a MPI Subarray as ftype for file view call MPI_Type_contiguous( me%nElems, etype, ftype, iError ) call check_mpi_error(iError,'type ftype in load_treelmesh') call MPI_Type_commit( ftype, iError ) call check_mpi_error(iError,'commit ftype in load_treelmesh') ! Calculate displacement for file view displacement = me%elemOffset * typesize * 1_MPI_OFFSET_KIND ! Set the view for each process on the file above ! Open the binary file for MPI I/O (Write) call MPI_File_open( me%global%comm, trim(ElemFileName), & & MPI_MODE_RDONLY, MPI_INFO_NULL, fh, iError ) call check_mpi_error(iError,'file_open in load_treelmesh') call MPI_File_set_view( fh, displacement, etype, ftype, "native", & & MPI_INFO_NULL, iError ) call check_mpi_error(iError,'set_view in load_treelmesh') ! Read data from the file call MPI_File_read_all( fh, buffer, me%nElems, etype, iostatus, iError ) call check_mpi_error(iError,'read_all in load_treelmesh') !Free the MPI_Datatypes which were created and close the file call MPI_Type_free (etype, iError) call check_mpi_error(iError,'free etype in load_treelmesh') call MPI_Type_free (ftype, iError) call check_mpi_error(iError,'free ftype in load_treelmesh') call MPI_File_close(fh, iError) call check_mpi_error(iError,'close file in load_treelmesh') ! END IO-part else !> The mesh is so small, it probably is better to read it on one process !! and distribute the resulting data via the network. write(logUnit(1),*) 'using one process to read and broadcast' allocate( globbuffer(me%global%nElems*2) ) if (me%global%myPart == 0) then !! single proc read inquire(iolength = file_rec_len) globbuffer call tem_open( newunit = fh, & & file = trim(ElemFileName), & & recl = file_rec_len, & & action = 'read', & & access = 'direct', & & form = 'unformatted' ) read(fh, rec=1) globbuffer close(fh) end if !! single proc read call MPI_Bcast( globbuffer, int(2*me%global%nElems), long_k_mpi, 0, & & me%global%comm, iError ) buffer = globbuffer(me%elemOffset*2+1:(me%ElemOffset+me%nElems)*2) deallocate(globbuffer) end if ! END IO-part ! Fill the tree with what was read in buffer from disk do iElem = 1, me%nElems me%treeID(iElem) = buffer( (iElem-1)*2 + 1 ) me%ElemPropertyBits(iElem) = buffer( (iElem-1)*2 + 2 ) end do deallocate( buffer ) ! clear the prp_sendHalo property since the distribution of elements ! might have changed do iElem=1,me%nElems me%ElemPropertyBits(iElem) = ibclr( me%ElemPropertyBits(iElem), & & prp_sendHalo ) end do call MPI_allgather( me%treeID(1), 1, long_k_mpi, & & me%Part_First, 1, long_k_mpi, & & me%global%comm, iError ) call MPI_allgather( me%treeID(me%nElems), 1, long_k_mpi, & & me%Part_Last, 1, long_k_mpi, & & me%global%comm, iError ) do iProp=1,me%global%nProperties call gather_Property( Property = me%Property(iProp), & & Header = me%global%Property(iProp), & & BitField = me%ElemPropertyBits, & & comm = me%global%comm ) end do write(logUnit(1),*) 'Done, reading the mesh!' end subroutine load_treelmesh ! ************************************************************************ ! ! ************************************************************************ ! !> Write a given mesh to disk. it is stored to the directory given in the !! [[tem_global_module:tem_global_type]]. !! !! Dump treelmesh_type%global to header.lua !! Dump treeID and propertyBits to elemlist.lsb (little endian environments) !! or elemlist.msb (big endian environments) !! subroutine dump_treelmesh( me, root_only ) ! -------------------------------------------------------------------- ! !> Mesh to dump to disk type(treelmesh_type), intent(inout) :: me !> root dump global mesh when true and !! all process dump its own mesh when false logical, intent(in), optional :: root_only ! -------------------------------------------------------------------- ! integer(kind=long_k), allocatable :: buffer(:) integer :: fh, etype, ftype integer :: iostatus( MPI_STATUS_SIZE ) integer :: iElem integer :: iError integer :: typesize logical :: root_not_participating character(len=300) :: ElemFileName character(len=4) :: EndianSuffix integer(kind=long_k) :: file_offset integer(kind=MPI_OFFSET_KIND) :: displacement ! -------------------------------------------------------------------- ! if(present(root_only)) then root_not_participating = root_only else root_not_participating = .true. endif EndianSuffix = tem_create_EndianSuffix() ! Dump global info to header.lua ! if root_not_participating = .false. that each process dumps its own mesh ! set global%nElems to local%nElems to dump local nElems ! in header file of each process if (.not. root_not_participating) me%global%nElems = me%nElems call dump_tem_global(me%global, root_only) ElemFileName = trim(me%global%dirname) // 'elemlist' // EndianSuffix file_offset = 0 allocate( buffer(me%nElems * 2) ) ! Get the offset for MPI File view if(root_not_participating) then call MPI_Exscan( int(me%nElems,long_k), file_offset, 1, & & long_k_mpi, MPI_SUM, me%global%comm, iError) else file_offset = 0 endif ! Fill the buffer with data to be written on disk do iElem = 1, me%nElems buffer((iElem-1)*2+1) = me%treeID(iElem) buffer((iElem-1)*2+2) = me%ElemPropertyBits(iElem) end do write(logUnit(1),*) 'using MPI to write' ! Open the binary file for MPI I/O (Write) call MPI_File_open( me%global%comm, & & trim(ElemFileName), & & ior(MPI_MODE_WRONLY,MPI_MODE_CREATE), & & MPI_INFO_NULL, & & fh, iError ) ! Catch if there was an exception by MPI call check_mpi_error(iError,'file_open in dump_treelmesh') ! Create a contiguous type to describe the vector per element call MPI_Type_contiguous( 2, long_k_mpi, etype, iError ) call check_mpi_error(iError,'type etype in dump_treelmesh') call MPI_Type_commit( etype, iError ) ! Catch MPI Exception(s) call check_mpi_error(iError,'commit etype in dump_treelmesh') !get size of etype call MPI_Type_size(etype, typesize, iError ) call check_mpi_error(iError,'typesize in dump_treelmesh') ! Create a MPI Contiguous as ftype for file view call MPI_Type_contiguous( me%nElems, etype, ftype, iError ) call check_mpi_error(iError,'type ftype in dump_treelmesh') call MPI_Type_commit( ftype, iError ) ! Catch MPI Exception(s) call check_mpi_error(iError,'commit ftype in dump_treelmesh') ! calculate dsplacement for file view displacement = file_offset * typesize * 1_MPI_OFFSET_KIND ! Set the view for each process on the file above call MPI_File_set_view( fh, displacement, etype, ftype, "native", & & MPI_INFO_NULL, iError ) ! Catch MPI Exception(s) call check_mpi_error(iError,'set_view in dump_treelmesh') ! Write the Data to File call MPI_File_write_all( fh, buffer, me%nElems, etype, iostatus, & & iError ) ! Catch MPI Exception(s) call check_mpi_error(iError,'write_all in dump_treelmesh') !Free the MPI_Datatypes which were created and close the file call MPI_Type_free(etype, iError) call check_mpi_error(iError,'free etype in dump_treelmesh') call MPI_Type_free(ftype, iError) call check_mpi_error(iError,'free ftype in dump_treelmesh') call MPI_File_close(fh, iError) call check_mpi_error(iError,'close file in dump_treelmesh') deallocate( buffer ) end subroutine dump_treelmesh ! ************************************************************************ ! ! ************************************************************************ ! !> summary: Generate a predefined cube !! This serves as an simple grid generation for performance or scaling !! analysis without being obliged to use Seeder. !! You have to specify the generic grid parameters in the lua file insted of !! the mesh folder !!```lua !! mesh = { predefined='cube', !! origin = {0.,0.,0.}, !! length = 10., !! refinementLevel = 6 } !!``` !! You have to specify the shape 'cube', a bounding box origin, its length !! and also the refinement level, which results in different amount elements !! in the grid. !! The result of this routine is mainly the treeID list with the additional !! lists for saving the properties. They are all set to zero here, however. !! As we only have a simple cube which includes all the elements on this !! level, the treeID list just contains contiguously increasing integers !! subroutine generate_treelm_cube( me, origin, length, level, myPart, nParts, & & comm ) ! -------------------------------------------------------------------- ! !> Mesh to generate type(treelmesh_type), intent(out) :: me !> Corner of the cube real(kind=rk), intent(in) :: origin(3) !> Length of cube real(kind=rk), intent(in) :: length !> Resolution level integer, intent(in) :: level !> Partition of the caller (starts with 0) integer, intent(in) :: myPart !> Number of partitions integer, intent(in) :: nParts !> communicator to be used integer, intent(in) :: comm ! -------------------------------------------------------------------- ! integer(kind=long_k) :: firstID, lastID integer(kind=long_k) :: share integer :: remainder integer :: iPart, iElem ! -------------------------------------------------------------------- ! me%global%nParts = nParts me%global%myPart = myPart me%global%comm = comm me%global%origin = origin me%global%BoundingCubeLength = length me%global%minLevel = level me%global%maxLevel = level me%global%label = 'Generic_Cube' me%global%predefined = 'cube' write(me%global%comment,'(a15,i7,a16,i2,a1)') & & 'Generated with ', nParts, ' Parts on Level ', level, '.' me%global%dirname = './' ! No properties in this mesh me%global%nProperties = 0 if (associated(me%global%Property)) deallocate(me%global%Property) allocate(me%global%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) ! Compute the treeIDs of the mesh: firstID = tem_firstIdAtLevel(level) lastID = tem_lastIdAtLevel(level) ! Total number of elements in this mesh me%global%nElems = lastID - firstID + 1_long_k share = me%global%nElems / int(nParts, kind=long_k) remainder = int(mod(me%global%nElems, int(nParts, kind=long_k))) ! The first partition starts always with the firstID me%Part_First(1) = firstID ! Up to remainder partitions have share + 1 elements do iPart=2,remainder+1 me%Part_Last(iPart-1) = me%Part_First(iPart-1) + share me%Part_First(iPart) = me%Part_Last(iPart-1) + 1 end do ! The remaining elements get exactly the share elements: do iPart=remainder+2,nParts me%Part_Last(iPart-1) = me%Part_First(iPart-1) + share - 1 me%Part_First(iPart) = me%Part_Last(iPart-1) + 1 end do ! The last partition ends always with the lastID me%Part_Last(nParts) = lastID ! Local data: me%nElems = int(me%Part_Last(myPart+1) - me%Part_First(myPart+1)) + 1 me%elemOffset = me%Part_first(myPart+1) - firstID allocate(me%treeID(me%nElems)) allocate(me%ElemPropertyBits(me%nElems)) if (associated(me%property)) deallocate(me%property) allocate(me%Property(me%global%nProperties)) ! No Properties: me%ElemPropertyBits = 0_long_k ! Filling the treeIDs: do iElem=1,me%nElems me%treeID(iElem) = me%Part_First(myPart+1) + int(iElem - 1, kind=long_k) end do end subroutine generate_treelm_cube ! ************************************************************************ ! ! ************************************************************************ ! function tem_CoordOf_2d_Id(id2d, level) result(coord) ! -------------------------------------------------------------------- ! !> input element ID integer(kind=long_k), intent(in) :: Id2d integer, intent(in) :: level !> coordinate of element return value integer :: coord(4) ! -------------------------------------------------------------------- ! integer(kind=long_k) :: fak(2) integer :: bitlevel integer :: i integer(kind=long_k) :: telem ! -------------------------------------------------------------------- ! coord(:) = 0 coord(4) = level fak(1) = 1 fak(2) = 2 bitlevel = 1 tElem = id2d ! get x coordinate from ! ! x = sum(iLevel=0, inf) { 2**iLevel * mod(tElem / (4**iLevel) ),2) } ! y = sum(iLevel=0, inf) { 2**iLevel * mod(tElem /(2 * (4**iLevel) ),2) } ! do while ((tElem / fak(1)) > 0) do i=1,2 coord(i) = coord(i) + bitlevel * int(mod(tElem / fak(i), 2_long_k)) end do bitlevel = bitlevel * 2 fak = fak * 4 end do end function tem_CoordOf_2d_Id ! ************************************************************************ ! ! ************************************************************************ ! !> summary: Generate a predefined slice !! This serves as an simple grid generation for performance or scaling !! analysis without being obliged to use Seeder. !! You have to specify the generic grid parameters in the lua file instead of !! the mesh folder !!```lua !! mesh = { predefined='slice', !! origin = {0.,0.,0.}, !! length = 10., !! refinementLevel = 6 } !!``` !! You have to specify the shape 'slice', a bounding box origin, its length !! and also the refinement level, which results in different amount elements !! in the grid. !! The result of this routine is mainly the treeID list with the additional !! lists for saving the properties. !! The generated slice will be all elements in the first layer (XY plane) in !! Z-direction, with periodicity in Z direction. !! subroutine generate_treelm_slice( me, origin, length, level, myPart, & & nParts, comm ) ! -------------------------------------------------------------------- ! !> Mesh to generate type(treelmesh_type), intent(out) :: me !> Corner of the cube real(kind=rk), intent(in) :: origin(3) !> Length of cube real(kind=rk), intent(in) :: length !> Resolution level integer, intent(in) :: level !> Partition of the caller (starts with 0) integer, intent(in) :: myPart !> Number of partitions integer, intent(in) :: nParts !> communicator to be used integer, intent(in) :: comm ! -------------------------------------------------------------------- ! integer(kind=long_k) :: firstID, lastID integer(kind=long_k) :: share integer(kind=long_k) :: previous integer :: remainder integer :: iPart, iElem integer :: coord(4) integer :: lastcoord ! -------------------------------------------------------------------- ! me%global%nParts = nParts me%global%myPart = myPart me%global%comm = comm me%global%origin = origin me%global%BoundingCubeLength = length me%global%minLevel = level me%global%maxLevel = level me%global%label = 'Generic_Slice' me%global%predefined = 'slice' write(me%global%comment,'(a15,i7,a16,i2,a1)') & & 'Generated with ', nParts, ' Parts on Level ', level, '.' me%global%dirname = './' ! Boundary property to define periodic boundary me%global%nProperties = 1 if (associated(me%global%Property)) deallocate(me%global%property) if (associated(me%Property)) deallocate(me%property) allocate(me%global%Property(me%global%nProperties)) allocate(me%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) ! Compute the treeIDs of the mesh: firstID = tem_firstIdAtLevel(level) lastcoord = 2**level - 1 lastID = tem_IdOfCoord( coord = [lastcoord, lastcoord, 0, level], & & offset = firstID ) ! Total number of elements in this mesh me%global%nElems = 4_long_k**level share = me%global%nElems / int(nParts, kind=long_k) remainder = int(mod(me%global%nElems, int(nParts, kind=long_k))) ! The first partition starts always with the firstID me%Part_First(1) = firstID ! Up to remainder partitions have share + 1 elements previous = 0_long_k do iPart=2,remainder+1 previous = previous + int(share) coord = tem_coordof_2d_id(previous, level) me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) previous = previous + 1 coord = tem_coordof_2d_id(previous, level) me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The remaining elements get exactly the share elements: do iPart=remainder+2,nParts previous = previous + int(share) - 1 coord = tem_coordof_2d_id(previous, level) me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) previous = previous + 1 coord = tem_coordof_2d_id(previous, level) me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The last partition ends always with the lastID me%Part_Last(nParts) = lastID ! Local data: if (myPart < remainder) then me%nElems = int(share+1) me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k) else me%nElems = int(share) me%elemOffset = (int(myPart, kind=long_k) * share) & & + int(remainder, kind=long_k) end if ! All elements have (periodic) boundaries me%Property(1)%nElems = me%nElems me%Property(1)%offset = me%elemOffset allocate(me%Property(1)%ElemID(me%Property(1)%nElems)) ! Please note, that the tem_bc_prop_module will set boundary conditions ! based on this label accordingly! me%global%Property(1)%label = 'internal 2D BC' me%global%Property(1)%bitpos = prp_hasBnd me%global%Property(1)%nElems = me%Property(1)%nElems allocate(me%treeID(me%nElems)) allocate(me%ElemPropertyBits(me%nElems)) ! No Properties: me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd) ! Filling the treeIDs: do iElem=1,me%nElems me%Property(1)%ElemID(iElem) = iElem coord = tem_coordof_2d_id(me%elemOffset + iElem - 1, level) me%treeID(iElem) = tem_idofcoord(coord, offset=firstID) end do end subroutine generate_treelm_slice ! ************************************************************************ ! ! ************************************************************************ ! !> summary: Generate a predefined line !> This serves as an simple grid generation for performance or scaling !! analysis without being obliged to use Seeder. !! You have to specify the generic grid parameters in the lua file instead of !! the mesh folder !! ~~~~~~~~~~~~~~~~~~~~~{.lua} !! mesh = { predefined='line', !! origin = {0.,0.,0.}, !! length = 10., !! refinementLevel = 6 } !! ~~~~~~~~~~~~~~~~~~~~~ !! \n !! You have to specify the shape 'line', a bounding box origin, its length !! and also the refinement level, which results in different amount elements !! in the grid.\n !! The result of this routine is mainly the treeID list with the additional !! lists for saving the properties. !! The generated line will be all elements along the X-Axis !! with periodicity in Y and Z direction. !! subroutine generate_treelm_line( me, origin, length, level, myPart, & & nParts, comm ) ! -------------------------------------------------------------------- ! !> Mesh to generate type(treelmesh_type), intent(out) :: me !> Corner of the cube real(kind=rk), intent(in) :: origin(3) !> Length of cube real(kind=rk), intent(in) :: length !> Resolution level integer, intent(in) :: level !> Partition of the caller (starts with 0) integer, intent(in) :: myPart !> Number of partitions integer, intent(in) :: nParts !> communicator to be used integer, intent(in) :: comm ! -------------------------------------------------------------------- ! integer(kind=long_k) :: firstID, lastID integer(kind=long_k) :: share integer :: remainder integer :: iPart, iElem integer :: coord(4) integer :: lastcoord ! -------------------------------------------------------------------- ! me%global%nParts = nParts me%global%myPart = myPart me%global%comm = comm me%global%origin = origin me%global%BoundingCubeLength = length me%global%minLevel = level me%global%maxLevel = level me%global%label = 'Generic_Line' me%global%predefined = 'line' write(me%global%comment,'(a15,i7,a16,i2,a1)') & & 'Generated with ', nParts, ' Parts on Level ', level, '.' me%global%dirname = './' ! Boundary property to define periodic boundary me%global%nProperties = 1 if (associated(me%global%property)) deallocate(me%global%property) if (associated(me%property)) deallocate(me%property) allocate(me%global%Property(me%global%nProperties)) allocate(me%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) ! Compute the treeIDs of the mesh: firstID = tem_firstIdAtLevel(level) lastcoord = 2**level - 1 lastID = tem_IdOfCoord( coord = [lastcoord, 0, 0, level], offset = firstID ) ! Total number of elements in this mesh me%global%nElems = 2_long_k**level share = me%global%nElems / int(nParts, kind=long_k) remainder = int(mod(me%global%nElems, int(nParts, kind=long_k))) ! The first partition starts always with the firstID me%Part_First(1) = firstID ! Up to remainder partitions have share + 1 elements coord = 0 coord(4) = level do iPart=2,remainder+1 coord(1) = coord(1) + int(share) me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) coord(1) = coord(1) + 1 me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The remaining elements get exactly the share elements: do iPart=remainder+2,nParts coord(1) = coord(1) + int(share) - 1 me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) coord(1) = coord(1) + 1 me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The last partition ends always with the lastID me%Part_Last(nParts) = lastID ! Local data: if (myPart < remainder) then me%nElems = int(share+1) me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k) else me%nElems = int(share) me%elemOffset = (int(myPart, kind=long_k) * share) & & + int(remainder, kind=long_k) end if ! All elements have (periodic) boundaries me%Property(1)%nElems = me%nElems me%Property(1)%offset = me%elemOffset allocate(me%Property(1)%ElemID(me%Property(1)%nElems)) ! Please note, that the tem_bc_prop_module will set boundary conditions ! based on this label accordingly! me%global%Property(1)%label = 'internal 1D BC' me%global%Property(1)%bitpos = prp_hasBnd me%global%Property(1)%nElems = me%Property(1)%nElems allocate(me%treeID(me%nElems)) allocate(me%ElemPropertyBits(me%nElems)) ! Only has boundary property: me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd) ! Filling the treeIDs: do iElem=1,me%nElems me%Property(1)%ElemID(iElem) = iElem ! We can only have 2**20 elements per dimension, so in this case where ! we create a line, we won't exceed the integer limit. Thus we safely can ! cast the long integer elemOffset to a normal integer. coord(1) = int(me%elemOffset) + iElem - 1 me%treeID(iElem) = tem_idofcoord(coord, offset=firstID) end do end subroutine generate_treelm_line ! ************************************************************************ ! ! ************************************************************************ ! !> Generate a predefined line with a given number of elements !! !! This serves as an simple grid generation for performance or scaling !! analysis without being obliged to use Seeder. !! You have to specify the generic grid parameters in the lua file instead of !! the mesh folder !! ~~~~~~~~~~~~~~~~~~~~~{.lua} !! mesh = { predefined='line', -- or: 'line_bounded' !! origin = {0.,0.,0.}, !! length = 10., !! elementcount = 6 } !! ~~~~~~~~~~~~~~~~~~~~~ !! !! You have to specify the shape 'line', a bounding box origin, its length !! and the number of elements, which results in different amount elements !! in the grid.\n !! The result of this routine is mainly the treeID list with the additional !! lists for saving the properties. !! The generated line will be a line of elementcount elements along the !! X-Axis with periodicity in all directions. !! !! It is also possible to generate a mesh with boundary conditions in the X !! direction (west and east), by using the predefined 'line_bounded'. !! In this case these two boundary conditions need to be provided in the !! solver configuration. subroutine generate_treelm_elements( me, origin, length, elementcount, & & myPart, nParts, comm, predefined, & & bclabel ) ! -------------------------------------------------------------------- ! !> Mesh to generate type(treelmesh_type), intent(out) :: me !> Corner of the cube real(kind=rk), intent(in) :: origin(3) !> Length of cube real(kind=rk), intent(in) :: length !> Number of elements in the line integer, intent(in) :: elementcount !> Partition of the caller (starts with 0) integer, intent(in) :: myPart !> Number of partitions integer, intent(in) :: nParts !> communicator to be used integer, intent(in) :: comm !> Label describing the internal mesh. character(len=*), intent(in) :: predefined !> Label describing the boundary conditions to set for this mesh. character(len=*), intent(in) :: bclabel ! -------------------------------------------------------------------- ! integer :: level integer(kind=long_k) :: firstID, lastID integer(kind=long_k) :: share integer :: remainder integer :: iPart, iElem integer :: coord(4) integer :: lastcoord integer :: xbound_pad ! -------------------------------------------------------------------- ! me%global%nParts = nParts me%global%myPart = myPart me%global%comm = comm ! Face definitions make the first and last face coincide. ! To allow boundary definitions on both sides, we need to make sure, there ! is at least one additional element. if (predefined == 'line_bounded') then xbound_pad = 1 else xbound_pad = 0 end if ! Find an appropriate level of at least 1. level = max( ceiling(log(real(elementcount+xbound_pad,kind=rk)) & & / log(2.0_rk)), 1 ) me%global%origin = origin me%global%BoundingCubeLength = length*(real(2**level, kind=rk) & & /real(elementcount, kind=rk)) me%global%minLevel = level me%global%maxLevel = level me%global%label = 'Generic_Line' me%global%predefined = predefined write(me%global%comment,'(a15,i7,a16,i2,a1)') & & 'Generated with ', nParts, ' Parts and ', elementcount, & & ' elements.' me%global%dirname = './' ! Boundary property to define periodic boundary me%global%nProperties = 1 if (associated(me%global%property)) deallocate(me%global%property) if (associated(me%property)) deallocate(me%property) allocate(me%global%Property(me%global%nProperties)) allocate(me%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) ! Compute the treeIDs of the mesh: firstID = tem_firstIdAtLevel(level) lastcoord = elementcount - 1 lastID = tem_IdOfCoord(coord = [lastcoord, 0, 0, level], offset = firstID ) ! Total number of elements in this mesh me%global%nElems = elementcount share = me%global%nElems / int(nParts, kind=long_k) remainder = int(mod(me%global%nElems, int(nParts, kind=long_k))) ! The first partition starts always with the firstID me%Part_First(1) = firstID ! Up to remainder partitions have share + 1 elements coord = 0 coord(4) = level do iPart=2,remainder+1 coord(1) = coord(1) + int(share) me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) coord(1) = coord(1) + 1 me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The remaining elements get exactly the share elements: do iPart=remainder+2,nParts coord(1) = coord(1) + int(share) - 1 me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID) coord(1) = coord(1) + 1 me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID) end do ! The last partition ends always with the lastID me%Part_Last(nParts) = lastID ! Local data: if (myPart < remainder) then me%nElems = int(share+1) me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k) else me%nElems = int(share) me%elemOffset = (int(myPart, kind=long_k) * share) & & + int(remainder, kind=long_k) end if ! All elements have (periodic) boundaries me%Property(1)%nElems = me%nElems me%Property(1)%offset = me%elemOffset allocate(me%Property(1)%ElemID(me%Property(1)%nElems)) ! Please note, that the tem_bc_prop_module will set boundary conditions ! based on this label accordingly! me%global%Property(1)%label = trim(bclabel) me%global%Property(1)%bitpos = prp_hasBnd me%global%Property(1)%nElems = me%Property(1)%nElems allocate(me%treeID(me%nElems)) allocate(me%ElemPropertyBits(me%nElems)) ! Only has Boundary Property: me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd) ! Filling the treeIDs: do iElem = 1, me%nElems me%Property(1)%ElemID(iElem) = iElem ! We can only have 2**20 elements per dimension, so in this case where ! we create a line, we won't exceed the integer limit. Thus we safely can ! cast the long integer elemOffset to a normal integer. coord(1) = int(me%elemOffset) + iElem - 1 me%treeID(iElem) = tem_idofcoord(coord, offset=firstID) end do end subroutine generate_treelm_elements ! ************************************************************************ ! ! ************************************************************************ ! !> summary: Generate a predefined single element !> This serves as an simple grid generation for performance or scaling !! analysis without being obliged to use Seeder. !! You have to specify the generic grid parameters in the lua file instead of !! the mesh folder !!```lua !! mesh = { predefined='single', !! origin = {0.,0.,0.}, !! length = 10.} !!``` !! You have to specify the shape 'single', a bounding box origin and its !! length. The generated element will have treeID 1 and is periodic in !! all directions. !! The given length will be the length of the element, and the bounding !! cube will be adapted accordingly, to map this to treeID 1. !! subroutine generate_treelm_single( me, origin, length, myPart, & & nParts, comm ) ! -------------------------------------------------------------------- ! !> Mesh to generate type(treelmesh_type), intent(out) :: me !> Corner of the cube real(kind=rk), intent(in) :: origin(3) !> Length of cube real(kind=rk), intent(in) :: length !> Partition of the caller (starts with 0) integer, intent(in) :: myPart !> Number of partitions integer, intent(in) :: nParts !> communicator to be used integer, intent(in) :: comm ! -------------------------------------------------------------------- ! integer(kind=long_k) :: firstID, lastID integer :: iPart ! -------------------------------------------------------------------- ! me%global%nParts = nParts me%global%myPart = myPart me%global%comm = comm me%global%origin = origin me%global%BoundingCubeLength = length me%global%minLevel = 1 me%global%maxLevel = 1 me%global%label = 'Generic_Single' me%global%predefined = 'single' write(me%global%comment,'(a15,i7,a16,i2,a1)') & & 'Generated with ', nParts, ' Parts on Level 1.' me%global%dirname = './' ! Boundary property to define periodic boundary me%global%nProperties = 1 if (associated(me%global%property)) deallocate(me%global%property) if (associated(me%property)) deallocate(me%property) allocate(me%global%Property(me%global%nProperties)) allocate(me%Property(me%global%nProperties)) allocate(me%Part_First(nParts)) allocate(me%Part_Last(nParts)) ! Compute the treeIDs of the mesh: firstID = 1_long_k lastID = firstID ! Total number of elements in this mesh me%global%nElems = 1_long_k ! The first partition starts always with the firstID me%Part_First(1) = firstID me%Part_Last(1) = lastID ! There is only one element, all partitions but the first one are empty. do iPart=2,nParts me%Part_First(iPart) = lastID me%Part_Last(iPart) = lastID - 1 end do ! Local data: if (myPart == 0) then me%nElems = 1 me%elemOffset = 0_long_k else me%nElems = 0 me%elemOffset = 1_long_k end if ! All elements have (periodic) boundaries me%Property(1)%nElems = me%nElems me%Property(1)%offset = me%elemOffset allocate(me%Property(1)%ElemID(me%Property(1)%nElems)) ! Please note, that the tem_bc_prop_module will set boundary conditions ! based on this label accordingly! me%global%Property(1)%label = 'internal 0D BC' me%global%Property(1)%bitpos = prp_hasBnd me%global%Property(1)%nElems = me%Property(1)%nElems allocate(me%treeID(me%nElems)) allocate(me%ElemPropertyBits(me%nElems)) ! No Properties: me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd) ! There is only one element: me%treeID = firstID end subroutine generate_treelm_single ! ************************************************************************ ! ! ************************************************************************ ! !> This subroutine creates a simple mesh based on a treeID list !! requires the setting of me%global before hand !! Works on the communicator specified in me%global%comm !! subroutine treelmesh_fromList(me, tIDlist, comm, dirname) ! -------------------------------------------------------------------- ! !> Structure to load the mesh to type(treelmesh_type), intent(inout) :: me !> treeID list to generate the global list in the mesh from integer(kind=long_k), intent(in) :: tIDlist(:) !> MPI Communicator to use, defaults to the one in me%global%comm if not !! specified integer, intent(in), optional :: comm !> Directory to store the mesh in. Is taken to be me%global%dirname if not !! specified character(len=*), intent(in), optional :: dirname ! -------------------------------------------------------------------- ! integer :: commLoc integer(kind=long_k) :: offset, nElems integer :: iError integer :: nElemsList ! mpi variables integer :: commsize, rank ! -------------------------------------------------------------------- ! if (present(dirname)) then me%global%dirname = dirname end if if (present(comm)) then commLoc = comm else commLoc = me%global%comm end if ! determine number of treeIDs in list nElemsList = size(tIDlist) ! Get size of communicator and my rank in it call mpi_comm_size(commLoc, commsize, iError) call mpi_comm_rank(commLoc, rank, iError) me%global%nParts = commsize me%global%myPart = rank if (associated(me%property)) deallocate(me%property) allocate(me%Property(me%global%nProperties)) allocate(me%treeID(nElemsList)) allocate(me%ElemPropertyBits(nElemsList)) allocate(me%Part_First(commsize)) allocate(me%Part_Last(commsize)) ! Assign treeID list to the new mesh. Sort first offset = 0 nElems = int(nElemsList, kind=long_k) ! Determine offset for each process (get nElems of ranks before) call MPI_Exscan( nElems, offset, 1, MPI_INTEGER8, MPI_SUM, commLoc, iError ) me%elemOffset = offset me%nElems = nElemsList ! Last process offset + its nElems is total number of elements nElems = offset + nElems ! Distribute among new communicator call MPI_Bcast( nElems, 1, MPI_INTEGER8, commsize-1, commLoc, iError ) ! Store in global number of elements me%global%nElems = nElems ! Copy treeID list to treeID in the mesh me%treeID = tIDlist ! Reset properties. Not needed me%ElemPropertyBits = 0 ! Exchange Part_First and Part_Last call mpi_allgather( me%treeID(1), & & 1, MPI_INTEGER8, & & me%Part_First, 1, MPI_INTEGER8, & & commLoc, iError ) call mpi_allgather( me%treeID(nElemsList), & & 1, MPI_INTEGER8, & & me%Part_Last, 1, MPI_INTEGER8, & & commLoc, iError ) end subroutine treelmesh_fromList ! ************************************************************************ ! ! ************************************************************************ ! !> write mesh information subroutine dump_meshHeader( me ) ! -------------------------------------------------------------------- ! !> fluid tree from mesh type( treelmesh_type ), intent(in) :: me ! -------------------------------------------------------------------- ! integer :: iProp ! -------------------------------------------------------------------- ! write(logunit(1),*) 'Got a mesh with following properties:' write(logUnit(1), "(' Mesh name: ',A)") & & trim(me%global%label) if (trim(me%global%predefined) /= '') then write(logunit(1), "(' This is a predefined mesh: ', A)") & & trim(me%global%predefined) end if write(logunit(1), "(' Total number of elements: ', I0)") me%global%nElems write(logunit(1), "(' Local number of elements: ', I0)") me%nElems write(logUnit(1), "(' Number of properties: ', I0)") & & me%global%nProperties write(logunit(1), "(' Number of partitions: ', I0)") me%global%nParts write(logunit(1), "(' Minimum level: ', I0)") me%global%minlevel write(logunit(1), "(' Maximum level: ', I0)") me%global%maxlevel write(logunit(1), "(' Bounding Box Origin: ', 3F10.4 )") & & me%global%origin write(logunit(1), "(' Bounding Box Length: ', F10.4 )") & & real(me%global%BoundingCubeLength) write(logUnit(1), "(' Coarsest element size: ', F10.4 )") & & real( me%global%BoundingCubeLength / & & real( 2**me%global%minlevel, kind=rk ) ) write(logUnit(1), "(' Finest element size: ', F10.4 )") & & real( me%global%BoundingCubeLength / & & real( 2**me%global%maxlevel, kind=rk ) ) do iProp=1,me%global%nProperties if ( me%global%Property(iProp)%BitPos == prp_hasBnd ) then write(logUnit(1), "(' Number of boundary elements: ', I0 )") & & me%Property(iProp)%nElems end if if ( me%global%Property(iProp)%BitPos == prp_hasQVal ) then write(logUnit(1), "(' Number of qVal elements: ', I0 )") & & me%Property(iProp)%nElems end if if ( me%global%Property(iProp)%BitPos == prp_hasNormal ) then write(logUnit(1), "('Number of elements with wall normals: ', I0 )") & & me%Property(iProp)%nElems end if end do end subroutine dump_meshHeader ! ************************************************************************ ! ! ************************************************************************ ! subroutine tem_load_weights( me, weights, success ) ! -------------------------------------------------------------------- ! type(treelmesh_type), intent(in) :: me real(kind=rk), intent(out) :: weights(me%nElems) logical, intent(out) :: success ! -------------------------------------------------------------------- ! integer(kind=MPI_OFFSET_KIND) :: displacement integer(kind=MPI_OFFSET_KIND) :: filebytes integer :: fh, ftype, iError, typesize integer :: iostatus( MPI_STATUS_SIZE ) integer :: level integer :: iElem logical :: ex character(len=PathLen) :: filename character(len=4) :: EndianSuffix ! -------------------------------------------------------------------- ! EndianSuffix = tem_create_endianSuffix() if (trim(me%weights_file) /= '') then filename = trim(me%weights_file)//EndianSuffix ! check if a corresponding weight file exists if (me%global%myPart == 0) then inquire(file=trim(filename), exist=ex) end if call MPI_Bcast(ex, 1, MPI_LOGICAL, 0, me%global%comm, iError) else ex = .false. end if if (ex) then ! Found a weights file, which is used to read a weight for each ! element. write(logUnit(3),*) 'Loading Weights from file: '//trim(filename) ! Open the binary file for MPI I/O (Write) call MPI_File_open( me%global%comm, trim(filename), & & MPI_MODE_RDONLY, MPI_INFO_NULL, fh, iError ) call check_mpi_error(iError,'file_open in load_weights') ! Create a MPI Subarray as ftype for file view call MPI_Type_contiguous( me%nElems, rk_mpi , ftype, iError ) call check_mpi_error(iError,'type ftype in load_weights') call MPI_Type_commit( ftype, iError ) call check_mpi_error(iError,'commit ftype in load_weights') !get size of etype call MPI_Type_size(rk_mpi, typesize, iError ) call check_mpi_error(iError,'typesize in load_weights') call MPI_File_get_size(fh, filebytes, iError) ! Check whether the weight file has as many values as there are ! elements in the mesh. if (filebytes /= typesize*me%global%nElems) then write(logunit(3),*) 'Mesh file has wrong number of elements!' write(logunit(3),*) 'NO BALANCING!' success = .false. else ! File has correct number of elements matching the mesh. ! calculate displacement displacement = me%elemoffset * typesize * 1_MPI_OFFSET_KIND ! Set the view for each process on the file above call MPI_File_set_view( fh, displacement, rk_mpi, ftype, & & "native", MPI_INFO_NULL, iError ) call check_mpi_error(iError,'set_view in load_weights') ! Read data from the file call MPI_File_read_all( fh, weights, me%nElems, rk_mpi, iostatus, & & iError ) call check_mpi_error(iError,'read_all in load_weights') success = .true. end if !Free the MPI_Datatypes which were created and close the file call MPI_Type_free(ftype, iError) call check_mpi_error(iError,'free ftype in load_weights') call MPI_File_close(fh, iError) call check_mpi_error(iError,'close file in load_weights') else if (maxval(me%levelweight) > 0.0_rk) then write(logunit(3),*) 'Using levelwise weights.' do iElem=1,me%nElems level = tem_levelOf( me%treeID(iElem) ) weights(iElem) = me%levelWeight(level) end do success = .true. else write(logunit(3),*) 'NO BALANCING!' success = .false. end if write(logUnit(3),*) 'Done loading weights.' end subroutine tem_load_weights ! ************************************************************************ ! ! ************************************************************************ ! subroutine tem_dump_weights( me, filename, weights ) ! -------------------------------------------------------------------- ! type(treelmesh_type), intent(in) :: me !> Weights file name character(len=*), intent(in) :: filename real(kind=rk), intent(in) :: weights(me%nElems) ! -------------------------------------------------------------------- ! integer :: nElems integer :: comm integer(kind=long_k) :: offset, globElems integer(kind=long_k) :: filesize integer(kind=MPI_OFFSET_KIND) :: displacement integer :: fh, ftype, iError integer :: iostatus( MPI_STATUS_SIZE ) integer :: typesize ! -------------------------------------------------------------------- ! nElems = me%nElems comm = me%global%comm offset = me%elemOffset globElems = me%global%nElems ! Found a weights file, which is used to read a weight for each ! element. write(logUnit(1),*) 'Dumping Weights to file: '//trim(filename) ! Open the binary file for MPI I/O (Write) call MPI_File_open( comm, trim(filename), & & ior(MPI_MODE_WRONLY,MPI_MODE_CREATE), & & MPI_INFO_NULL, fh, iError ) call check_mpi_error(iError,'file_open in dump_weights') ! Create a MPI Subarray as ftype for file view call MPI_Type_contiguous( nElems, rk_mpi , ftype, iError ) call check_mpi_error(iError,'type ftype in dump_weights') call MPI_Type_commit( ftype, iError ) call check_mpi_error(iError,'commit ftype in dump_weights') !get size of etype call MPI_Type_size(rk_mpi, typesize, iError ) call check_mpi_error(iError,'typesize in dump_weights') ! calculate displacement displacement= offset * typesize * 1_MPI_OFFSET_KIND ! set filesize filesize = globElems * typesize call MPI_File_set_size( fh, filesize, iError ) ! Set the view for each process on the file above call MPI_File_set_view( fh, displacement, rk_mpi, ftype, & & "native", MPI_INFO_NULL, iError ) call check_mpi_error(iError,'set_view in dump_weights') ! Read data from the file call MPI_File_write_all( fh, weights, nElems, rk_mpi, iostatus, iError ) call check_mpi_error(iError,'read_all in dump_weights') !Free the MPI_Datatypes which were created and close the file call MPI_Type_free(ftype, iError) call check_mpi_error(iError,'free ftype in dump_weights') call MPI_File_close(fh, iError) call check_mpi_error(iError,'close file in dump_weights') write(logUnit(1),*) 'Done dumping weights.' end subroutine tem_dump_weights ! ************************************************************************ ! ! ************************************************************************ ! !> Exchange mesh elements with other process !! The following data within treelmesh_type is updated in this routine: !! treeID, ElemPropertyBits, Part_First, Part_Last !! Data pathList is deallocated. subroutine exchange_elements( me, sparta ) ! -------------------------------------------------------------------- ! !> Tree Mesh type !! nElems and offset should be ready before this routine type( treelmesh_type ), intent(inout) :: me !> How many elements to exchange with other processes type( tem_sparta_type ), intent(in) :: sparta ! -------------------------------------------------------------------- ! integer :: iError, iProp ! -------------------------------------------------------------------- ! ! Exchange treeID and propertyBits --------------------------------- call tem_exchange_sparta( me = sparta, & & val = me%treeID, & & nComponents = 1, & & comm = me%global%comm ) call tem_exchange_sparta( me = sparta, & & val = me%ElemPropertyBits, & & nComponents = 1, & & comm = me%global%comm ) ! Exchange treeID and propertyBits --------------------------------- ! Update Part_first and Part_Last ------------------------ call MPI_allgather( me%treeID(1), 1, long_k_mpi, & & me%Part_First, 1, long_k_mpi, & & me%global%comm, iError ) call MPI_allgather( me%treeID(me%nElems), 1, long_k_mpi, & & me%Part_Last, 1, long_k_mpi, & & me%global%comm, iError ) ! Update Part_first and Part_Last ------------------------ do iProp=1,me%global%nProperties call gather_Property( Property = me%Property(iProp), & & Header = me%global%Property(iProp), & & BitField = me%ElemPropertyBits, & & comm = me%global%comm ) end do if (allocated(me%pathList)) deallocate( me%pathList ) end subroutine exchange_elements ! ************************************************************************ ! end module treelmesh_module ! **************************************************************************** ! !> \page motivation Motivation for the TreElm Library !! !! Mesh-based algorithms build one of the most important discretization schemes !! for the numerical modeling of field problems. !! They are also well suited for parallel computations, as the meshes can be !! partitioned and the smaller problems solved concurrently. !! Their basic idea is the discretization of the computational domain by !! covering it fully with smaller, well defined elements. !! For complex geometries, it is in general necessary to use an unstructured !! mesh to discretize the computational domain. !! In unstructured meshes each element has typically a different size and form. !! Therefore, all information of element geometry and neighbor relations have to !! be described explicitly. !! This fully explicit description however, poses a severe hurdle for parallel !! processing, as potentially any part of the mesh might be required in the !! neighborhood of any process. !! Due to the limited memory, providing the complete mesh information locally to !! each partition is not an option, as the ability to process larger meshes is !! often the motivation to switch from one to many computing units in the first !! place. !! Although, it is possible to move this limitation out of the solver itself, !! Tu et al. pointed out the importance to avoid such bottlenecks throughout !! the complete chain of tools in the simulation. !! Thus, not only the main computational loop, but the entire simulation process !! has to be considered, including mesh generation and post-processing. !! The usage of hierarchical subdivisions in a fixed space resulting in and !! octree representation offers an option to overcome the hurdle for parallel !! processing of complex geometries with mesh-based methods. !! !! TreElM provides a basis for massively parallel mesh-based simulations, using !! an octree discretization in combination with a space-filling curve. !! The usage of octree meshes for flow simulations has long been proposed, for !! example by Flaherty et al. in 1997. !! By using octree meshes the main bottleneck of unstructured !! meshes is avoided, while the flexibility of local refinement and resolution !! of arbitrary complex geometries is maintained. !! With the known topology of the octree, the needed information about other !! partitions can be kept at a minimum and the connectivity can be computed !! locally. !! Direct usage of cubic octree elements has the advantage that they can be !! computed efficiently in most schemes, or are even required as for example !! in the lattice Boltzmann method. !! \image html clumps_sm.png !! The figure above shows an example for meshes, that are obtained by the !! octree discretization. !! It shows the filling of a channel with many small cylindric obstacles, !! indicated by transparent objects. !! The complete space covered by the octree is outlined by the blue grid lines !! and shows the refinement towards the channel with obstacles. !! Inside the channel, there is the actual computational domain, shown by the !! green cubical elements. !! !! The basic concept in TreElM is to exploit the known hierarchical !! topology, provided by the octree. !! Another concept that is supported by the octree and beneficial for the !! parallel operation is the strictly elementwise view on the mesh data. !! As each element will be uniquely associated with a certain partition during !! the computation, such a view enables the parallel treatment straight from !! reading or creating the mesh onwards. !! TreElM therefore maintains an elemental view on the mesh data in the data layout !! and the operations on the data. !! !! \section toolchain The Simulation Tool Chain !! !! The tools in the framework range from the mesh generator *Seeder* !! over a finite volume solver for compressible and a !! lattice Boltzmann solver for incompressible flow to the post-processing tool !! *Harvester*, that produces visualization files from simulation data !! attached to the mesh. !! This interaction is sketched in the following Figure. !! \image html apes_schematic.png !! Since neither input nor output (I/O) can be avoided to the largest part, even !! with an integrated workflow, data structures are designed to be suitable for !! efficient parallel reading and writing. !! !! \subsection see_also See Also !! - \ref tem_distributed_octree "The distributed Octree" !! - \ref tem_node_identification "Node identification in the tree" !! !! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID" !! !! \page tem_distributed_octree Distributed Octree !! An [octree](http://en.wikipedia.org/wiki/Octree) !! is obtained by recursive bisection in each dimension of a !! three-dimensional cube. !! As such it is the generalization of a binary tree into three dimensions, the !! two dimensional case is usually referred to as quadtree. !! We use quadtrees and two dimensional meshes for illustrations. !! A three-dimensional picture of the actual octree meshes, obtained by this !! method is provided in the Figure. !! \image html clumps_sm.png !! This illustration shows the subdivision of coarser elements into eight !! smaller elements and the recursive repetition of this refinement where !! necessary. !! The outermost cube is also addressed as the universe, as it contains the !! complete representable space of the octree and all elements have to live !! inside of it. !! The basic idea in TreElM for a scalable parallel deployment of the octree !! discretization is a simple serialized tree representation and sketched in !! the following figure. !! \image html tree_structure_sm.png !! The serialization turns the hierarchical structure of the tree, where each !! node has several links to other nodes into a serial ordered list of one node !! after the other. !! This is achieved by a [space-filling curve](http://en.wikipedia.org/wiki/Space-filling_curve). !! Space filling curves are designed to fill the three dimensional space !! completely with a one-dimensional one. !! They thus map a one-dimensional unit interval onto a higher dimensional unit !! interval. !! This feature can be exploited to define a serialization of the three !! dimensional space. !! They are also defined recursively and complete coverage of the continuous !! space is achieved after infinite iterations. !! !! Finite iterations of the space filling curve definition can be interpreted as !! discrete curves and can be used to describe the ordering for finite intervals !! in the higher dimensional space. !! This recursive definition of the space filling curve ordering nicely fits the !! recursive definition of the spatial subdivisions by the octree. !! We can thus interpret the space filling curve as a rule how the children of a !! node in the octree should be numbered. !! Consecutively numbering all nodes on each refinement level according to the !! space filling curve results in a breadth first numbering of the complete !! tree. !! As these numbers identify the elements in the full tree, we refer to them as !! [treeid]s. !! Their properties are described in the following and the relation !! between mesh and tree is illustrated with the help of a quadtree in the Figure !! above. !! !! \section partitioning Partitioning of the Tree !! With this numbering scheme for all nodes in the full tree it is possible to !! deduce the geometrical position of each node purely by its assigned number. !! Just the [treeid]s of the leaf !! nodes, that are actually part of the !! computational domain are stored with the described ordering to represent the !! mesh. !! Thus, we obtain a sparse octree, that contains only the relevant leaf nodes !! for the computational nodes. !! In the Figure of the 3D mesh, these are highlighted as green cubes. !! None of the other tree nodes, including the parents of the elements in the !! computational domain, are actually stored. !! Due to the ordering, the serialized list of elements is efficiently !! searchable. !! The one dimensional list of elements can be easily partitioned by splitting !! it into individual chunks. !! This design is sketched in the Figure above for a two !! dimensional mesh as shown on the left side of the figure. !! There is an obstacle in the middle of the mesh, covering the four central !! fine elements. !! The mesh is refined towards this obstacle. !! A line is used to indicate the connection between consecutive elements !! according to the ordering with the space filling curve. !! On the right side the according quadtree is shown and the partitioning into !! 6 parts is indicated by regions with different colors. !! Thus, we obtain a simple ordered linear representation of the sparse octree. !! Each actually existing element in the mesh is numbered by position in !! addition to its [treeid] and the resulting serialized tree representation is !! given at the bottom of the Figure in the form of a simple [treeid] list. !! !! \subsection see_also See also !! - \ref treelm_fileformat Treelm File Format !! - \ref tem_balance_module "Load Balancing" !! !! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID" !! !! \page treelm_fileformat Treelm File Format !! As already introduced in the previous sections, the !! \ref tem_distributed_octree "octree mesh" is represented by the sparse !! storage of the leaf nodes in the computational domain. !! The storage of the mesh on disk is basically the sequence of [treeid]s !! building the computational domain. !! However, this is in general insufficient to fully describe the mesh and an !! additional elemental information is introduced to attach further data to the !! elements. !! TreElM uses 64-bit signed integers to encode the [treeid]s for !! greatest portability. !! In addition, a second 64-bit signed integer is used to build a bit-mask for !! properties like boundary conditions, that might be attached to each element. !! This data is written in native binary format to disk, resulting in a file !! with 16 bytes per element. !! Due to this simple format, the data might be accessed on an elemental basis !! with Fortran's direct IO, it can be easily converted between big and little !! endian representations and allows fast reading and writing. !! !! This binary data is accompanied by a header file with descriptive information !! about the mesh, like the total number of elements, the physical extent and !! origin of the universal cube and descriptions for the attached properties. !! The data in this header is provided in the form of a [Lua](http://www.lua.org/) !! script, to allow for a flexible handling of additional data and future evolution !! of the format. !! With the help of the \ref aotus_module "Aotus library" data can be easily !! retrieved from and written to such Lua scripts. !! !! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID" !! !! !! \page implementation Implementation in the APES Framework !! In distributed computations it is important to avoid any kind of bottleneck, !! as the local computing resources in each computing unit are !! usually quite limited in comparison to the total computing power. !! For large scale computations that aim to use the full system for their !! simulation, it therefore is essential that the complete processing chain is !! distributable. !! Otherwise any step from the mesh generation over the simulation to the !! post-processing might prohibit a successful simulation. !! The *TreElM* library builds the basis for a completely !! integrated framework that allows the distributed processing of all necessary !! simulation steps. !! To enable a flexible and convenient user configuration, the \ref aotus_module !! "Aotus library" !! is not only used in *TreElM* header files, but also exposed to the !! calling applications. !! The layout of the overall framework with the central *TreElM* library is !! shown by the schematic figure in the \ref motivation "Motivation". !! !! The solvers read the mesh data, using *TreElM* functions and process the !! mesh into a solver specific data structure with attached simulation data. !! Results are written in the form of restart files, that follow the elemental !! design of the *TreElM* mesh by writing the elemental solution following the !! same ordering to disk. !! With the help of tracking objects it is possible to write subsets of the mesh !! with attached simulation data in arbitrary intervals to disk. !! The *TreElM* library not only provides the means to identify the halos to !! be exchanged between partitions but also provides various communication !! pattern that can be used for the actual exchange. !! Due to this encapsulation of the communication, the communication layout can !! be easily exchanged. !! It is even possible to replace the complete parallel paradigm in the !! deployment with only few changes to the code. !! !! Aside from these basic functionalities the library provides several auxiliary !! routines, that provide for example logging and debugging facilities. !! Finally the post processing tool *Harvester* is used to visualize the !! results of the simulation. !! Since *Harvester* is a stand alone tool, it can be deployed on specialized !! machines, that are more suited for visualization. !! Additionally, the separated design ensures the independence of the solvers !! from any third party libraries, that are typically required for visualization !! outputs. !! !! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID" !! !! !! \page tem_requirements Requirements !! !! Below are the most important requirements for getting, compiling and running !! solvers developed with the TreElm framework. !! !! Compilation !! ----------- !! - [Mercurial](mercurial.selenic.com): Getting sources !! - [Python](www.python.org): the build tool waf relies on python and most of !! the scripts inside the package are written with Python !! - Compiler: the vast source files are written in Fortran 95/2003. !! We recommend [GNU Gfortran 4.6+](http://gcc.gnu.org/wiki/GFortran) !! and [Intel 12.1+](http://software.intel.com/en-us/articles/ !! intel-compilers/) !! - MPI 2.0 enabled environment !! We recommend [OpenMPI 1.5+](http://www.open-mpi.org/) !! !! Post-Processing and Visualization !! --------------------------------- !! - [Gnuplot](http://www.gnuplot.info/) !! - [Paraview](http://www.paraview.org/) !! - [Tecplot](http://www.tecplot.com)