module multigrid_m !! Multigrid for solving the Helmholtz problem use netcdf use precision_m, only : FP use screen_io_m, only : get_stderr use comm_handling_m, only : is_master use descriptors_m, only : BND_TYPE_NEUMANN use csrmat_m, only : csrmat_t, csrmat_data_t, csr_times_vec use mesh_cart_m, only : mesh_cart_t, mesh_cart_data_t use equilibrium_m, only : equilibrium_t use boundaries_perp_m, only : set_perpbnds use, intrinsic :: iso_c_binding, only : c_int32_t, c_ptr, c_loc implicit none type, public, bind(C) :: multigrid_data_t !! Data type used to expose `multigrid_t` data to C/C++ !! Please see `multigrid_t` for descriptions of individual !! properties/members. !! `nlvls` value is copied directly. !! Arrays are exposed with pointers. integer(c_int32_t) :: nlvls type(c_ptr) :: mesh_finest_ptr type(c_ptr) :: mesh_coarse_array type(c_ptr) :: mrcsr type(c_ptr) :: mrcsr_inner type(c_ptr) :: mpcsr type(c_ptr) :: np_lvl type(c_ptr) :: np_boundary_lvl type(c_ptr) :: np_inner_lvl type(c_ptr) :: first_idx_lvl type(c_ptr) :: first_idx_boundary_lvl type(c_ptr) :: first_idx_inner_lvl end type multigrid_data_t type, public :: multigrid_intermediate_data_t !! Data type needed to expose `multigrid_t` data to C/C++. !! Objects of this type shall hold arrays of types !! `mesh_cart_data_t` and `csrmat_data_t`. !! These are strictly required for the exposure (we cannot use !! the arrays of type `mesh_cart_t` and `csrmat_t` that exist within !! `multigrid_t`). !! In princple this data type could be merged into `multigrid_t`, !! but keeping it separate should simplify maintenance. type(mesh_cart_data_t), public :: mesh_finest_ptr_data type(mesh_cart_data_t), public, allocatable, dimension(:) :: & mesh_coarse_array_data type(csrmat_data_t), public, allocatable, dimension(:) :: mrcsr_data type(csrmat_data_t), public, allocatable, dimension(:) :: mrcsr_inner_data type(csrmat_data_t), public, allocatable, dimension(:) :: mpcsr_data contains final :: destructor_multigrid_intermediate_data end type multigrid_intermediate_data_t type, public :: multigrid_t !! Object which contains the multigrid meshes for solving Helmholtz !! problems. This includes nlvl number of meshes of increasing !! coarseness, restriction and prolongation matrices. integer, private :: nlvls !! Number of multigrid levels type(mesh_cart_t), private, pointer :: mesh_finest_ptr => null() !! Pointer to the finest mesh (lvl=1) type(mesh_cart_t), private, dimension(:), allocatable :: & mesh_coarse_array !! The remaining coarser meshes (indexed lvl=2:nlvls) type(csrmat_t), private, allocatable, dimension(:) :: mrcsr !! Restriction matrices (trivial) type(csrmat_t), private, allocatable, dimension(:) :: mrcsr_inner !! Restriction matrices for inner points type(csrmat_t), private, allocatable, dimension(:) :: mpcsr !! Prolongation matrices integer, private, allocatable, dimension(:) :: np_lvl !! Number of total mesh points at each level integer, private, allocatable, dimension(:) :: np_boundary_lvl !! Number of boundary mesh points at each level integer, private, allocatable, dimension(:) :: np_inner_lvl !! Number of inner mesh points at each level integer, private, allocatable, dimension(:) :: first_idx_lvl !! For mesh points stored consecutively from different levels. The size !! of first_idx_lvl is nlvls + 1 !! first_idx_lvl(lvl) gives first index of mesh points from lvl !! first_idx_lvl(lvl + 1) - 1 gives last index of mesh points from lvl integer, private, allocatable, dimension(:) :: first_idx_boundary_lvl !! For boundary points stored consecutively from different levels integer, private, allocatable, dimension(:) :: first_idx_inner_lvl !! For inner points stored consecutively from different levels contains procedure, public :: init => initialize_multigrid procedure, public :: write_netcdf => multigrid_write_netcdf procedure, public :: read_netcdf => multigrid_read_netcdf procedure, public :: set_phi_coarse procedure, public :: copy procedure, public :: get_nlvls procedure, public :: get_np procedure, public :: get_np_boundary procedure, public :: get_np_inner procedure, public :: get_mesh_pointer procedure, public :: get_first_idx_lvl_pointer procedure, public :: get_first_idx procedure, public :: get_first_idx_boundary_lvl_pointer procedure, public :: get_first_idx_boundary procedure, public :: get_first_idx_inner_lvl_pointer procedure, public :: get_first_idx_inner procedure, public :: get_np_total procedure, public :: get_np_inner_total procedure, public :: restrict procedure, public :: restrict_inner procedure, public :: prolong procedure, public :: extrapolate_boundary_points procedure, public :: expose_data procedure, private :: initialize_np_and_idx_arrays procedure, private :: create_mrcsr procedure, private :: create_mrcsr_inner procedure, private :: create_mpcsr procedure, public :: get_mrcsr_pointer procedure, public :: get_mpcsr_pointer final :: destructor_multigrid end type interface module subroutine initialize_multigrid(self, equi, phi, nlvls, & spacing_f, size_neighbor, & size_ghost_layer, mesh_finest, & reorder_size, extend_beyond_wall, & dbgout) !! Sets up multigrids, i.e. !! - creates meshes on all levels !! - creates prolongation and restriction operators class(multigrid_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium (not changed) real(FP), intent(in) :: phi !! Toroidal angle integer, intent(in) :: nlvls !! Number of multigrid levels real(FP), intent(inout) :: spacing_f !! Grid spacing on finest level (only adapted possibly for slab !! geometry) integer, intent(in) :: size_neighbor !! Number of neighbor points determined in each direction (for !! finest level, on coarser levels always 1) integer, intent(in) :: size_ghost_layer !! Depth of ghost point layer (for finest level, on coarser levels !! always 0) type(mesh_cart_t), target, intent(out) :: mesh_finest !! Finest mesh (lvl=1) that was created integer, intent(in), dimension(nlvls), optional :: reorder_size !! Grid reordering for each level logical, intent(in) :: extend_beyond_wall !! If true, meshes are extend beyond the wall integer, intent(in), optional :: dbgout !! Debug output level, default = 0, none end subroutine module subroutine initialize_np_and_idx_arrays(self) !! Creates arrays containing information on number of points !! (total, boundary, and inner) for all mesh levels. The np_lvl !! arrays can be used for defining explicit-shape arrays; the !! first_idx_lvl arrays can be used for consecutively-stored mesh !! data. class(multigrid_t), target, intent(inout) :: self !! Instance of the type end subroutine module subroutine multigrid_write_netcdf(self, fgid) !! Writes information of multigrid to netcdf file class(multigrid_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 multigrid_read_netcdf(self, fgid, mesh_finest) !! Reads information of multigrid from netcdf file class(multigrid_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: fgid !! File or group id number of existing Netcdf4 file type(mesh_cart_t), target, intent(out) :: mesh_finest !! Finest mesh (lvl=1) that was read from file end subroutine module subroutine restrict(self, lvl_f, uf, lvl_c, uc) !! Restricts quantity from fine grid to coarse grid class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl_f !! Level of fine grid real(FP), dimension(:), intent(in) :: uf !! Quantity on fine grid integer, intent(in) :: lvl_c !! Level of coarse grid real(FP), dimension(:), intent(out) :: uc !! Quantity on coarse grid end subroutine module subroutine restrict_inner(self, lvl_f, uf, lvl_c, uc) !! Restricts quantity defined only on inner grid points from fine !! grid to coarse grid class(multigrid_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: lvl_f !! Level of fine grid real(FP), dimension(:), intent(in) :: uf !! Quantity on fine grid integer, intent(in) :: lvl_c !! Level of coarse grid real(FP), dimension(:), intent(out) :: uc !! Quantity on coarse grid end subroutine module subroutine prolong(self, lvl_c, uc, lvl_f, uf) !! Prolongs quantity from coarse grid to fine grid class(multigrid_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: lvl_c !! Level of coarse grid real(FP), dimension(:), intent(in) :: uc !! Quantity on coarse grid ! Should have dimension self%get_np(lvl_c) integer, intent(in) :: lvl_f !! Level of fine grid real(FP), dimension(:), intent(out) :: uf !! Quantity on fine grid ! Should have dimension self%get_np(lvl_f) end subroutine module subroutine extrapolate_boundary_points(self, lvl, u) !! Extrapolates boundary points according to homegeneous Neumann !! boundary condition class(multigrid_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: lvl !! Level of mesh real(FP), dimension(:), intent(inout) :: u !! Quantity end subroutine module subroutine expose_data(self, intermediate_object, data_object) !! Exposes multigrid data through multigrid_data object !! boundary condition class(multigrid_t), intent(in), target :: self !! Instance of the type type(multigrid_intermediate_data_t), intent(out), target :: intermediate_object !! Intermediate object, holds arrays of mesh_cart and csrmat "data" type(multigrid_data_t), intent(out) :: data_object !! Destination multigrid_data object end subroutine module subroutine create_mrcsr(self, order, dbgout) !! Creates the restriction matrix mappings mrcsr between all !! consecutive mesh levels !! - Projection method based on wighted 9 point average is used !! - Only inner grid points are taken into account class(multigrid_t), target, intent(inout) :: self !! Instance of the type integer, intent(in) :: order !! Order of restriction !! 0: trivial restriction !! 1: weighted restriction (9-point stencil) integer, intent(in) :: dbgout !! Debug output level end subroutine module subroutine create_mrcsr_inner(self, dbgout) !! Creates restriction matrix mrcsr_inner purely based on inner !! grid points !! - Based on conventional restriction matrix previously created by !! create_mrcsr class(multigrid_t), target, intent(inout) :: self !! Instance of the type integer, intent(in) :: dbgout !! Debug output level end subroutine module subroutine create_mpcsr(self, dbgout) !! Creates the prolongation matrix mappings mpcsr between all !! consecutive mesh levels !! - Based on bilinear interpolation !! - Only inner grid points can be set via resulting matrix class(multigrid_t), target, intent(inout) :: self !! Instance of the type integer, intent(in) :: dbgout !! Debug output level end subroutine end interface contains subroutine set_phi_coarse(self, new_phi) !! Sets phi value on coarse meshes to a new value !! On finest mesh that is pointed to, needs to be set separately class(multigrid_t), target, intent(inout) :: self !! Instance of the type real(FP), intent(in) :: new_phi !! New phi value integer :: lvl if (allocated(self%mesh_coarse_array)) then do lvl = 2, self%nlvls call self%mesh_coarse_array(lvl)%set_phi(new_phi) enddo else write(get_stderr(), *) & 'error(multigrid%set_phi_coarse):: meshes not yet created' error stop endif end subroutine subroutine copy(self, multigrid_copy, mesh_finest) !! Copies the type to multigrid_copy, setting the internal pointer to !! mesh_finest class(multigrid_t), intent(in) :: self !! Instance of the type type(multigrid_t), intent(out) :: multigrid_copy !! Copied multigrid type(mesh_cart_t), target, intent(in) :: mesh_finest !! Finest (lvl=1) mesh to point to multigrid_copy = self multigrid_copy%mesh_finest_ptr => mesh_finest end subroutine function get_mesh_pointer(self, lvl) result(ptr) !! Returns mesh on desired level class(multigrid_t), target, intent(inout) :: self !! Instance of the type integer, intent(in) :: lvl !! Multigrid level of desired mesh type(mesh_cart_t), pointer :: ptr !! Pointer to desired mesh if (allocated(self%mesh_coarse_array) .and. & associated(self%mesh_finest_ptr)) then if (lvl == 1) then ptr => self%mesh_finest_ptr else ptr => self%mesh_coarse_array(lvl) end if else write(get_stderr(), *) & 'error(multigrid%get_mesh_pointer):: meshes not yet created' error stop endif end function pure integer function get_nlvls(self) !! Returns nlvls class(multigrid_t), intent(in) :: self !! Instance of the type get_nlvls = self%nlvls end function pure integer function get_np(self, lvl) !! Returns the number of total mesh points at the given level class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Desired multigrid level get_np = self%np_lvl(lvl) end function pure integer function get_np_boundary(self, lvl) !! Returns the number of boundary mesh points at the given level class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Desired multigrid level get_np_boundary = self%np_boundary_lvl(lvl) end function pure integer function get_np_inner(self, lvl) !! Returns the number of inner mesh points at the given level class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Desired multigrid level get_np_inner = self%np_inner_lvl(lvl) end function function get_first_idx_lvl_pointer(self) result(ptr) !! Returns a pointer to the first_idx_lvl array class(multigrid_t), target, intent(in) :: self !! Instance of the type integer, pointer, dimension(:) :: ptr ptr => self%first_idx_lvl end function pure function get_first_idx(self, lvl) result(res) !! Returns first_idx at lvl class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Level (up to nlvl+1) integer :: res res = self%first_idx_lvl(lvl) end function function get_first_idx_boundary_lvl_pointer(self) result(ptr) !! Returns a pointer to the first_idx_boundary_lvl array class(multigrid_t), target, intent(in) :: self !! Instance of the type integer, pointer, dimension(:) :: ptr ptr => self%first_idx_boundary_lvl end function pure function get_first_idx_boundary(self, lvl) result(res) !! Returns first_idx_boundary at lvl class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Level (up to nlvl+1) integer :: res res = self%first_idx_boundary_lvl(lvl) end function function get_first_idx_inner_lvl_pointer(self) result(ptr) !! Returns a pointer to the first_idx_inner_lvl array class(multigrid_t), target, intent(in) :: self !! Instance of the type integer, pointer, dimension(:) :: ptr ptr => self%first_idx_inner_lvl end function pure function get_first_idx_inner(self, lvl) result(res) !! Returns first_idx_inner at lvl class(multigrid_t), intent(in) :: self !! Instance of the type integer, intent(in) :: lvl !! Level (up to nlvl+1) integer :: res res = self%first_idx_inner_lvl(lvl) end function pure function get_np_total(self) result(res) !! Returns number of mesh points sumed over all levels class(multigrid_t), intent(in) :: self !! Instance of the type integer :: res res = self%get_first_idx(self%nlvls + 1) - 1 end function pure function get_np_inner_total(self) result(res) !! Returns number of inner mesh points sumed over all levels class(multigrid_t), intent(in) :: self !! Instance of the type integer :: res res = self%get_first_idx_inner(self%nlvls + 1) - 1 end function function get_mrcsr_pointer(self, lvl_f) result(res) !! Returns pointer to restriction matrix class(multigrid_t), intent(in), target :: self !! Instance of the type integer, intent(in) :: lvl_f !! Level of restriction matrix (from lvl_f to lvl_f+1) type(csrmat_t), pointer :: res !! Helmholtz matrix for level lvl res => self%mrcsr(lvl_f) end function function get_mpcsr_pointer(self, lvl_f) result(res) !! Returns pointer to prolongation matrix class(multigrid_t), intent(in), target :: self !! Instance of the type integer, intent(in) :: lvl_f !! Level of prolongation matrix (from lvl_f+1 to lvl_f) type(csrmat_t), pointer :: res !! Helmholtz matrix for level lvl res => self%mpcsr(lvl_f) end function subroutine destructor_multigrid(self) !! Frees memory associated with multigrid type(multigrid_t), intent(inout) :: self !! Instance of the type ! Actual data if (associated(self%mesh_finest_ptr)) self%mesh_finest_ptr => null() if (allocated(self%mesh_coarse_array)) & deallocate(self%mesh_coarse_array) if (allocated(self%mrcsr)) deallocate(self%mrcsr) if (allocated(self%mrcsr_inner)) deallocate(self%mrcsr_inner) if (allocated(self%mpcsr)) deallocate(self%mpcsr) if (allocated(self%np_inner_lvl)) deallocate(self%np_inner_lvl) if (allocated(self%np_boundary_lvl)) deallocate(self%np_boundary_lvl) if (allocated(self%np_lvl)) deallocate(self%np_lvl) if (allocated(self%first_idx_inner_lvl)) & deallocate(self%first_idx_inner_lvl) if (allocated(self%first_idx_boundary_lvl)) & deallocate(self%first_idx_boundary_lvl) if (allocated(self%first_idx_lvl)) deallocate(self%first_idx_lvl) self%nlvls = 0 end subroutine subroutine destructor_multigrid_intermediate_data(self) !! Frees memory associated with multigrid_intermediate_data type(multigrid_intermediate_data_t), intent(inout) :: self !! Instance of the type if (allocated(self%mesh_coarse_array_data)) & deallocate(self%mesh_coarse_array_data) if (allocated(self%mrcsr_data)) deallocate(self%mrcsr_data) if (allocated(self%mrcsr_inner_data)) deallocate(self%mrcsr_inner_data) if (allocated(self%mpcsr_data)) deallocate(self%mpcsr_data) end subroutine end module