mesh_cart_m.f90 Source File

module mesh_cart_m
    !! Module implementing the 2d mesh
    use netcdf
    use precision_m, only: FP, FP_LARGEST
    use comm_handling_m, only: is_master
    use equilibrium_m,    only: equilibrium_t
    use slab_equilibrium_m, only: slab_t
    use descriptors_m, only: DISTRICT_CORE, DISTRICT_CLOSED, DISTRICT_SOL, &
                             DISTRICT_WALL, DISTRICT_PRIVFLUX, DISTRICT_DOME, &
                             DISTRICT_OUT
    use error_handling_m, only: handle_error, error_info_t
    use status_codes_m, only: PARALLAX_ERR_MESH, PARALLAX_WRN_MESH
    use status_codes_m, only: PARALLAX_ERR_MESH
    use, intrinsic :: iso_c_binding, only: c_int32_t, c_bool, c_ptr, c_loc, &
                                           c_null_ptr
    implicit none

    integer, parameter, private :: PINFO_INNER      = 1
    !! Descriptor (internal) for inner grid point
    integer, parameter, private :: PINFO_BOUNDARY   = 2
    !! Descriptor (internal) for boundary point
    integer, parameter, private :: PINFO_GHOST      = 3
    !! Descriptor (internal) for ghost point

    type, public, bind(C) :: mesh_cart_data_t
        !! Data type used to expose `mesh_cart_t` data to C/C++
        !! Please see `mesh_cart_t` for descriptions of individual
        !! properties/members.
        !! Values are copied directly.
        !! Arrays are exposed with pointers.
        logical(c_bool)    :: extend_beyond_wall
        integer(c_int32_t) :: lvl
        integer(c_int32_t) :: lvst
        real(FP)           :: spacing_f
        real(FP)           :: spacing_c
        integer(c_int32_t) :: n_points

        real(FP)           :: phi
        real(FP)           :: xmin
        real(FP)           :: ymin
        integer(c_int32_t) :: nx_f
        integer(c_int32_t) :: ny_f
        type(c_ptr)        :: cart_i
        type(c_ptr)        :: cart_j

        integer(c_int32_t) :: size_neighbor
        type(c_ptr)        :: index_neighbor
        integer(c_int32_t) :: size_ghost_layer

        integer(c_int32_t) :: n_points_inner
        type(c_ptr)        :: inner_indices
        integer(c_int32_t) :: n_points_boundary
        type(c_ptr)        :: boundary_indices
        integer(c_int32_t) :: n_points_ghost
        type(c_ptr)        :: ghost_indices
        type(c_ptr)        :: pinfo

        type(c_ptr)        :: district

        integer(c_int32_t) :: n_points_red
        integer(c_int32_t) :: n_points_black
        type(c_ptr)        :: redblack_indices
    end type mesh_cart_data_t

    type, public :: mesh_cart_t
        !! Datatype for a mesh_2d
        logical, private :: extend_beyond_wall
        !! If true, mesh is extended beyond the wall
        integer, private :: lvl
        !! Level of coarseness of grid, required for multigrid solver lvl=1 is finest level.
        !! lvl=2,3,4.. correspond to a increased coarseness by a factor of 2,4,8,... respectively.
        integer, private :: lvst
        !! Discrete grid distance
        real(FP), private :: spacing_f
        !! Grid spacing on the finest level
        real(FP), public :: spacing_c
        !! Grid spacing on the actual level (coarse)
        !! NOTE: This variable is set public only temporarily and must not be
        !!  changed from outside of this module.
        !!  In normal circumstances use the routine get_spacing_c.
        !!  (It needs to be accessed directly in helmholtz_operators,
        !!  to circumvent the known issue of OMP scaling with intel compiler)
        integer, public :: n_points
        !! Number of total mesh points
        !! NOTE: This variable is set public only temporarily and must not be
        !!  changed from outside of this module.
        !!  In normal circumstances use the routine get_n_points.
        !!  (It needs to be accessed directly in helmholtz_operators,
        !!  to circumvent the known issue of OMP scaling with intel compiler)
        real(FP), private :: phi
        !! Toroidal angle of mesh, normalized to range [0, 2pi)
        real(FP), private :: xmin
        !! Reference point to compute x-coordinate
        real(FP), private :: ymin
        !! Reference point to compute y-coordinate
        integer, private :: nx_f
        !! x-dimension of rectangular Cartesian unmasked grid (finest level)
        integer, private :: ny_f
        !! y-dimension of rectangular Cartesian unmasked grid (finest level)
        integer, private, allocatable, dimension(:) :: cart_i
        !! Indices in Cartesian logical space in x direction (x_i, y_j)
        integer, private, allocatable, dimension(:) :: cart_j
        !! Indices in Cartesian logical space in y direction (x_i, y_j)
        integer, private, allocatable, dimension(:,:) :: cart_shaped_l
        !! A 2d array where the (i, j)^th element gives the index l
        !! such that cart_i(l) == i and cart_j(l) == j
        integer, private :: min_i, max_i, min_j, max_j
        !! Limits of the cart_i and cart_j arrays

        integer, private :: size_neighbor
        !! Number of neighbor points in each direction (left, right, bottom, top)
        integer, public, allocatable, dimension(:,:,:) :: index_neighbor
        !! Index of neighbor points
        integer, private :: size_ghost_layer
        !! Depth of ghost layer

        integer, private :: n_points_inner
        !! Number of inner points
        integer, public, allocatable, dimension(:) :: inner_indices
        !! Indices of inner points
        integer, private :: n_points_boundary
        !! Number of boundary points
        integer, public, allocatable, dimension(:) :: boundary_indices
        !! Indices of boundary points
        integer, private :: n_points_ghost
        !! Number of ghost points
        integer, public, allocatable, dimension(:) :: ghost_indices
        !! Indices of ghost points
        integer, private, allocatable, dimension(:) :: pinfo
        !! Stores information for mesh point

        integer, private, allocatable, dimension(:) :: district
        !! Stores information on logical location of grid point

        logical, private :: yperiodic
        !! Only relevant for SLAB equilibrium, indicates periodicity in y-direction

        integer, private :: n_points_red
        !! Number of red points (chequerboard)
        integer, private :: n_points_black
        !! Number of black points (chequerboard)
        integer, public, allocatable, dimension(:) :: redblack_indices
        !! Index list of redblack points, eventually needed for red-black gauss seidel smoother
        !! First nred entries are indices of red points
        !! Following nblack entries are indices of black points
        !! Final nboundary+nghost entries are indices of boundary and ghost points

    contains
        ! Getter routines
        procedure, public :: get_phi
        procedure, public :: set_phi
        procedure, public :: get_lvl
        procedure, public :: get_lvst
        procedure, public :: get_n_points
        procedure, public :: get_spacing_f
        procedure, public :: get_spacing_c
        procedure, public :: get_size_neighbor
        procedure, public :: get_cart_i
        procedure, public :: get_cart_j
        procedure, public :: get_x
        procedure, public :: get_y
        procedure, public :: get_n_points_inner
        procedure, public :: get_n_points_boundary
        procedure, public :: get_n_points_ghost
        procedure, public :: get_xmin
        procedure, public :: get_ymin
        procedure, public :: get_n_points_red
        procedure, public :: get_n_points_black

        procedure, public :: is_inner_point
        procedure, public :: is_boundary_point
        procedure, public :: is_ghost_point
        procedure, public :: get_district

        procedure, public :: cart_to_mesh_index
        procedure, public :: build_shaped_cart_arrays
        procedure, public :: deallocate_shaped_cart_arrays
        procedure, public :: get_index_neighbor

        procedure, public :: determine_interpolation_stencil
        procedure, public :: get_surrounding_indices

        procedure, public :: normal_to_boundary

        ! I/O rotines
        procedure, public :: display => display_mesh_cart
        procedure, public :: write_netcdf => write_netcdf_mesh_cart
        procedure, public :: read_netcdf => read_netcdf_mesh_cart

        ! Build routines
        procedure, public :: init => init_mesh_cart
        procedure, private :: build_cart
        procedure, private :: build_cart_slab
        procedure, private :: build_connectivity
        procedure, private :: build_boundary
        procedure, private :: build_ghost_layer
        procedure, private :: build_pinfo
        procedure, private :: build_district
        procedure, private :: build_patch
        procedure, private :: build_redblack
        procedure, private :: point_inside
        procedure, public :: reorder

        ! Routines for MPI comminication of mesh
        procedure, public :: exchange_mesh_mpi

        ! C interoperability routines
        procedure, public :: expose_data
    end type

    interface

        pure real(FP) module function get_phi(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        module subroutine set_phi(self, new_phi)
            !! Sets phi value of mesh to a new value
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            real(FP), intent(in) :: new_phi
            !! New phi value
        end subroutine

        pure integer module function get_lvl(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_lvst(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_n_points(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure real(FP) module function get_spacing_f(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure real(FP) module function get_spacing_c(self)
            !! Returns mesh distance of level, spacing_c = spacing_f*lvst
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_size_neighbor(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_cart_i(self, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure integer module function get_cart_j(self, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure real(FP) module function get_x(self, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure real(FP) module function get_y(self, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure integer module function get_n_points_inner(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_n_points_boundary(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_n_points_ghost(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure real(FP) module function get_xmin(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure real(FP) module function get_ymin(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_n_points_red(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_n_points_black(self)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end function

        pure integer module function get_index_neighbor(self, ishift, jshift, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ishift
            !! Shift in horizontal direction w.r.t. point ind
            integer, intent(in) :: jshift
            !! Shift in vertical direction w.r.t. point ind
            integer, intent(in) :: ind
            !! Mesh index
        end function

        integer module function cart_to_mesh_index(self, i, j)
            !! Given Cartesian indizes (i,j) returns mesh index
            !! if no point could be found 0 is returned
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: i
            !! Horizontal Cartesian index
            integer, intent(in) :: j
            !! Vertical Cartesian index
        end function

        module subroutine build_shaped_cart_arrays(self)
            !! Builds a 2D array where shaped_cart_l(i,j) = l gives the index
            !! where cart_i(l) == i and cart_j(l) == j, or zero if no point
            !! is found. Used to accelerate cart_to_mesh_index
            class(mesh_cart_t), intent(inout) :: self
        end subroutine

        module subroutine deallocate_shaped_cart_arrays(self)
            !! Deallocates the shaped_cart_l array to reduce the memory
            !! requirement for the grid
            class(mesh_cart_t), intent(inout) :: self
        end subroutine

        pure logical module function is_inner_point(self, ind)
            !! Indicates if point is inner grid point
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure logical module function is_boundary_point(self, ind)
            !! Indicates if point is boundary point
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure logical module function is_ghost_point(self, ind)
            !! Indicates if point is ghost point
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        pure integer module function get_district(self, ind)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index
        end function

        module subroutine determine_interpolation_stencil(self, xp, yp, &
            intorder, ns_max_nearest, ns_stencil, intorder_actual, &
            ind_interpol_stencil)
            !! Finds ns_stencil x ns_stencil indices of interpolation stencil
            !! around (xp, yp), if possible of order intorder
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            real(FP), intent(in) :: xp
            !! x-coordinate of point
            real(FP), intent(in) :: yp
            !! y-coordinate of point
            integer, intent(in) :: intorder
            !! Requested stencil order
            integer, intent(in) :: ns_max_nearest
            !! Size of area that is searched for nearest point
            integer, intent(out) :: ns_stencil
            !! Size of determined stencil
            integer, intent(out) :: intorder_actual
            !! Order of determined stencil
            integer, allocatable, dimension(:,:), intent(out) :: ind_interpol_stencil
            !! Indices of determined stencil
        end subroutine

        module subroutine get_surrounding_indices(self, xp, yp, ns, ind_sur, ns_complete, ind_nearest)
            !! Finds ns x ns indices of points surrounding (xp, yp)
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            real(FP), intent(in) :: xp
            !! x-coordinate of point
            real(FP), intent(in) :: yp
            !! y-coordinate of point
            integer, intent(in) :: ns
            !! Dimension of surrounding local grid
            integer, intent(out), optional, dimension(ns, ns) :: ind_sur
            !! Indices of surrounding grid points
            integer, intent(out), optional :: ns_complete
            !! Dimension of surrounding local grid, that is complete (w/o missing points)
            integer, intent(out), optional :: ind_nearest
            !! Index of nearest point on surrounding grid
        end subroutine

        module subroutine normal_to_boundary(self, ind, nbx, nby)
            !! Returns the discrete normal vector to the boundary
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: ind
            !! Mesh index, must be index of a boundary point
            integer, intent(out) :: nbx
            !! x-component of discrete normal vector
            integer, intent(out) :: nby
            !! y-component of discrete normal vector
        end subroutine

        module subroutine display_mesh_cart(self)
            !! Displays information for mesh
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
        end subroutine

        module subroutine write_netcdf_mesh_cart(self, fgid)
            !! Writes the mesh to netcdf
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in) :: fgid
            !! File or group id number of existing Netcdf4 file
        end subroutine

        module subroutine read_netcdf_mesh_cart(self, fgid)
            !! Reads mesh from netcdf
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in) :: fgid
            !! File or group id number of existing Netcdf4 file
        end subroutine

        module subroutine init_mesh_cart(self, equi, phi, lvl, spacing_f, &
                                         size_neighbor, size_ghost_layer, &
                                         extend_beyond_wall, dbgout)
            !! Builds complete mesh
            class(mesh_cart_t), intent(out) :: self
            !! Instance of the type
            class(equilibrium_t), intent(inout) :: equi
            !! Equilibrium
            real(FP), intent(in) :: phi
            !! Toroidal angle
            integer, intent(in) :: lvl
            !! Coarseness level
            real(FP), intent(inout) :: spacing_f
            !! Grid spacing on finest level
            integer, intent(in) :: size_neighbor
            !! Number of neighbor points determined in each direction
            integer, intent(in) :: size_ghost_layer
            !! Depth of ghost point layer
            logical, intent(in) ::extend_beyond_wall
            !! If true, mesh is extended beyond the wall
            integer, intent(in), optional ::dbgout
            !! Debug output level
        end subroutine

        module subroutine build_cart(self, equi, lvl, spacing_f)
            !! Builds Cartesian grid
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            class(equilibrium_t), intent(inout) :: equi
            !! Equilibrium
            integer, intent(in) :: lvl
            !! Coarseness level
            real(FP), intent(in) :: spacing_f
            !! Grid spacing on finest level
        end subroutine

        module subroutine build_connectivity(self, size_neighbor)
            !! Establishes connectivity among grid points
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in) :: size_neighbor
            !! Number of neighbor points in each direction (left, right, bottom, top)
        end subroutine

        module subroutine build_boundary(self)
            !! Builds discrete boundary, i.e. finds point that are located on boundary
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
        end subroutine

        module subroutine build_ghost_layer(self, size_ghost_layer)
            !! Builds ghost layer surrounding inner grid
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in) :: size_ghost_layer
            !! Depth of ghost layer
        end subroutine

        module subroutine build_pinfo(self)
            !! Build information for individual mesh points
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
        end subroutine

        module subroutine build_district(self, equi)
            !! Build district (logical location) for individual mesh points
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            class(equilibrium_t), intent(in) :: equi
            !! Equilibrium
        end subroutine

        module subroutine build_patch(self)
            !! Applies patch to mesh
            !! boundary points without connection to inner points are moved the ghost points
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
        end subroutine

        module subroutine build_redblack(self)
            !! Builds index list identifying red and black points (chequerboard) of mesh
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
        end subroutine

        module subroutine build_cart_slab(self, equi, lvl, npwr)
            !! Builds Cartesian grid for slab geometry
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            class(equilibrium_t), intent(inout) :: equi
            !! Equilibrium
            integer, intent(in) :: lvl
            !! Coarseness level
            integer, intent(in) :: npwr
            !! Power exponent for number of grid points in x direction nxf = 2^npwr+1
        end subroutine

        module subroutine reorder(self, grid_reorder_block, dbgout)
            !! Reorders mesh indices for cache optimisation
            class(mesh_cart_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in)  :: grid_reorder_block
            !! Size of basic unshuffled grid square
            integer, intent(in), optional ::dbgout
            !! Debug output level
        end subroutine

        logical module function point_inside(self, district)
            !! Definition which district (point location in logical sense)
            !! is inside the mesh
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: district
            !! District descriptor
        end function

        module subroutine exchange_mesh_mpi(self, comm, step, mesh_received)
            !! MPI communication of mesh
            !! Sends mesh to rank-step and receives from rank+step
            class(mesh_cart_t), intent(in) :: self
            !! Instance of the type
            integer, intent(in) :: comm
            !! MPI communicator
            integer, intent(in) :: step
            !! Step size of MPI communication
            type(mesh_cart_t), intent(out) :: mesh_received
            !! Mesh received from rank+step
        end subroutine

        module subroutine expose_data(self, data_object)
            !! Exposes mesh data through mesh_cart_data_t object
            class(mesh_cart_t), intent(in), target :: self
            !! Instance of the type
            type(mesh_cart_data_t), intent(out)  :: data_object
            !! Destination mesh cart data object.
            !! **Warning**: the purpose of this method is to allow
            !! access to the data of the instance.
        end subroutine

    end interface

end module