submodule (multigrid_m) multigrid_s use screen_io_m, only : get_stdout, nf90_handle_err implicit none contains module subroutine initialize_multigrid(self, equi, phi, nlvls, & spacing_f, size_neighbor, & size_ghost_layer, mesh_finest, & reorder_size, extend_beyond_wall, & dbgout) class(multigrid_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi real(FP), intent(in) :: phi integer, intent(in) :: nlvls real(FP), intent(inout) :: spacing_f integer, intent(in) :: size_neighbor integer, intent(in) :: size_ghost_layer type(mesh_cart_t), target, intent(out) :: mesh_finest integer, intent(in), dimension(nlvls), optional :: reorder_size logical, intent(in) :: extend_beyond_wall integer, intent(in), optional :: dbgout integer :: outi integer :: lvl, sngb, sghl, lvl_f, lvl_c ! Set output level outi = 0 if (present(dbgout)) then if (is_master()) then outi = dbgout endif endif if (outi >= 1) then write(get_stdout(), *) '' write(get_stdout(), *) & '------------------------------------------------' write(get_stdout(), 201) nlvls, spacing_f, & size_neighbor, size_ghost_layer,extend_beyond_wall 201 FORMAT('Setting up multigrid' & 3X, 'nlvls = ', I3, / & 3X, 'spacing_f = ', ES14.6E3, / & 3X, 'size_neighbor = ', I3, / & 3X, 'size_ghost_layer = ', I3, / & 3X, 'extend_beyond_wall = ',L1 ) endif ! Set number of grid levels self%nlvls = nlvls ! Create meshes on all levels allocate(self%mesh_coarse_array(2:self%nlvls)) do lvl = 1, self%nlvls if (lvl == 1) then sngb = size_neighbor sghl = size_ghost_layer call mesh_finest%init(equi, phi, lvl, spacing_f, & sngb, sghl, extend_beyond_wall, dbgout) if (present(reorder_size)) then call mesh_finest%reorder(reorder_size(lvl), dbgout) endif self%mesh_finest_ptr => mesh_finest else sngb = 1 sghl = 0 call self%mesh_coarse_array(lvl)%init(equi, phi, lvl, & spacing_f, sngb, & sghl, & extend_beyond_wall, & dbgout) if (present(reorder_size)) then call self%mesh_coarse_array(lvl)%reorder( & reorder_size(lvl), dbgout) endif endif enddo ! Create n_points arrays (regular and consecutive) call self%initialize_np_and_idx_arrays() ! Create prolongation and restriction matrices allocate(self%mrcsr(self%nlvls - 1)) allocate(self%mrcsr_inner(self%nlvls - 1)) allocate(self%mpcsr(self%nlvls - 1)) call self%create_mrcsr(order=1, dbgout=outi) call self%create_mrcsr_inner(dbgout=outi) call self%create_mpcsr(dbgout=outi) do lvl = 2, self%nlvls ! Deallocate the array used to accelerate cart_to_mesh_index ! on the coarse grids, since these won't be used again call self%mesh_coarse_array(lvl)%deallocate_shaped_cart_arrays() enddo if (outi >= 1) then write(get_stdout(), *) 'Multigrid succesfully set up:' write(get_stdout(), *) & '------------------------------------------------' write(get_stdout(), *) '' endif end subroutine module subroutine initialize_np_and_idx_arrays(self) class(multigrid_t), target, intent(inout) :: self integer :: lvl type(mesh_cart_t), pointer :: mmesh integer, dimension(:), pointer :: np, np_b, np_i, fi, fi_b, fi_i allocate(self%np_lvl(self%nlvls)) allocate(self%np_boundary_lvl(self%nlvls)) allocate(self%np_inner_lvl(self%nlvls)) allocate(self%first_idx_lvl(self%nlvls + 1)) allocate(self%first_idx_boundary_lvl(self%nlvls + 1)) allocate(self%first_idx_inner_lvl(self%nlvls + 1)) np => self%np_lvl np_b => self%np_boundary_lvl np_i => self%np_inner_lvl fi => self%first_idx_lvl fi_b => self%first_idx_boundary_lvl fi_i => self%first_idx_inner_lvl fi(1) = 1 fi_b(1) = 1 fi_i(1) = 1 do lvl = 1, self%nlvls mmesh => self%get_mesh_pointer(lvl) np(lvl) = mmesh%get_n_points() np_b(lvl) = mmesh%get_n_points_boundary() np_i(lvl) = mmesh%get_n_points_inner() fi(lvl + 1) = fi(lvl) + np(lvl) fi_b(lvl + 1) = fi_b(lvl) + np_b(lvl) fi_i(lvl + 1) = fi_i(lvl) + np_i(lvl) enddo end subroutine module subroutine multigrid_write_netcdf(self, fgid) class(multigrid_t), intent(inout) :: self integer, intent(in) :: fgid type(mesh_cart_t), pointer :: mmesh integer :: nf90_stat, lvl, mesh_id, rest_id, rest_inner_id, prol_id character(len=3) :: lvl_char, lvl_char_c, lvl_char_f character(len=12) :: mesh_group_name character(len=26) :: rest_group_name character(len=32) :: rest_inner_group_name character(len=27) :: prol_group_name character(len=*), parameter :: cname = "multigrid_write_netcdf" ! Name used for logging errors ! Write global attributes nf90_stat = nf90_put_att(fgid, NF90_GLOBAL, 'nlvls', self%nlvls) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Write meshes (each level in its own group) do lvl = 1, self%nlvls write(lvl_char, '(I3.3)')lvl mesh_group_name = 'mesh_lvl_' // lvl_char ! Create group nf90_stat = nf90_def_grp(fgid, mesh_group_name, mesh_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) mmesh => self%get_mesh_pointer(lvl) call mmesh%write_netcdf(mesh_id) enddo ! Write prolongation and restriction matrices (each level its own ! group) do lvl = 1, self%nlvls - 1 write(lvl_char_f, '(I3.3)')lvl write(lvl_char_c, '(I3.3)')lvl + 1 rest_group_name = 'restriction_lvl_' // lvl_char_f // & &'_to_' // lvl_char_c rest_inner_group_name = 'restriction_inner_lvl_' // lvl_char_f // & &'_to_' // lvl_char_c prol_group_name = 'prolongation_lvl_' // lvl_char_c // & &'_to_' // lvl_char_f ! Create group nf90_stat = nf90_def_grp(fgid, rest_group_name, rest_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_grp(fgid, rest_inner_group_name, & rest_inner_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_def_grp(fgid, prol_group_name, prol_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) call self%mrcsr(lvl)%write_netcdf(rest_id) call self%mrcsr_inner(lvl)%write_netcdf(rest_inner_id) call self%mpcsr(lvl)%write_netcdf(prol_id) enddo end subroutine module subroutine multigrid_read_netcdf(self, fgid, mesh_finest) class(multigrid_t), intent(inout) :: self integer, intent(in) :: fgid type(mesh_cart_t), target, intent(out) :: mesh_finest type(mesh_cart_t), pointer :: mmesh integer :: nf90_stat integer :: lvl, mesh_id, rest_id, rest_inner_id, prol_id character(len=3) :: lvl_char, lvl_char_c, lvl_char_f character(len=12) :: mesh_group_name character(len=26) :: rest_group_name character(len=32) :: rest_inner_group_name character(len=27) :: prol_group_name character(len=*), parameter :: cname = "multigrid_read_netcdf" ! Name used for logging errors ! Associate pointer to finest mesh if (.not. associated(self%mesh_finest_ptr)) then self%mesh_finest_ptr => mesh_finest else write(get_stderr(), *) & 'error(multigrid%read_netcdf): & &self%mesh_finest_ptr must not be associated' stop end if ! Read global attributes nf90_stat = nf90_get_att(fgid, NF90_GLOBAL, 'nlvls', self%nlvls) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) ! Allocate structures if (.not.allocated(self%mesh_coarse_array)) then allocate(self%mesh_coarse_array(2:self%nlvls)) else write(get_stderr(), *) & 'error(multigrid%read_netcdf): & &structures must not be allocated' stop endif if (.not. allocated(self%mrcsr)) then allocate(self%mrcsr(self%nlvls - 1)) else write(get_stderr(), *) & 'error(multigrid%read_netcdf): & &structures must not be allocated' stop endif if (.not. allocated(self%mrcsr_inner)) then allocate(self%mrcsr_inner(self%nlvls - 1)) else write(get_stderr(), *) & 'error(multigrid%read_netcdf): & &structures must not be allocated' stop endif if (.not. allocated(self%mpcsr)) then allocate(self%mpcsr(self%nlvls - 1)) else write(get_stderr(), *) & 'error(multigrid%read_netcdf): & &structures must not be allocated' stop endif ! Read meshes (each level in its own group) do lvl = 1, self%nlvls write(lvl_char, '(I3.3)')lvl mesh_group_name = 'mesh_lvl_' // lvl_char nf90_stat = nf90_inq_grp_ncid(fgid, mesh_group_name, mesh_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) mmesh => self%get_mesh_pointer(lvl) call mmesh%read_netcdf(mesh_id) enddo ! Read prolongation and restriction matrices (each level in its own ! group) do lvl = 1, self%nlvls - 1 write(lvl_char_f, '(I3.3)')lvl write(lvl_char_c, '(I3.3)')lvl + 1 rest_group_name = 'restriction_lvl_' // lvl_char_f & &// '_to_' // lvl_char_c rest_inner_group_name = 'restriction_inner_lvl_' // lvl_char_f & &// '_to_' // lvl_char_c prol_group_name = 'prolongation_lvl_' // lvl_char_c & &// '_to_' // lvl_char_f nf90_stat = nf90_inq_grp_ncid(fgid, rest_group_name, rest_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_inq_grp_ncid(fgid, rest_inner_group_name, & rest_inner_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) nf90_stat = nf90_inq_grp_ncid(fgid, prol_group_name, prol_id) if (nf90_stat /= 0) call nf90_handle_err(nf90_stat, cname) call self%mrcsr(lvl)%read_netcdf(rest_id) call self%mrcsr_inner(lvl)%read_netcdf(rest_inner_id) call self%mpcsr(lvl)%read_netcdf(prol_id) enddo ! Create n_points arrays (regular and consecutive) call self%initialize_np_and_idx_arrays() end subroutine module subroutine restrict(self, lvl_f, uf, lvl_c, uc) class(multigrid_t), intent(in) :: self integer, intent(in) :: lvl_f real(FP), dimension(:), intent(in) :: uf integer, intent(in) :: lvl_c real(FP), dimension(:), intent(out) :: uc if ( ((lvl_f + 1) /= lvl_c ) .or. (lvl_c > self%nlvls) ) then write(get_stderr(), *) & 'error(multigrid%restrict): levels not of compatible coarseness' error stop endif ! Check dimension of arrays (for compiler issue with explicit-shape ! arrays uf and uc) if ( (size(uf) /= self%get_np(lvl_f)) .or. & (size(uc) /= self%get_np(lvl_c)) ) then write(get_stderr(), *) & 'error(multigrid%restrict): size of uf or uc wrong' error stop endif call csr_times_vec(self%mrcsr(lvl_f), uf, uc) end subroutine module subroutine restrict_inner(self, lvl_f, uf, lvl_c, uc) class(multigrid_t), intent(inout) :: self integer, intent(in) :: lvl_f real(FP), dimension(:), intent(in) :: uf integer, intent(in) :: lvl_c real(FP), dimension(:), intent(out) :: uc if ( ((lvl_f + 1) /= lvl_c ) .or. (lvl_c > self%nlvls) ) then write(get_stderr(), *) & 'error(multigrid%restrict_inner): & &levels not of compatible coarseness' error stop endif ! Check dimension of arrays (for compiler issue with explicit-shape ! arrays uf and uc) if ( (size(uf) /= self%get_np_inner(lvl_f)) .or. & (size(uc) /= self%get_np_inner(lvl_c)) ) then write(get_stderr(), *) & 'error(multigrid%restrict_inner): size of uf or uc wrong' error stop endif call csr_times_vec(self%mrcsr_inner(lvl_f), uf, uc) end subroutine module subroutine prolong(self, lvl_c, uc, lvl_f, uf) class(multigrid_t), intent(inout) :: self integer, intent(in) :: lvl_c real(FP), dimension(:), intent(in) :: uc integer, intent(in) :: lvl_f real(FP), dimension(:), intent(out) :: uf if ( ((lvl_f + 1) /= lvl_c ) .or. (lvl_c > self%nlvls) ) then write(get_stderr(), *) & 'error(multigrid%prolong): levels not of compatible coarseness' error stop endif ! Check dimension of arrays (for compiler issue with explicit-shape ! arrays uf and uc) if ( (size(uc) /= self%get_np(lvl_c)) .or. & (size(uf) /= self%get_np(lvl_f)) ) then write(get_stderr(), *) & 'error(multigrid%prolong): size of uc or uf wrong' error stop endif call csr_times_vec(self%mpcsr(lvl_f), uc, uf) end subroutine module subroutine extrapolate_boundary_points(self, lvl, u) class(multigrid_t), intent(inout) :: self integer, intent(in) :: lvl real(FP), dimension(:), intent(inout) :: u type(mesh_cart_t), pointer :: mmesh integer :: kb integer, dimension(self%get_np_boundary(lvl)) :: bnd_descrs !$omp do do kb = 1, self%get_np_boundary(lvl) bnd_descrs(kb) = BND_TYPE_NEUMANN enddo !$omp end do mmesh => self%get_mesh_pointer(lvl) call set_perpbnds(mmesh, bnd_descrs, u) end subroutine module subroutine expose_data(self, intermediate_object, data_object) class(multigrid_t), intent(in), target :: self type(multigrid_intermediate_data_t), intent(out), target :: intermediate_object type(multigrid_data_t), intent(out) :: data_object integer ::lvl ! Allocate intermediate object arrays allocate(intermediate_object%mesh_coarse_array_data(2:self%nlvls)) allocate(intermediate_object%mrcsr_data(self%nlvls-1)) allocate(intermediate_object%mrcsr_inner_data(self%nlvls-1)) allocate(intermediate_object%mpcsr_data(self%nlvls-1)) ! Populate intermediate object call self%mesh_finest_ptr%expose_data(intermediate_object%mesh_finest_ptr_data) do lvl = 2, self%nlvls call self%mesh_coarse_array(lvl)%expose_data( & intermediate_object%mesh_coarse_array_data(lvl)) enddo do lvl = 1, self%nlvls - 1 call self%mrcsr(lvl)%expose_data( & intermediate_object%mrcsr_data(lvl)) call self%mrcsr_inner(lvl)%expose_data( & intermediate_object%mrcsr_inner_data(lvl)) call self%mpcsr(lvl)%expose_data( & intermediate_object%mpcsr_data(lvl)) enddo data_object%nlvls = self%nlvls data_object%mesh_finest_ptr = c_loc(intermediate_object%mesh_finest_ptr_data) data_object%mesh_coarse_array = c_loc(intermediate_object%mesh_coarse_array_data) data_object%mrcsr = c_loc(intermediate_object%mrcsr_data) data_object%mrcsr_inner = c_loc(intermediate_object%mrcsr_inner_data) data_object%mpcsr = c_loc(intermediate_object%mpcsr_data) data_object%np_lvl = c_loc(self%np_lvl) data_object%np_boundary_lvl = c_loc(self%np_boundary_lvl) data_object%np_inner_lvl = c_loc(self%np_inner_lvl) data_object%first_idx_lvl = c_loc(self%first_idx_lvl) data_object%first_idx_boundary_lvl = c_loc(self%first_idx_boundary_lvl) data_object%first_idx_inner_lvl = c_loc(self%first_idx_inner_lvl) end subroutine module subroutine create_mrcsr(self, order, dbgout) class(multigrid_t), target, intent(inout) :: self integer, intent(in) :: order integer, intent(in) :: dbgout integer :: lvl_f ! Level of fine mesh integer :: lvl_c ! Level of coarse mesh type(csrmat_t), pointer :: rcsr_fc ! Pointer to restriction matrix between lvl_f and lvl_c type(mesh_cart_t), pointer :: mmesh_f, mmesh_c real(FP), dimension(-order:order, -order:order) :: rcoeff integer :: n_nonzero, lc, lf, i, j, i_s, j_s, lfn integer, dimension(:,:,:), allocatable :: work_adjacents do lvl_f = 1, self%nlvls - 1 lvl_c = lvl_f + 1 if (dbgout >= 1) then write(get_stdout(), '(A,I3,A,I3)') & ' Creating restriction matrix, & &lvl_f = ', lvl_f, '; lvl_c = ', lvl_c endif rcsr_fc => self%mrcsr(lvl_f) mmesh_f => self%get_mesh_pointer(lvl_f) mmesh_c => self%get_mesh_pointer(lvl_c) allocate(work_adjacents(self%get_np(lvl_c), & -order:order, -order:order)) ! Coefficients for restriction matrix if (order == 0) then rcoeff(0, 0) = 1.0_FP elseif (order == 1) then rcoeff(-1, -1) = 1.0_FP / 16.0_FP rcoeff( 0, -1) = 1.0_FP / 8.0_FP rcoeff( 1, -1) = 1.0_FP / 16.0_FP rcoeff(-1, 0) = 1.0_FP / 8.0_FP rcoeff( 0, 0) = 1.0_FP / 4.0_FP rcoeff( 1, 0) = 1.0_FP / 8.0_FP rcoeff(-1, 1) = 1.0_FP / 16.0_FP rcoeff( 0, 1) = 1.0_FP / 8.0_FP rcoeff( 1, 1) = 1.0_FP / 16.0_FP else write(get_stderr(), *) & 'error(multigrid%create_mrcsr): order not valid', order error stop endif ! Parallel region for doing expensive work !$omp parallel default(none) private(lc, i_s, j_s, i, j, lf, lfn) & !$omp shared(self, lvl_c, lvl_f, mmesh_c, mmesh_f, order, & !$omp work_adjacents) !$omp do do lc = 1, self%get_np(lvl_c) ! Initialise all adjacent points to zero do i_s = -order, order do j_s = -order, order work_adjacents(lc, i_s, j_s) = 0 end do end do if (mmesh_c%is_inner_point(lc)) then ! Perform weighted 9 point average for inner points ! Restriction of boundary points and ghost points to zero i = mmesh_c%get_cart_i(lc) j = mmesh_c%get_cart_j(lc) lf = mmesh_f%cart_to_mesh_index(i, j) if (lf == 0) then write(get_stderr(), *) & 'error(multigrid%create_mrcsr): fine & &point could not be found for coarse point', lc error stop else do i_s = -order, order do j_s = -order, order lfn = mmesh_f%get_index_neighbor(i_s, j_s, lf) if (lfn == 0) then ! If point cannot be found for restriction, ! use fine point itself lfn = lf elseif (mmesh_f%is_ghost_point(lfn)) then ! Also, if the neighbour point is a ghost ! point, use the fine point lfn = lf endif work_adjacents(lc, i_s, j_s) = lfn enddo enddo endif endif enddo !$omp end do !$omp end parallel ! Allocate matrix rcsr_fc%ndim = self%get_np(lvl_c) rcsr_fc%ncol = self%get_np(lvl_f) rcsr_fc%nnz = self%get_np_inner(lvl_c) & * (2 * order + 1)**2 allocate(rcsr_fc%i(rcsr_fc%ndim + 1), source=0) allocate(rcsr_fc%j(rcsr_fc%nnz), source=0) allocate(rcsr_fc%val(rcsr_fc%nnz), source=0.0_FP) ! Serial region for filling matrix n_nonzero = 0 do lc = 1, self%get_np(lvl_c) rcsr_fc%i(lc) = n_nonzero + 1 do i_s = -order, order do j_s = -order, order lfn = work_adjacents(lc, i_s, j_s) if (lfn /= 0) then n_nonzero = n_nonzero + 1 rcsr_fc%j(n_nonzero) = lfn rcsr_fc%val(n_nonzero) = rcoeff(i_s, j_s) endif enddo enddo enddo rcsr_fc%i(rcsr_fc%ndim + 1) = n_nonzero + 1 rcsr_fc%nnz = n_nonzero deallocate(work_adjacents) if (dbgout >= 1) call rcsr_fc%display() end do end subroutine module subroutine create_mrcsr_inner(self, dbgout) class(multigrid_t), target, intent(inout) :: self integer, intent(in) :: dbgout integer :: lvl_f ! Level of fine mesh integer :: lvl_c ! Level of coarse mesh type(csrmat_t), pointer :: rcsr_fc, rcsr_inner_fc ! Pointer to restriction matrices between lvl_f and lvl_c integer, parameter :: max_nz = 9 ! Maximum number of non-zero elements in a row of rcsr_fc type(mesh_cart_t), pointer :: mmesh_f, mmesh_c integer, dimension(:,:), allocatable :: j_store real(FP), dimension(:,:), allocatable :: val_store integer, dimension(:), allocatable :: full_to_inner_f integer :: n_nonzero, kic, lc, s, s_shift, lf, kif, np real(FP) :: deficit, supplement do lvl_f = 1, self%nlvls - 1 lvl_c = lvl_f + 1 if (dbgout >= 1) then write(get_stdout(), '(A,I3,A,I3)') & ' Creating inner restriction matrix, & &lvl_f = ', lvl_f, '; lvl_c = ', lvl_c endif rcsr_fc => self%mrcsr(lvl_f) rcsr_inner_fc => self%mrcsr_inner(lvl_f) mmesh_f => self%get_mesh_pointer(lvl_f) mmesh_c => self%get_mesh_pointer(lvl_c) allocate(j_store(self%get_np_inner(lvl_c), max_nz)) allocate(val_store(self%get_np_inner(lvl_c), max_nz)) allocate(full_to_inner_f(self%get_np(lvl_f))) ! Allocate matrix rcsr_inner_fc%ndim = self%get_np_inner(lvl_c) rcsr_inner_fc%ncol = self%get_np_inner(lvl_f) rcsr_inner_fc%nnz = rcsr_fc%nnz allocate(rcsr_inner_fc%i(rcsr_inner_fc%ndim + 1), source=0) allocate(rcsr_inner_fc%j(rcsr_inner_fc%nnz), source=0) allocate(rcsr_inner_fc%val(rcsr_inner_fc%nnz), source=0.0_FP) ! Expensive, do in parallel !$omp parallel default(none) private(kic, lc, deficit, np, s, lf, & !$omp kif, supplement, s_shift) & !$omp shared(self, lvl_c, lvl_f, rcsr_inner_fc, rcsr_fc, mmesh_f, & !$omp mmesh_c, val_store, j_store, full_to_inner_f) !$omp do do lf = 1, self%get_np(lvl_f) full_to_inner_f(lf) = 0 enddo !$omp end do ! Requires thread synchronisation at this point (don't use nowait) !$omp do do kif = 1, self%get_np_inner(lvl_f) lf = mmesh_f%inner_indices(kif) full_to_inner_f(lf) = kif enddo !$omp end do ! Requires thread synchronisation at this point (don't use nowait) !$omp do do kic = 1, self%get_np_inner(lvl_c) do s = 1, max_nz ! Initialise the storage arrays to zero j_store(kic, s) = 0 val_store(kic, s) = 0.0_FP enddo lc = mmesh_c%inner_indices(kic) ! Compute deficit/supplement of restriction for points near ! boundary, where conventional restriction includes boundary ! points deficit = 1.0_FP np = 0 do s = rcsr_fc%i(lc), rcsr_fc%i(lc + 1) - 1 lf = rcsr_fc%j(s) ! Find inner index corresponding to lf kif = full_to_inner_f(lf) if (kif /= 0) then deficit = deficit - rcsr_fc%val(s) np = np + 1 endif enddo ! Compute supplement to restriction that will be homogeneously ! added to remaining points involved in restriction if (np /= 0) then supplement = deficit / (1.0_FP*np) else write(get_stderr(), *) & 'error(multigrid%create_mrcsr_inner): & &inner restriction matrix could not be built', lc error stop endif ! Build row of inner restriction matrix do s = rcsr_fc%i(lc), rcsr_fc%i(lc + 1) - 1 lf = rcsr_fc%j(s) ! Find inner index corresponding to lf kif = full_to_inner_f(lf) s_shift = s - rcsr_fc%i(lc) + 1 if (s_shift > max_nz) then write(get_stderr(), *) & "error(multigrid%create_mrcsr_inner): & &s_shift larger than allocated array. & &Increase max_nz" error stop endif if (kif /= 0) then j_store(kic, s_shift) = kif val_store(kic, s_shift) = rcsr_fc%val(s) + supplement endif enddo enddo !$omp end do !$omp end parallel ! Serial, fill CSR matrix n_nonzero = 0 do kic = 1, self%get_np_inner(lvl_c) rcsr_inner_fc%i(kic) = n_nonzero + 1 lc = mmesh_c%inner_indices(kic) do s = rcsr_fc%i(lc), rcsr_fc%i(lc + 1) - 1 s_shift = s - rcsr_fc%i(lc) + 1 kif = j_store(kic, s_shift) if (kif /= 0) then n_nonzero = n_nonzero + 1 rcsr_inner_fc%j(n_nonzero) = kif rcsr_inner_fc%val(n_nonzero) = val_store(kic, s_shift) endif enddo enddo rcsr_inner_fc%i(rcsr_inner_fc%ndim + 1) = n_nonzero + 1 rcsr_inner_fc%nnz = n_nonzero deallocate(j_store) deallocate(val_store) deallocate(full_to_inner_f) if (dbgout >= 1) call rcsr_inner_fc%display() end do end subroutine module subroutine create_mpcsr(self, dbgout) class(multigrid_t), target, intent(inout) :: self integer, intent(in) :: dbgout integer :: lvl_f ! Level of fine mesh integer :: lvl_c ! Level of coarse mesh type(csrmat_t), pointer :: pcsr_fc ! Pointer to prolongation matrix between lvl_f and lvl_c type(mesh_cart_t), pointer :: mmesh_f, mmesh_c real(FP), dimension(-1:1, -1:1) :: pcoeff integer :: n_nonzero, lf, i, j, i_s, j_s, i_c, j_c, lcn integer :: lvstf integer, dimension(:,:,:), allocatable :: work_adjacents do lvl_f = 1, self%nlvls - 1 lvl_c = lvl_f + 1 if (dbgout >= 1) then write(get_stdout(), '(A,I3,A,I3)') & ' Creating prolongation matrix, & &lvl_f = ', lvl_f, '; lvl_c = ', lvl_c endif pcsr_fc => self%mpcsr(lvl_f) mmesh_f => self%get_mesh_pointer(lvl_f) mmesh_c => self%get_mesh_pointer(lvl_c) allocate(work_adjacents(self%get_np(lvl_f), & -1:1, -1:1)) ! Coefficients for prolongation matrix pcoeff(-1, -1) = 1.0_FP / 4.0_FP pcoeff( 0, -1) = 1.0_FP / 2.0_FP pcoeff( 1, -1) = 1.0_FP / 4.0_FP pcoeff(-1, 0) = 1.0_FP / 2.0_FP pcoeff( 0, 0) = 1.0_FP pcoeff( 1, 0) = 1.0_FP / 2.0_FP pcoeff(-1, 1) = 1.0_FP / 4.0_FP pcoeff( 0, 1) = 1.0_FP / 2.0_FP pcoeff( 1, 1) = 1.0_FP / 4.0_FP lvstf = mmesh_f%get_lvst() ! Parallel region for expensive part !$omp parallel do schedule(static) default(none) & !$omp private(lf, i, j, i_s, j_s, i_c, j_c) & !$omp shared(self, work_adjacents, lvl_f, lvl_c, & !$omp mmesh_f, mmesh_c, lvstf) do lf = 1, self%get_np(lvl_f) i = mmesh_f%get_cart_i(lf) j = mmesh_f%get_cart_j(lf) do i_s = -1, 1 do j_s = -1, 1 i_c = i + i_s * lvstf j_c = j + j_s * lvstf work_adjacents(lf, i_s, j_s) = & mmesh_c%cart_to_mesh_index(i_c, j_c) enddo enddo enddo !$omp end parallel do ! Allocate matrix pcsr_fc%ndim = self%get_np(lvl_f) pcsr_fc%ncol = self%get_np(lvl_c) pcsr_fc%nnz = count(work_adjacents /= 0) allocate(pcsr_fc%i(pcsr_fc%ndim + 1), source=0) allocate(pcsr_fc%j(pcsr_fc%nnz), source=0) allocate(pcsr_fc%val(pcsr_fc%nnz), source=0.0_FP) ! Serial part for CSR matrix generation n_nonzero = 0 do lf = 1, self%get_np(lvl_f) pcsr_fc%i(lf) = n_nonzero + 1 do i_s = -1, 1 do j_s = -1, 1 lcn = work_adjacents(lf, i_s, j_s) if (lcn /= 0) then n_nonzero = n_nonzero + 1 pcsr_fc%j(n_nonzero) = lcn pcsr_fc%val(n_nonzero) = pcoeff(i_s, j_s) endif enddo enddo enddo pcsr_fc%i(pcsr_fc%ndim + 1) = n_nonzero + 1 pcsr_fc%nnz = n_nonzero deallocate(work_adjacents) if (dbgout >= 1) call pcsr_fc%display() end do end subroutine end submodule