module PDEelliptic !! Various tools for solving elliptic partial differential equations - uses MUMPS, scalapack, lapack, openmpi, and blas use, intrinsic:: iso_fortran_env, only: stderr=>error_unit, stdout=>output_unit use mumps_interface, only : mumps_struc, mumps_exec use mpimod, only: mpi_comm_world use phys_consts, only: wp, debug implicit none private public :: elliptic3D_cart,elliptic2D_cart,elliptic2D_polarization,elliptic2D_polarization_periodic,& elliptic_workers, check_mumps_status, quiet_mumps interface ! elliptic2d.f90 module function elliptic2D_polarization(srcterm,SigP2,SigP3,SigH,gradSigH2,gradSigH3,Cm,v2,v3,Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,dx1, & dx1i,dx2all,dx2iall,dx3all,dx3iall,Phi0,perflag,it) real(wp), dimension(:,:), intent(in) :: srcterm,SigP2,SigP3,SigH,Cm,gradSigH2,gradSigH3,v2,v3 !! ZZZ - THESE WILL NEED TO BE MODIFIED CONDUCTIVITIES, AND WE'LL NEED THREE OF THEM real(wp), dimension(:), intent(in) :: Vminx2,Vmaxx2 real(wp), dimension(:), intent(in) :: Vminx3,Vmaxx3 real(wp), intent(in) :: dt real(wp), dimension(0:), intent(in) :: dx1 !backward diffs start at index zero due to ghost cells real(wp), dimension(:), intent(in) :: dx1i !centered diffs do not include any ghost cells real(wp), dimension(0:), intent(in) :: dx2all real(wp), dimension(:), intent(in) :: dx2iall real(wp), dimension(0:), intent(in) :: dx3all real(wp), dimension(:), intent(in) :: dx3iall real(wp), dimension(:,:), intent(in) :: Phi0 logical, intent(in) :: perflag integer, intent(in) :: it real(wp), dimension(size(SigP2,1),size(SigP2,2)) :: elliptic2D_polarization end function elliptic2D_polarization module function elliptic2D_polarization_periodic(srcterm,SigP,SigH,gradSigH2,gradSigH3,Cm,v2,v3,Vminx2,Vmaxx2, & Vminx3,Vmaxx3,dt,dx1,dx1i,dx2all,dx2iall,dx3all,dx3iall,Phi0,perflag,it) real(wp), dimension(:,:), intent(in) :: srcterm,SigP,SigH,gradSigH2,gradSigH3,Cm,v2,v3 real(wp), dimension(:), intent(in) :: Vminx2,Vmaxx2 real(wp), dimension(:), intent(in) :: Vminx3,Vmaxx3 real(wp), intent(in) :: dt real(wp), dimension(0:), intent(in) :: dx1 !backward diffs start at index zero due to ghost cells real(wp), dimension(:), intent(in) :: dx1i !centered diffs do not include any ghost cells real(wp), dimension(0:), intent(in) :: dx2all real(wp), dimension(:), intent(in) :: dx2iall real(wp), dimension(0:), intent(in) :: dx3all real(wp), dimension(:), intent(in) :: dx3iall real(wp), dimension(:,:), intent(in) :: Phi0 logical, intent(in) :: perflag integer, intent(in) :: it real(wp), dimension(size(SigP,1),size(SigP,2)) :: elliptic2D_polarization_periodic end function elliptic2D_polarization_periodic module function elliptic2D_cart(srcterm,sig0,sigP,Vminx1,Vmaxx1,Vminx3,Vmaxx3,& dx1,dx1i,dx3all,dx3iall,flagdirich,perflag,gridflag,it) real(wp), dimension(:,:,:), intent(in) :: srcterm,sig0,sigP !arrays passed in will still have full rank 3 real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1 real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3 real(wp), dimension(0:), intent(in) :: dx1 !backward diffs start at index zero due to ghost cells real(wp), dimension(:), intent(in) :: dx1i !centered diffs do not include any ghost cells real(wp), dimension(0:), intent(in) :: dx3all real(wp), dimension(:), intent(in) :: dx3iall integer, intent(in) :: flagdirich logical, intent(in) :: perflag integer, intent(in) :: gridflag integer, intent(in) :: it real(wp), dimension(size(sig0,1),1,size(sig0,3)) :: elliptic2D_cart end function elliptic2D_cart end interface interface ! elliptic3d module function elliptic3D_cart(srcterm,Ac,Bc,Cc,Dc,Ec,Fc,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, & dx1,dx1i,dx2all,dx2iall,dx3all,dx3iall,flagdirich,perflag,it) real(wp), dimension(:,:,:), intent(in) :: srcterm,Ac,Bc,Cc,Dc,Ec,Fc real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1 real(wp), dimension(:,:), intent(in) :: Vminx2,Vmaxx2 real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3 real(wp), dimension(0:), intent(in) :: dx1 !backweard diffs start at index zero due to ghost cells real(wp), dimension(:), intent(in) :: dx1i !centered diffs do not include any ghost cells real(wp), dimension(0:), intent(in) :: dx2all real(wp), dimension(:), intent(in) :: dx2iall real(wp), dimension(0:), intent(in) :: dx3all real(wp), dimension(:), intent(in) :: dx3iall integer, intent(in) :: flagdirich logical, intent(in) :: perflag integer, intent(in) :: it real(wp), dimension(size(srcterm,1),size(srcterm,2),size(srcterm,3)) :: elliptic3D_cart end function elliptic3D_cart end interface integer, dimension(:), pointer, protected, save :: mumps_perm !cached permutation, unclear whether save is necessary... contains subroutine quiet_mumps(obj) !! this must be called AFTER the first mumps call that had job=-1 !! Needs Mumps >= 5.2 to actually take effect, does nothing on older MUMPS !! it stops the 100's of megabytes of logging console text, probably speeding up as well type(MUMPS_STRUC), intent(inout) :: obj obj%icntl(1) = stderr ! error messages obj%icntl(2) = stdout ! diagnosic, statistics, and warning messages obj%icntl(3) = stdout! ! global info, for the host (myid==0) obj%icntl(4) = 1 ! default is 2, this reduces verbosity end subroutine quiet_mumps subroutine elliptic_workers() !! ALLOWS WORKERS TO ENTER MUMPS SOLVES type(MUMPS_STRUC) :: mumps_par !FIRE UP MUMPS mumps_par%COMM = MPI_COMM_WORLD mumps_par%JOB = -1 mumps_par%SYM = 0 mumps_par%PAR = 1 call MUMPS_exec(mumps_par) call quiet_mumps(mumps_par) !ROOT WILL LOAD OUR PROBLEM !SOLVE (ALL WORKERS NEED TO SEE THIS CALL) mumps_par%JOB = 6 call MUMPS_exec(mumps_par) call check_mumps_status(mumps_par, 'elliptic_workers') !DEALLOCATE STRUCTURES USED BY WORKERS DURING SOLVE mumps_par%JOB = -2 call MUMPS_exec(mumps_par) end subroutine elliptic_workers subroutine check_mumps_status(p, name) !! check if Mumps error occurred type(MUMPS_STRUC), intent(in) :: p character(*), intent(in) :: name if (p%info(1) < 0 .or. p%infog(1) < 0) then write(stderr, *) 'Gemini:PDEelliptic:' // name // ' MUMPS ERROR: details:' if (p%info(1) == -1) write(stderr,'(a,i4)') 'the error was reported by processor #',p%info(2) write(stderr, *) 'Mumps Error: info(1,2,8):', p%info(1:2), p%info(8) write(stderr, *) 'Mumps Error: infoG(1,2)', p%infoG(1:2) write(stderr, *) 'for error number meaning, see "8 Error Diagnostics" of MUMPS manual' error stop endif end subroutine check_mumps_status end module PDEelliptic