module solver3d_m !! 3D solvers, involving MPI communication use MPI use screen_io_m, only: get_stdout, get_stderr use comm_handling_m, only: is_master use precision_m, only : FP, MPI_FP, FP_NAN use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_SOLVER3D implicit none abstract interface subroutine matvec_interface(u, v, ipar) !! Interface for matrix vector multiplication v = A*u import FP real(FP), dimension(*) :: u real(FP), dimension(*) :: v integer :: ipar(*) end subroutine end interface abstract interface subroutine precond_interface(u, v, ipar) !! Interface for preconditioner import FP real(FP), dimension(*) :: u real(FP), dimension(*) :: v integer :: ipar(*) end subroutine end interface type, public, abstract :: solver3d_t !! Abstract class for 3D solver integer, private :: ndim_loc !! Local dimension of problem procedure(matvec_interface), pointer, nopass, private :: matvec => err_no_matvec !! Function pointer to matrix-vector multiplication routine procedure(precond_interface), pointer, nopass, private :: precondl => err_no_prec !! Function pointer to left preconditioning routine procedure(precond_interface), pointer, nopass, private :: precondr => err_no_prec !! Function pointer to right preconditioning routine contains procedure(create_abstract), deferred :: create procedure(solve_abstract), deferred :: solve end type private :: err_no_matvec private :: err_no_prec interface module subroutine create_abstract(self, comm, krylov_method, & ndim_loc, resmax, & matvec, precondl, precondr, & maxiter, nrestart, & dbgout) !! Initialisation of 3d solver class(solver3d_t), intent(inout) :: self !! Instance of class integer, intent(in) :: comm !! MPI communicator character(len=*), intent(in) :: krylov_method !! Krylov method to be used integer, intent(in) :: ndim_loc !! Local dimension (on rank) of matrix real(FP), intent(in) :: resmax !! Maximum allowed residual of solution procedure(matvec_interface) :: matvec !! Routine performing matrix vector multiplication procedure(precond_interface), optional :: precondl !! Routine performing left preconditioning procedure(precond_interface), optional :: precondr !! Routine performing right preconditioning integer, intent(in), optional :: maxiter !! Maximum number of outer iterations (default=10) integer, intent(in), optional :: nrestart !! Number of iterations before restart (default=10) integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine solve_abstract(self, comm, rhs, sol, res, info, & res_true) !! Solves 3d problem class(solver3d_t), intent(inout) :: self !! Instance of class integer, intent(in) :: comm !! MPI communicator real(FP), dimension(self%ndim_loc), intent(in) :: rhs !! Right hand side real(FP), dimension(self%ndim_loc), intent(inout) :: sol !! Solution on output, on input initial guess real(FP), intent(out) :: res !! Relative pseudo-residuum after solve integer, intent(out) :: info !! Info parameter, successful solve if >= 0 real(FP), intent(out), optional :: res_true !! True relative residuum after solve end subroutine end interface type, public, extends(solver3d_t) :: solver3d_PIM_t !! Iterative 3D solver based on PIM library integer, private :: outi !! Debug output level character(len=:), allocatable, private :: krylov_method real(FP), private :: resmax !! Maximum allowed tolerance of solution x on normalised residual integer, public, dimension(13) :: ipar !! Integer parameters of PIM solver real(FP), private, dimension(6) :: dpar !! Real parameters of PIM solver contains procedure, public :: create => create_PIM procedure, public :: solve => solve_PIM final :: destructor_PIM end type interface module subroutine create_PIM(self, comm, krylov_method, & ndim_loc, resmax, & matvec, precondl, precondr, & maxiter, nrestart, & dbgout) !! Initialisation of 3d solver class(solver3d_PIM_t), intent(inout) :: self !! PIM solver integer, intent(in) :: comm !! MPI communicator character(len=*), intent(in) :: krylov_method !! Krylov method to be used, available options are (see PIM manual): !! - CG !! - RGMRES !! - BICGSTAB !! - RBICGSTAB integer, intent(in) :: ndim_loc !! Local dimension (on rank) of matrix real(FP), intent(in) :: resmax !! Maximum allowed tolerance of solution x on normalised residual procedure(matvec_interface) :: matvec !! Routine performing matrix vector multiplication procedure(precond_interface), optional :: precondl !! Routine performing left preconditioning procedure(precond_interface), optional :: precondr !! Routine performing right preconditioning integer, intent(in), optional :: maxiter !! Maximum number of outer iterations (default=10) integer, intent(in), optional :: nrestart !! Number of iterations before restart (default=10) integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine solve_PIM(self, comm, rhs, sol, res, info, res_true) !! Solves 3d problem class(solver3d_PIM_t), intent(inout) :: self !! PIM solver integer, intent(in) :: comm !! Not used, !! only the communicator from the create_PIM routine is used real(FP), dimension(self%ndim_loc), intent(in) :: rhs !! Right hand side real(FP), dimension(self%ndim_loc), intent(inout) :: sol !! Solution on output, on input initial guess real(FP), intent(out) :: res !! Relative pseudo-residuum after solve !! see PIM manual table 2 integer, intent(out) :: info !! On succesful solve (info>0) the number of iterations/restarts, !! for negative numbers see PIM manual real(FP), intent(out), optional :: res_true !! True relative residuum after solve end subroutine module subroutine destructor_PIM(self) !! Clears solver type(solver3d_PIM_t), intent(inout) :: self !! PIM solver end subroutine end interface #ifdef ENABLE_CERFACS type, public, extends(solver3d_t) :: solver3d_CERFACS_t !! Iterative 3D solver based on CERFACS library integer, private :: outi !! Debug output level integer, dimension(8), private :: icntl !! Integer solver parameters (see CERFACS manual) real(FP), dimension(5), private :: cntl !! Real solver parameters (see CERFACS manual) integer, private :: ndim_glob !! Global dimension of problem integer, private :: nrestart !! Restart parameter integer, private :: lwork !! Dimension of work array contains procedure, public :: create => create_CERFACS procedure, public :: solve => solve_CERFACS final :: destructor_CERFACS end type interface module subroutine create_CERFACS(self, comm, krylov_method, & ndim_loc, resmax, & matvec, precondl, precondr, & maxiter, nrestart, & dbgout) !! Initialisation of 3d solver class(solver3d_CERFACS_t), intent(inout) :: self !! CERFACS solver integer, intent(in) :: comm !! MPI communicator character(len=*), intent(in) :: krylov_method !! Krylov method to be used, available options are: !! - RGMRES integer, intent(in) :: ndim_loc !! Local dimension (on rank) of matrix real(FP), intent(in) :: resmax !! Maximum allowed tolerance of solution x on normalised residual procedure(matvec_interface) :: matvec !! Routine performing matrix vector multiplication procedure(precond_interface), optional :: precondl !! Routine performing left preconditioning procedure(precond_interface), optional :: precondr !! Routine performing right preconditioning integer, intent(in), optional :: maxiter !! TODO integer, intent(in), optional :: nrestart !! TODO integer, intent(in), optional :: dbgout !! Debug output level end subroutine module subroutine solve_CERFACS(self, comm, rhs, sol, res, info, & res_true) !! Solves 3d problem class(solver3d_CERFACS_t), intent(inout) :: self !! CERFACS solver integer, intent(in) :: comm !! TODO real(FP), dimension(self%ndim_loc), intent(in) :: rhs !! Right hand side real(FP), dimension(self%ndim_loc), intent(inout) :: sol !! Solution on output, on input initial guess real(FP), intent(out) :: res !! Residuum of preconditioned system integer, intent(out) :: info !! On succesful solve (info>0) the number of iterations real(FP), intent(out), optional :: res_true !! Residuum of unpreconditioned system end subroutine module subroutine destructor_CERFACS(self) !! Clears solver type(solver3d_CERFACS_t), intent(inout) :: self !! PIM solver end subroutine end interface #endif contains subroutine err_no_matvec(u, v, ipar) !! Error routine, called if matvec routine is not correctly associated real(FP), dimension(*) :: u real(FP), dimension(*) :: v integer :: ipar(*) call handle_error('matvec routine not correctly associated', & PARALLAX_ERR_SOLVER3D, __LINE__, __FILE__) end subroutine subroutine err_no_prec(u, v, ipar) !! Error routine, called if precondl or precondr routine !! is not correctly associated real(FP), dimension(*) :: u real(FP), dimension(*) :: v integer :: ipar(*) call handle_error('precond routine not correctly associated', & PARALLAX_ERR_SOLVER3D, __LINE__, __FILE__) end subroutine end module