module multigrid_solver_m !! Multigrid solver for the Helmholtz problem use, intrinsic :: iso_c_binding, only: c_ptr, c_loc use netcdf use precision_m, only : FP use descriptors_m, only : DISTRICT_CORE, DISTRICT_WALL, & DISTRICT_DOME, DISTRICT_OUT use csrmat_m, only : csrmat_t, csrmat_data_t, csr_times_vec, csr_residuum use mesh_cart_m, only : mesh_cart_t use multigrid_m, only : multigrid_t, multigrid_intermediate_data_t, multigrid_data_t use direct_solver_m, only : direct_solver_t use boundaries_perp_m, only : set_perpbnds, extrapolate_ghost_points use helmholtz_operator_m, only : build_helmholtz_csr use splitting_m, only : splitting_t, splitting_jacobi_cpu_t, & splitting_gauss_seidel_cpu_t, & splitting_gauss_seidel_redblack_cpu_t type, public :: multigrid_solver_t !! Multigrid solver object type(multigrid_t), private, pointer :: multigrid => null() !! Multigrid on which to solve the Helmholtz problem type(csrmat_t), allocatable, dimension(:), private :: hcsrs !! Helmholtz matrices real(FP), allocatable, dimension(:), private :: hdiags_inv !! Inverse of diagonals of Helmholtz matrices integer, allocatable, dimension(:), private :: bnd_descrs !! Boundary descriptors integer, private :: ncycle !! Number of cycles (1: V-cycle, 2: W-cycle) class(splitting_t), allocatable, private, dimension(:) :: msmoother !! Direct solver for coarsest level integer, private :: npresmooth !! Number of presmoothing steps integer, private :: npostsmooth !! Dumber of postsmoothing steps class(direct_solver_t), allocatable, private :: coarsest_solver !! Direct solver for coarsest level contains procedure, public :: create => create_multigrid_solver procedure, public :: update => update_multigrid_solver procedure, public :: init => initialize_multigrid_solver procedure, public :: cycle procedure, public :: residuum procedure, public :: expose_multigrid_data procedure, public :: expose_hcsr procedure, public :: expose_hdiags_inv procedure, public :: get_hcsr_finest_pointer final :: destructor_multigrid_solver end type interface module subroutine create_multigrid_solver(self, multigrid, & bnd_type_core, bnd_type_wall, bnd_type_dome, & bnd_type_out, co, lambda, xi ) !! Creates multigrid solver class(multigrid_solver_t), intent(inout) :: self !! Instance of the type type(multigrid_t), target, intent(in) :: multigrid !! Multigrid on which to solve the Helmholtz problem integer, intent(in) :: bnd_type_core !! Boundary descriptor for core boundary integer, intent(in) :: bnd_type_wall !! Boundary descriptor for wall boundary integer, intent(in) :: bnd_type_dome !! Boundary descriptor for dome boundary integer, intent(in) :: bnd_type_out !! Boundary descriptor for outer(mask) boundary real(FP), dimension(multigrid%get_np(1)), intent(in) :: co !! Coefficient co real(FP), dimension(multigrid%get_np_inner(1)), intent(in) :: lambda !! Coefficient lambda real(FP), dimension(multigrid%get_np_inner(1)), intent(in) :: xi !! Coefficient xi end subroutine module subroutine initialize_multigrid_solver(self, ncycle, smoother, & npresmooth, npostsmooth, dirsolver) !! Initialises multigrid solver with parameters and solvers class(multigrid_solver_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: ncycle !! Number of multigrid cycles 1=V, 2=W class(splitting_t), allocatable, intent(inout) :: smoother !! Smoother !! On output the object is deallocated integer, intent(in) :: npresmooth !! Number of presmoothing steps integer, intent(in) :: npostsmooth !! Number of postsmoothing steps class(direct_solver_t), allocatable, intent(inout) :: dirsolver !! Direct solver, used on coarsest level !! On output the object is deallocated end subroutine module subroutine update_multigrid_solver(self, co, lambda, xi, & bnd_type_core, & bnd_type_wall, & bnd_type_dome, bnd_type_out) !! Updates multigrid solver with new values of coefficients class(multigrid_solver_t), intent(inout) :: self !! Instance of the type real(FP), dimension(self%multigrid%get_np(1)), intent(in) :: co !! Coefficient co real(FP), dimension(self%multigrid%get_np_inner(1)), & intent(in) :: lambda !! Coefficient lambda real(FP), dimension(self%multigrid%get_np_inner(1)), & intent(in) :: xi !! Coefficient xi integer, intent(in), optional :: bnd_type_core !! Boundary descriptor for core boundary integer, intent(in), optional :: bnd_type_wall !! Boundary descriptor for wall boundary integer, intent(in), optional :: bnd_type_dome !! Boundary descriptor for dome boundary integer, intent(in), optional :: bnd_type_out !! Boundary descriptor for outer(mask) boundary end subroutine recursive module subroutine cycle(self, u, b, lvl) !! Performs a multigrid cycle class(multigrid_solver_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: lvl !! Level real(FP), dimension(:), intent(inout) :: u !! On input initial guess, on output result after one cycle real(FP), dimension(:), intent(inout) :: b !! Right hand side end subroutine module subroutine residuum(self, u, b, res, resmax) !! Computes the residuum ||Hu-b|| / ||b|| class(multigrid_solver_t), intent(in) :: self !! Instance of the type real(FP), dimension(:), intent(inout) :: u !! Guess for solution real(FP), dimension(:), intent(inout) :: b !! Right hand side real(FP), intent(out) :: res !! Residuum (root mean square) : ||Hu-b|| / ||b||+eps real(FP), intent(out), optional :: resmax !! Residuum maximum deviation max(Hu-b) end subroutine module subroutine expose_multigrid_data(self, intermediate_object, data_object) !! Calls multigrid expose data class(multigrid_solver_t), intent(in) :: self !! Instance of the type type(multigrid_intermediate_data_t), intent(out) :: intermediate_object !! Intermediate object type(multigrid_data_t), intent(out) :: data_object !! Output object end subroutine module subroutine expose_hcsr(self, lvl, res) !! Exposes the Helmholtz matrix for given level !! Required for the GPU solver. class(multigrid_solver_t), intent(in):: self !! Instance of the type integer, intent(in) :: lvl !! Level for which matrix shall be returned type(csrmat_data_t), intent(out) :: res !! Helmholtz matrix for level lvl end subroutine module subroutine expose_hdiags_inv(self, res) !! Exposes the inverse diagonal of Helmholtz matrices !! Required for the GPU solver. class(multigrid_solver_t), intent(in), target :: self !! Instance of the type type(c_ptr), intent(out) :: res !! Pointer to inverse diagonal of Helmhotz matrices end subroutine end interface contains function get_hcsr_finest_pointer(self) result(ptr) !! Returns pointer to Helmholtz matrix on finest level class(multigrid_solver_t), intent(in), target :: self !! Instance of the type type(csrmat_t), pointer :: ptr !! Pointer to Helmholtz matrix on finest level ptr => self%hcsrs(1) end function subroutine destructor_multigrid_solver(self) !! Frees memory associated with multigrid solver type(multigrid_solver_t), intent(inout) :: self !! Instance of the type if (allocated(self%hcsrs)) deallocate(self%hcsrs) if (allocated(self%coarsest_solver)) deallocate(self%coarsest_solver) if (allocated(self%hdiags_inv)) deallocate(self%hdiags_inv) if (allocated(self%bnd_descrs)) deallocate(self%bnd_descrs) if (allocated(self%msmoother)) deallocate(self%msmoother) if (associated(self%multigrid)) self%multigrid => null() self%npresmooth = 0 self%npostsmooth = 0 self%ncycle = 0 end subroutine end module