Mesh implementation
The data structures:
The mesh descriptor is defined in
grid/mesh.F90
:
type, extends(basis_set_abst_t) :: mesh_t
! Components are public by default
class(box_t), pointer :: box !< simulation box; polymorphic pointer
class(coordinate_system_t), pointer :: coord_system !< the underlying coordinate system (affine, cubic, curvilinear, etc.)
type(index_t) :: idx !< the indexing scheme.
logical :: use_curvilinear
real(real64), allocatable :: spacing(:) !< the (constant) spacing between the points
! When running serially, the local number of points is
! equal to the global number of points.
! Otherwise, the next two are different on each node.
integer :: np !< Local number of points in mesh
integer :: np_part !< Local points plus ghost points plus boundary points.
integer(int64) :: np_global !< Global number of points in mesh.
integer(int64) :: np_part_global !< Global number of inner points and boundary points.
logical :: parallel_in_domains !< will I run parallel in domains?
type(mpi_grp_t) :: mpi_grp !< the mpi group describing parallelization in domains
type(par_vec_t) :: pv !< describes parallel vectors defined on the mesh.
type(partition_t) :: partition !< describes how the inner points are assigned to the domains
real(real64), allocatable :: x(:,:) !< The (local) \b points
real(real64) :: volume_element !< The global volume element.
real(real64), allocatable :: vol_pp(:) !< Element of volume for curvilinear coordinates.
logical :: masked_periodic_boundaries !< flag whether to use the mask defined below
character(len=256) :: periodic_boundary_mask !< mask for which the periodic boundary conditions
!! are replaced by zero boundary conditions
contains
procedure :: end => mesh_end !< @copydoc mesh_oct_m::mesh_end
procedure :: init => mesh_init !< @copydoc mesh_oct_m::mesh_init
procedure :: write_info => mesh_write_info !< @copydoc mesh_oct_m::mesh_write_info
procedure :: dmesh_allreduce_0, zmesh_allreduce_0, imesh_allreduce_0, lmesh_allreduce_0
procedure :: dmesh_allreduce_1, zmesh_allreduce_1, imesh_allreduce_1, lmesh_allreduce_1
procedure :: dmesh_allreduce_2, zmesh_allreduce_2, imesh_allreduce_2, lmesh_allreduce_2
procedure :: dmesh_allreduce_3, zmesh_allreduce_3, imesh_allreduce_3, lmesh_allreduce_3
procedure :: dmesh_allreduce_4, zmesh_allreduce_4, imesh_allreduce_4, lmesh_allreduce_4
procedure :: dmesh_allreduce_5, zmesh_allreduce_5, imesh_allreduce_5, lmesh_allreduce_5
generic :: allreduce => dmesh_allreduce_0, zmesh_allreduce_0, imesh_allreduce_0, lmesh_allreduce_0
generic :: allreduce => dmesh_allreduce_1, zmesh_allreduce_1, imesh_allreduce_1, lmesh_allreduce_1
generic :: allreduce => dmesh_allreduce_2, zmesh_allreduce_2, imesh_allreduce_2, lmesh_allreduce_2
generic :: allreduce => dmesh_allreduce_3, zmesh_allreduce_3, imesh_allreduce_3, lmesh_allreduce_3
generic :: allreduce => dmesh_allreduce_4, zmesh_allreduce_4, imesh_allreduce_4, lmesh_allreduce_4
generic :: allreduce => dmesh_allreduce_5, zmesh_allreduce_5, imesh_allreduce_5, lmesh_allreduce_5
end type mesh_t
which uses the structure index_t
, containing the range and the mapping arrays:
type index_t
! Components are public by default
integer :: dim !< the dimension
integer, allocatable :: nr(:,:) !< dimensions of the box where the points are contained
integer, allocatable :: ll(:) !< literally nr(2,:) - nr(1,:) + 1 - 2*enlarge(:)
integer, allocatable :: enlarge(:) !< number of points to add for boundary conditions
integer(int64) :: checksum
integer(int64), pointer :: grid_to_spatial_global(:) !< map: global grid index -> spatial index
integer(int64), pointer :: spatial_to_grid_global(:) !< inverse map: spatial index -> global grid index
integer :: type !< index type: IDX_CUBIC or IDX_HILBERT
integer :: bits !< bits per dimension for Hilbert index
integer, allocatable :: offset(:) !< offset for getting the indices from the spatial index
integer, allocatable :: stride(:)
type(MPI_Win) :: window_grid_to_spatial !< handle to shared memory window for map
type(MPI_Win) :: window_spatial_to_grid !< handle to shared memory window for inverse map
end type index_t
About lxyz and lxyz_inv maps:
The direct map:
lxyz(d, i) = R(d)i
where d denotes the dimension (i.e. x, y or z) and i is the index of the point.
The inverse map:
lxyz_inv( Rxi , Ryi , Rzi ) = i
The points defined by lxyz define a rectangular box of d dimensions. The real mesh vectors are related to Ri by multiplication with spacing and possibly distortion for curvilinear coordinates.
Note that this index array is the same in all parallel domains, and relates the integer coordinates to the global index of a given point.
type mesh_cube_map_t
! Components are public by default
integer :: nmap !< The number of maps
integer, allocatable :: map(:, :)
type(accel_mem_t) :: map_buffer
end type mesh_cube_map_t
Note on packed states:
Note on curvilinear meshes:
The mesh::x(:,:)
array always contains a regular mesh, which gets ‘distorded’ to a curvilinear mesh by additional function calls.