submodule (multigrid_solver_m) multigrid_solver_s use screen_io_m, only : get_stderr use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_HELMHOLTZ use boundaries_perp_m, only : bnd_types_to_bnd_descrs implicit none contains module subroutine create_multigrid_solver(self, multigrid, & bnd_type_core, bnd_type_wall, bnd_type_dome, & bnd_type_out, co, lambda, xi ) class(multigrid_solver_t), intent(inout) :: self type(multigrid_t), target, intent(in) :: multigrid integer, intent(in) :: bnd_type_core integer, intent(in) :: bnd_type_wall integer, intent(in) :: bnd_type_dome integer, intent(in) :: bnd_type_out real(FP), dimension(multigrid%get_np(1)), intent(in) :: co real(FP), dimension(multigrid%get_np_inner(1)), intent(in) :: lambda real(FP), dimension(multigrid%get_np_inner(1)), intent(in) :: xi integer :: nlvls, lvl, kb, l, district integer :: pt_all_i, pt_boundary_i, pt_inner_i, ptb integer, dimension(:), pointer :: pfi, pfi_i, pfi_b => null() type(mesh_cart_t), pointer :: mmesh => null() real(FP), dimension(multigrid%get_np_total()) :: mco real(FP), dimension(multigrid%get_np_inner_total()) :: mlambda, mxi self%multigrid => multigrid nlvls = self%multigrid%get_nlvls() pfi => self%multigrid%get_first_idx_lvl_pointer() pfi_i => self%multigrid%get_first_idx_inner_lvl_pointer() pfi_b => self%multigrid%get_first_idx_boundary_lvl_pointer() ! Allocation allocate(self%hcsrs(nlvls)) allocate(self%hdiags_inv(pfi(nlvls + 1) - 1)) allocate(self%bnd_descrs(pfi_b(nlvls + 1) - 1)) ! Build boundary descriptor for quantity on each level do lvl = 1, nlvls mmesh => self%multigrid%get_mesh_pointer(lvl) call bnd_types_to_bnd_descrs(mmesh, & bnd_type_core, bnd_type_wall, bnd_type_dome, bnd_type_out, & self%bnd_descrs(pfi_b(lvl):pfi_b(lvl+1)-1) ) enddo call restrict_coefficients(self%multigrid, co, lambda, xi, & mco, mlambda, mxi) ! Build Helmholtz matrix on each level do lvl = 1, nlvls mmesh => self%multigrid%get_mesh_pointer(lvl) pt_all_i = pfi(lvl) pt_inner_i = pfi_i(lvl) pt_boundary_i = pfi_b(lvl) call build_helmholtz_csr(mmesh, & self%hcsrs(lvl), & self%hdiags_inv(pt_all_i), & mco(pt_all_i), & mxi(pt_inner_i), & mlambda(pt_inner_i), & self%bnd_descrs(pt_boundary_i)) enddo mmesh => null() pfi => null() pfi_i => null() pfi_b => null() end subroutine module subroutine initialize_multigrid_solver(self, ncycle, smoother, & npresmooth, npostsmooth, dirsolver) class(multigrid_solver_t), intent(inout) :: self integer, intent(in) :: ncycle class(splitting_t), allocatable, intent(inout) :: smoother integer, intent(in) :: npresmooth integer, intent(in) :: npostsmooth class(direct_solver_t), allocatable, intent(inout) :: dirsolver integer :: nlvls, lvl type(mesh_cart_t), pointer :: mmesh nlvls = self%multigrid%get_nlvls() ! Set parameters self%ncycle = ncycle self%npresmooth = npresmooth self%npostsmooth = npostsmooth ! Create smoothers on each level if (.not.allocated(smoother)) then call handle_error('smoother provided not allocated', & PARALLAX_ERR_HELMHOLTZ, __LINE__, __FILE__) endif allocate(self%msmoother(nlvls), mold=smoother) deallocate(smoother) do lvl = 1, nlvls mmesh => self%multigrid%get_mesh_pointer(lvl) call self%msmoother(lvl)%create(mmesh) enddo ! Create direct solver on coarsest level if (.not. allocated(dirsolver)) then call handle_error('direct solver provided not allocated', & PARALLAX_ERR_HELMHOLTZ, __LINE__, __FILE__) endif call move_alloc(dirsolver, self%coarsest_solver) call self%coarsest_solver%create(self%hcsrs(nlvls)) call self%coarsest_solver%factorise(self%hcsrs(nlvls),& 'ANALYSE_AND_NUMFAC') end subroutine module subroutine update_multigrid_solver(self, co, lambda, xi, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out) class(multigrid_solver_t), intent(inout) :: self real(FP), dimension(self%multigrid%get_np(1)), intent(in) :: co real(FP), dimension(self%multigrid%get_np_inner(1)), & intent(in) :: lambda real(FP), dimension(self%multigrid%get_np_inner(1)), & intent(in) :: xi integer, intent(in), optional :: bnd_type_core integer, intent(in), optional :: bnd_type_wall integer, intent(in), optional :: bnd_type_dome integer, intent(in), optional :: bnd_type_out logical :: update_with_bnds integer :: nlvls, lvl, pt_all_i, pt_inner_i, pt_bndy_i, pt_bndy_i_next type(mesh_cart_t), pointer :: mmesh => null() real(FP), dimension(self%multigrid%get_np_total()) :: mco real(FP), dimension(self%multigrid%get_np_inner_total()) :: mlambda, mxi update_with_bnds = present(bnd_type_core) & .and. present(bnd_type_wall) & .and. present(bnd_type_dome) & .and. present(bnd_type_out) nlvls = self%multigrid%get_nlvls() call restrict_coefficients(self%multigrid, co, lambda, xi, & mco, mlambda, mxi) if (update_with_bnds) then ! Build boundary descriptor for quantity on each level do lvl = 1, nlvls mmesh => self%multigrid%get_mesh_pointer(lvl) pt_bndy_i = self%multigrid%get_first_idx_boundary(lvl) pt_bndy_i_next = self%multigrid%get_first_idx_boundary(lvl+1) call bnd_types_to_bnd_descrs(mmesh, & bnd_type_core, bnd_type_wall, bnd_type_dome, bnd_type_out, & self%bnd_descrs(pt_bndy_i:pt_bndy_i_next-1) ) enddo endif ! Build Helmholtz matrix on each level do lvl = 1, nlvls mmesh => self%multigrid%get_mesh_pointer(lvl) pt_all_i = self%multigrid%get_first_idx(lvl) pt_inner_i = self%multigrid%get_first_idx_inner(lvl) pt_bndy_i = self%multigrid%get_first_idx_boundary(lvl) if (update_with_bnds) then call build_helmholtz_csr(mmesh, & self%hcsrs(lvl), & self%hdiags_inv(pt_all_i), & mco(pt_all_i), & mxi(pt_inner_i), & mlambda(pt_inner_i), & self%bnd_descrs(pt_bndy_i) ) else call build_helmholtz_csr(mmesh, & self%hcsrs(lvl), & self%hdiags_inv(pt_all_i), & mco(pt_all_i), & mxi(pt_inner_i), & mlambda(pt_inner_i) ) endif enddo ! Factorise matrix on coarsest level if (allocated(self%coarsest_solver)) then ! This if clause is useful if the update routine is called without ! having called the init routine previously. This is useful, if ! one wants just to access the matrices, e.g for the PETSc solver. call self%coarsest_solver%factorise(self%hcsrs(nlvls), 'NUMFAC') endif mmesh => null() end subroutine recursive module subroutine cycle(self, u, b, lvl) class(multigrid_solver_t), intent(inout) :: self integer, intent(in) :: lvl real(FP), dimension(:), intent(inout) :: u real(FP), dimension(:), intent(inout) :: b integer :: l, k, nlvls, lvl_c, ismooth, pt_all_i, pt_boundary_i integer, dimension(:), pointer :: pfi, pfi_b type(mesh_cart_t), pointer :: mmesh real(FP), dimension(:), allocatable :: df, ef, dc, ec nlvls = self%multigrid%get_nlvls() pfi => self%multigrid%get_first_idx_lvl_pointer() pfi_b => self%multigrid%get_first_idx_boundary_lvl_pointer() mmesh => self%multigrid%get_mesh_pointer(lvl) ! Check dimension of arrays (for stupid issue with intel/19 compiler) if ( (size(u) /= self%multigrid%get_np(lvl)) .or. & (size(b) /= self%multigrid%get_np(lvl)) ) then write(get_stderr(), *) & 'error(multigrid_solver%cycle): size of u or b wrong' stop endif allocate(df(self%multigrid%get_np(lvl))) allocate(ef, mold=df) allocate(dc(self%multigrid%get_np(min(lvl + 1, nlvls)))) allocate(ec, mold=dc) if (lvl == nlvls) then ! Apply direct solver on coarsest level call self%coarsest_solver%backsub(self%hcsrs(lvl), b, u) else lvl_c = lvl + 1 pt_all_i = pfi(lvl) pt_boundary_i = pfi_b(lvl) ! Presmoothing !$omp parallel private(l, ismooth) do ismooth = 1, self%npresmooth call self%msmoother(lvl)%apply(a=self%hcsrs(lvl), & dinv=self%hdiags_inv(pt_all_i), & x=u, & b=b ) enddo ! Compute defect call csr_times_vec(self%hcsrs(lvl), u, df) !$omp do do l = 1, self%multigrid%get_np(lvl) df(l) = df(l) - b(l) enddo !$omp end do nowait !$omp do do l = 1, self%multigrid%get_np(lvl_c) ec(l) = 0.0_FP enddo !$omp end do ! Restrict defect call self%multigrid%restrict(lvl, df, lvl_c, dc) !$omp end parallel ! Recursive multigrid call on defect do k = 1, self%ncycle call self%cycle(u=ec, b=dc, lvl=lvl_c) enddo !$omp parallel ! Prolongate error call self%multigrid%prolong(lvl_c, ec, lvl, ef) call set_perpbnds(mmesh, self%bnd_descrs(pt_boundary_i), ef) ! Correct solution !$omp do do l = 1, self%multigrid%get_np(lvl) u(l) = u(l) - ef(l) enddo !$omp end do ! Postmoothing do ismooth = 1, self%npostsmooth call self%msmoother(lvl)%apply(a=self%hcsrs(lvl), & dinv=self%hdiags_inv(pt_all_i), & x=u, & b=b ) enddo !$omp end parallel endif end subroutine module subroutine residuum(self, u, b, res, resmax) class(multigrid_solver_t), intent(in) :: self real(FP), dimension(:), intent(inout) :: u real(FP), dimension(:), intent(inout) :: b real(FP), intent(out) :: res real(FP), intent(out), optional :: resmax real(FP) :: resmax_work ! Check dimension of arrays (for stupid issue with intel/19 compiler) if ( (size(u) /= self%multigrid%get_np(1)) .or. & (size(b) /= self%multigrid%get_np(1)) ) then write(get_stderr(), *) & 'error(multigrid_solver%residuum): size of u or b wrong' stop endif call csr_residuum(self%hcsrs(1), u, b, res, resmax_work) if (present(resmax)) then resmax = resmax_work endif end subroutine module subroutine expose_multigrid_data(self, intermediate_object, data_object) class(multigrid_solver_t), intent(in) :: self type(multigrid_intermediate_data_t), intent(out) :: intermediate_object type(multigrid_data_t), intent(out) :: data_object call self%multigrid%expose_data(intermediate_object, data_object) end subroutine module subroutine expose_hcsr(self, lvl, res) class(multigrid_solver_t), intent(in):: self integer, intent(in) :: lvl type(csrmat_data_t), intent(out) :: res call self%hcsrs(lvl)%expose_data(res) end subroutine module subroutine expose_hdiags_inv(self, res) class(multigrid_solver_t), intent(in), target :: self type(c_ptr), intent(out) :: res res = c_loc(self%hdiags_inv) end subroutine subroutine restrict_coefficients(multigrid, co, lambda, xi, & mco, mlambda, mxi) type(multigrid_t), intent(inout) :: multigrid !! Multigrid real(FP), intent(in), dimension(multigrid%get_np(1)) :: co !! Coefficient co (on finest mesh) real(FP), intent(in), dimension(multigrid%get_np_inner(1)) :: lambda !! Coefficient lambda (on finest mesh) real(FP), intent(in), dimension(multigrid%get_np_inner(1)) :: xi !! Coefficient xi (on finest mesh) real(FP), & dimension(multigrid%get_first_idx(multigrid%get_nlvls()+1)-1), & intent(out) :: mco !! Coefficient co restricted to all levels real(FP), & dimension( & multigrid%get_first_idx_inner(multigrid%get_nlvls()+1)-1), & intent(out) :: mlambda !! Coefficient lambda restricted to all levels real(FP), & dimension(& multigrid%get_first_idx_inner(multigrid%get_nlvls()+1)-1), & intent(out) :: mxi !! Coefficient xi restricted to all levels integer :: lvl, nlvls, pt_all_i, pt_inner_i, l, ki, pta, pti type(mesh_cart_t), pointer :: mmesh => null() integer, dimension(:), pointer :: pfi, pfi_i, pfi_b nlvls = multigrid%get_nlvls() pfi => multigrid%get_first_idx_lvl_pointer() pfi_i => multigrid%get_first_idx_inner_lvl_pointer() pfi_b => multigrid%get_first_idx_boundary_lvl_pointer() do lvl = 1, nlvls mmesh => multigrid%get_mesh_pointer(lvl) pt_all_i = pfi(lvl) pt_inner_i = pfi_i(lvl) if (lvl == 1) then !$omp parallel default(none) private(l, ki, pta, pti) & !$omp shared(multigrid, lvl, pt_all_i, pt_inner_i, & !$omp co, mco, lambda, mlambda, xi, mxi) !$omp do do l = 1, multigrid%get_np(lvl) pta = l - 1 + pt_all_i mco(pta) = co(l) enddo !$omp end do !$omp do do ki = 1, multigrid%get_np_inner(lvl) pti = ki - 1 + pt_inner_i mlambda(pti) = lambda(ki) mxi(pti) = xi(ki) enddo !$omp end do !$omp end parallel else call multigrid%restrict( & lvl-1, & mco(pfi(lvl-1):pfi(lvl)-1), & lvl, & mco(pfi(lvl):pfi(lvl+1)-1) ) ! Extrapolate boundary and ghost points call multigrid%extrapolate_boundary_points( & lvl, & mco(pfi(lvl):pfi(lvl+1)-1) ) call extrapolate_ghost_points( & mmesh, & mco(pfi(lvl):pfi(lvl+1)-1) ) call multigrid%restrict_inner( & lvl-1,& mlambda(pfi_i(lvl-1):pfi_i(lvl)-1), & lvl, & mlambda(pfi_i(lvl):pfi_i(lvl+1)-1) ) call multigrid%restrict_inner( & lvl-1, & mxi(pfi_i(lvl-1):pfi_i(lvl)-1), & lvl, & mxi(pfi_i(lvl):pfi_i(lvl+1)-1) ) endif enddo mmesh => null() end subroutine end submodule