submodule(solver_aligned3d_m) solve_aligned3d_s !! Contains routine that performs a 3D aligned solve for canonical elliptic problem implicit none contains module subroutine solve_aligned3d(comm_handler, & equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & map, & params_solver, solver_type, & rhs, sol, & niter, res, true_res, cnt_matvec, & lambda, xi, co_stag, fcx, & bnd_descrs, & flutter_fac_cano, flutter_fac_stag, & opsinplane_cano, opsinplane_stag, & apar_stag, apar_cano, & bnd_descrs_flux) type(comm_handler_t), intent(in) :: comm_handler type(mesh_cart_t), intent(in) :: mesh_cano type(mesh_cart_t), intent(in) :: mesh_stag type(equilibrium_storage_t), intent(in):: equi_on_cano type(equilibrium_storage_t), intent(in) :: equi_on_stag type(parallel_map_t), intent(in) :: map type(params_solver_aligned3d_t), intent(in) :: params_solver character(len=*), intent(in) :: solver_type real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: rhs real(GP), dimension(mesh_cano%get_n_points()), intent(inout) :: sol integer, intent(out) :: niter real(GP), intent(out) :: res real(GP), intent(out) :: true_res integer, intent(out), optional :: cnt_matvec real(GP), dimension(mesh_cano%get_n_points()), intent(in), optional :: lambda real(GP), dimension(mesh_cano%get_n_points()), intent(in), optional :: xi real(GP), dimension(mesh_stag%get_n_points()), intent(in), optional :: co_stag real(GP), dimension(mesh_cano%get_n_points()), intent(in), optional :: fcx integer, dimension(mesh_cano%get_n_points_boundary()), intent(in), optional :: bnd_descrs real(GP), intent(in), optional :: flutter_fac_cano real(GP), intent(in), optional :: flutter_fac_stag type(inplane_operators_t), intent(in), optional :: opsinplane_cano type(inplane_operators_t), intent(in), optional :: opsinplane_stag real(GP), dimension(mesh_stag%get_n_points()), intent(in), optional :: apar_stag real(GP), dimension(mesh_cano%get_n_points()), intent(in), optional :: apar_cano integer, dimension(mesh_stag%get_n_points_boundary()), intent(in), optional :: bnd_descrs_flux logical :: with_flutter, any_flutter_input integer :: np, np_stag, np_cano_fwd, np_stag_bwd, nrecv integer, allocatable, dimension(:) :: np_recv integer :: cnt_matvec_eff class(solver3d_t), allocatable :: solver3d real(GP), allocatable, dimension(:) :: ufc_hfwd, pgrad_hbwd ! Work arrays within matvec routine below real(GP), allocatable, dimension(:) :: diag_inv ! Inverse diagonal used for Jacobi preconditioning ! Check consistency of input with_flutter = present(flutter_fac_cano) & .and. present(flutter_fac_stag) & .and. present(opsinplane_cano) & .and. present(opsinplane_stag) & .and. present(apar_stag) & .and. present(apar_cano) any_flutter_input = present(flutter_fac_cano) & .or. present(flutter_fac_stag) & .or. present(opsinplane_cano) & .or. present(opsinplane_stag) & .or. present(apar_stag) & .or. present(apar_cano) if (.not.with_flutter .and. any_flutter_input) then call handle_error('Input of flutter quantities not consistent', & GRILLIX_ERR_SOLVER3D, __LINE__, __FILE__) endif ! Switch off flutter computation if flutter_facs=0 if (with_flutter) then if ( (abs(flutter_fac_cano) < GP_EPS) .and. (abs(flutter_fac_stag) < GP_EPS) ) then with_flutter = .false. endif endif ! Number of points of meshes on adjacent planes np = mesh_cano%get_n_points() np_stag = mesh_stag%get_n_points() call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, 1, [np_stag], nrecv, np_recv) np_stag_bwd = np_recv(1) deallocate(np_recv) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, 1, [np], nrecv, np_recv) np_cano_fwd = np_recv(1) deallocate(np_recv) ! Allocate work variables used within matvec operation allocate(ufc_hfwd(np_cano_fwd)) allocate(pgrad_hbwd(np_stag_bwd)) ! Compute inverse diagonal for Jacobi preconditioning if (params_solver%precond_type == 'JAC') then call compute_diag_inv_approx() endif ! Initialise 3D solver call solver3d_factory(solver_type, solver3d) select case(params_solver%precond_type) case('NONE') call solver3d%create(comm_handler%get_comm_planes(), & krylov_method=params_solver%krylov_method, & ndim_loc=mesh_cano%get_n_points(), & resmax=params_solver%resmax, & matvec=matvec, & maxiter=params_solver%maxiter, & nrestart=params_solver%nrestart, & dbgout=params_solver%dbgout) case('JAC') call solver3d%create(comm_handler%get_comm_planes(), & krylov_method=params_solver%krylov_method, & ndim_loc=mesh_cano%get_n_points(), & resmax=params_solver%resmax, & matvec=matvec, & precondl=precond_jacobi, & maxiter=params_solver%maxiter, & nrestart=params_solver%nrestart, & dbgout=params_solver%dbgout) case default call handle_error('Preconditioner not valid', & GRILLIX_ERR_SOLVER3D, __LINE__, __FILE__, & error_info_t(params_solver%precond_type)) end select ! Solve system cnt_matvec_eff = 0 call solver3d%solve(comm_handler%get_comm_planes(), & rhs, sol, res, niter, true_res) if (present(cnt_matvec)) then cnt_matvec = cnt_matvec_eff endif deallocate(solver3d) deallocate(pgrad_hbwd) deallocate(ufc_hfwd) if (allocated(diag_inv)) deallocate(diag_inv) contains subroutine matvec(u, v, ipar) !! Matrix Vector Multiplication for parallel helmholtz operator real(GP), dimension(*) :: u real(GP), dimension(*) :: v integer :: ipar(*) integer :: l, ki, kb, kg, jb, bndnm_nnz, request(2) integer, dimension(3) :: bndnmn_cols real(GP), dimension(3) :: bndnmn_vals real(GP) :: smb real(GP), dimension(mesh_cano%get_n_points()) :: ufc, pgrad_cano real(GP), dimension(mesh_stag%get_n_points()) :: pgrad, ufc_stag cnt_matvec_eff = cnt_matvec_eff + 1 if (.not.present(fcx)) then !$omp parallel default(none) & !$omp shared(mesh_cano, u, ufc) & !$omp private(l) !$omp do do l = 1, mesh_cano%get_n_points() ufc(l) = u(l) enddo !$omp end do !$omp end parallel else !$omp parallel default(none) & !$omp shared(mesh_cano, u, ufc, fcx) & !$omp private(l) !$omp do do l = 1, mesh_cano%get_n_points() ufc(l) = fcx(l) * u(l) enddo !$omp end do !$omp end parallel endif ! MPI communication with the forward plane for the parallel gradient call start_comm_planes(comm_handler%get_comm_planes(), & mesh_cano%get_n_points(), np_cano_fwd, ufc, +1, ufc_hfwd, request) call finalize_comm_planes(request) !$omp parallel default(none) & !$omp shared(mesh_stag, map, equi_on_stag, opsinplane_stag, & !$omp ufc, ufc_hfwd, ufc_stag, pgrad, co_stag, & !$omp with_flutter, flutter_fac_stag, apar_stag, bnd_descrs_flux) & !$omp private(l, ki, kb, kg, & !$omp bndnmn_cols, bndnmn_vals, bndnm_nnz, smb, jb) !$omp do do l = 1, mesh_stag%get_n_points() pgrad(l) = map%par_grad_single(ufc, ufc_hfwd, l) enddo !$omp end do if (with_flutter) then call map%lift_cano_to_stag(ufc, ufc_hfwd, ufc_stag) !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) pgrad(l) = pgrad(l) + flutter_fac_stag & * opsinplane_stag%flutter_grad(apar_stag, ufc_stag, l) enddo !$omp end do if (present(bnd_descrs_flux)) then ! Calculate boundary points using boundary condition ! Linear operation acting as a matrix-vector multiplication call set_perpbnds(mesh_stag, bnd_descrs_flux, pgrad) else ! Calculate boundary points using ghost points !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) pgrad(l) = pgrad(l) + flutter_fac_stag & * opsinplane_stag%flutter_grad(apar_stag, ufc_stag, l) enddo !$omp end do end if ! Extrapolation on ghost points ! Linear operation acting as a matrix-vector multiplication call extrapolate_ghost_points(mesh_stag, pgrad) endif if (present(co_stag)) then !$omp do do l = 1, mesh_stag%get_n_points() pgrad(l) = pgrad(l) * co_stag(l) enddo !$omp end do endif !$omp end parallel ! MPI communication with the backward plane for the parallel divergence call start_comm_planes(comm_handler%get_comm_planes(), & mesh_stag%get_n_points(), np_stag_bwd, pgrad, -1, pgrad_hbwd, request) call finalize_comm_planes(request) !$omp parallel default(none) & !$omp shared(mesh_cano, map, equi_on_cano, opsinplane_cano, & !$omp v, ufc, pgrad, pgrad_hbwd, pgrad_cano, xi, lambda, & !$omp with_flutter, flutter_fac_cano, apar_cano, bnd_descrs) & !$omp private(l, ki, kb, kg, & !$omp bndnmn_cols, bndnmn_vals, bndnm_nnz, smb, jb) !$omp do do l = 1, mesh_cano%get_n_points() v(l) = (-1.0_GP) * map%par_divb_single(pgrad, pgrad_hbwd, l) enddo !$omp end do if (with_flutter) then call map%flat_stag_to_cano(pgrad, pgrad_hbwd, pgrad_cano) ! Calculate flutter divergence on inner points !$omp do do ki = 1, mesh_cano%get_n_points_inner() l = mesh_cano%inner_indices(ki) v(l) = v(l) - flutter_fac_cano & * opsinplane_cano%flutter_div(apar_cano, pgrad_cano, l) enddo !$omp end do endif !$omp do do l = 1, mesh_cano%get_n_points() if (present(xi)) then v(l) = v(l) * xi(l) endif if (present(lambda)) then v(l) = lambda(l) * u(l) + v(l) endif enddo !$omp end do ! Adaption for boundary conditions if (present(bnd_descrs)) then !$omp do do kb = 1, mesh_cano%get_n_points_boundary() l = mesh_cano%boundary_indices(kb) if ( bnd_descrs(kb) == BND_TYPE_DIRICHLET_ZERO ) then v(l) = u(l) elseif ( bnd_descrs(kb) == BND_TYPE_NEUMANN ) then call compute_bndnmn_matrix_row(mesh_cano, l, bndnmn_cols, bndnmn_vals, bndnm_nnz) smb = 0.0_GP do jb = 1, bndnm_nnz smb = smb + bndnmn_vals(jb)*u(bndnmn_cols(jb)) enddo v(l) = smb else call handle_error('Matvec, boundary_descriptor not valid', & GRILLIX_ERR_SOLVER3D, __LINE__, __FILE__, & error_info_t('bnd_descr: ',[kb, bnd_descrs(kb)])) endif enddo !$omp end do !$omp do do kg = 1, mesh_cano%get_n_points_ghost() ! On ghost points identity is returned l = mesh_cano%ghost_indices(kg) v(l) = u(l) enddo !$omp end do endif !$omp end parallel end subroutine subroutine precond_jacobi(u, v, ipar) !! Jacobi preconditioner, i.e. diagonal part of matrix !! with inversion of perpendicular neumann boundary conditions !! and ghost points real(GP), dimension(*) :: u real(GP), dimension(*) :: v integer :: ipar(*) integer :: l, kb, kg real(GP), dimension(mesh_cano%get_n_points_boundary()) :: bnd_vals !$omp parallel default(none) & !$omp shared(mesh_cano, u, v, diag_inv, bnd_descrs, bnd_vals) & !$omp private(l, kb, kg) !$omp do do l = 1, mesh_cano%get_n_points() v(l) = diag_inv(l) * u(l) enddo !$omp end do ! Adaption for boundary conditions if (present(bnd_descrs)) then !$omp do do kb = 1, mesh_cano%get_n_points_boundary() l = mesh_cano%boundary_indices(kb) bnd_vals(kb) = u(l) enddo !$omp end do call set_perpbnds(mesh_cano, bnd_descrs, v , bnd_vals) !$omp do do kg = 1, mesh_cano%get_n_points_ghost() l = mesh_cano%ghost_indices(kg) v(l) = u(l) enddo !$omp end do endif !$omp end parallel end subroutine subroutine compute_diag_inv_approx() !! Computes an approximate for the inverse diagonal of the parallel helmholtz matrix !! \f[ !! D = \lambda + 2 * \xi *fcx * co_cano / (dpar_cano_cano_fwd * dpar_cano_cano_bwd) !! \f] !! !! Flutter terms are not taken into account integer :: l, request(2) real(GP), allocatable, dimension(:) :: co_stag_hbwd allocate(diag_inv(mesh_cano%get_n_points())) if (present(co_stag)) then allocate(co_stag_hbwd(np_stag_bwd)) ! MPI communication for co_stag with the backward plane call start_comm_planes(comm_handler%get_comm_planes(), & mesh_stag%get_n_points(), np_stag_bwd, & co_stag, -1, co_stag_hbwd, request) call finalize_comm_planes(request) endif !$omp parallel default(none) & !$omp shared(mesh_cano, map, diag_inv, & !$omp fcx, xi, lambda, co_stag, co_stag_hbwd) & !$omp private(l) !$omp do do l = 1, mesh_cano%get_n_points() if (present(fcx)) then diag_inv(l) = fcx(l) else diag_inv(l) = 1.0_GP endif if (present(xi)) then diag_inv(l) = diag_inv(l) * xi(l) endif if (present(co_stag)) then diag_inv(l) = diag_inv(l) * map%flat_stag_to_cano_single(co_stag, co_stag_hbwd, l) endif diag_inv(l) = diag_inv(l) * 2.0_GP / (map%dpar_cano_cano_fwd(l) * map%dpar_cano_cano_bwd(l)) if (present(lambda)) then diag_inv(l) = diag_inv(l) + lambda(l) endif diag_inv(l) = 1.0_GP / diag_inv(l) enddo !$omp end do !$omp end parallel if (allocated(co_stag_hbwd)) deallocate(co_stag_hbwd) end subroutine end subroutine end submodule