C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/code_cg3d_petsc/cg3d_petsc.F,v 1.3 2016/07/05 15:56:42 dgoldberg Exp $ C $Name: $ #include "CPP_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE CG3D_PETSC( U cg3d_x, ! solution vector U cg3d_b, ! rhs I tolerance, U maxIter, I myIter, I myThid ) C /============================================================\ C | SUBROUTINE | C | o | C |============================================================| C | | C \============================================================/ IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "CG3D.h" #ifdef ALLOW_PETSC #include "finclude/petsc.h" #include "CG3D_PETSC.h" ! UNCOMMENT IF V3.0 !#include "finclude/petscvec.h" !#include "finclude/petscmat.h" !#include "finclude/petscksp.h" !#include "finclude/petscpc.h" #endif C === Global variables === C !INPUT/OUTPUT ARGUMENTS C cg_Uin, cg_Vin - input and output velocities C cg_Bu, cg_Bv - driving stress INTEGER myThid INTEGER maxiter INTEGER myIter _RL tolerance _RL cg3d_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL cg3d_b (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) C LOCAL VARIABLES INTEGER i, j, bi, bj, cg_halo, conv_flag INTEGER iter, col_iter, k _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0 _RL dot_p1_tile (nSx,nSy) _RL dot_p2_tile (nSx,nSy) CHARACTER*(MAX_LEN_MBUF) msgBuf ! LOGICAL create_mat, destroy_mat #ifdef ALLOW_PETSC INTEGER indices(Nr*(snx*nsx*sny*nsy)) INTEGER cg3d_dofs_cum_sum (0:nPx*nPy*MAX_CG3D_PETSC_CPUINVERT-1), & idx(1) _RL rhs_values(Nr*(snx*nsx*sny*nsy)) _RL solution_values(Nr*(snx*nsx*sny*nsy)) _RL mat_values (7,1), mat_val_return(1) INTEGER indices_col(7) INTEGER local_dofs, global_dofs, dof_index, dof_index_col INTEGER local_offset INTEGER iters INTEGER color, rank PC subpc KSP subksp(1) Vec rhs Vec solution PetscErrorCode ierr #ifdef ALLOW_USE_MPI integer mpiRC, mpiMyWid #endif #endif ! print *, "GOT HERE CG3D", myByLo(mythid), myByHi(myThid), ! & myBxLo(mythid), myBxHi(myThid) ! RETURN #ifdef ALLOW_PETSC IF ((myIter.eq.1) .OR. (.not.cg3d_petsc_reuse_mat)) THEN cg3d_need2create_mat = .TRUE. ELSE cg3d_need2create_mat = .FALSE. ENDIF IF ((myIter.eq.nTimeSteps) .OR. (.not.cg3d_petsc_reuse_mat)) THEN cg3d_need2destroy_mat = .TRUE. ELSE cg3d_need2destroy_mat = .FALSE. ENDIF ! Print *, "GOT HERE CG3D", local_dofs,global_dofs, ! & cg3d_dofs_cum_sum #ifdef ALLOW_USE_MPI CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC ) local_dofs = cg3d_dofs_process (mpiMyWid) global_dofs = 0 cg3d_dofs_cum_sum(0) = 0 DO i=0,nPx*nPy*MAX_CG3D_PETSC_CPUINVERT-1 ! IF (i.le.cg3d_petsc_cpuInVert) THEN global_dofs = global_dofs + cg3d_dofs_process (i) if (i.ge.1) then cg3d_dofs_cum_sum(i) = cg3d_dofs_cum_sum(i-1)+ & cg3d_dofs_process(i-1) endif ! ENDIF ENDDO local_offset = cg3d_dofs_cum_sum(mpimywid) #else local_dofs = cg3d_dofs_process (0) global_dofs = local_dofs local_offset = 0 #endif Print *, "GOT HERE CG3D", local_dofs,global_dofs, & cg3d_dofs_cum_sum !---------------------- call VecCreate(PETSC_COMM_WORLD, rhs, ierr) call VecSetSizes(rhs, local_dofs, global_dofs, ierr) call VecSetType(rhs, VECMPI, ierr) call VecCreate(PETSC_COMM_WORLD, solution, ierr) call VecSetSizes(solution, local_dofs, global_dofs, ierr) call VecSetType(solution, VECMPI, ierr) do i=1,local_dofs indices(i) = i-1 + local_offset end do do i=1,Nr*nSx*nSy*sNx*sNy rhs_values (i) = 0. _d 0 solution_values (i) = 0. _d 0 enddo ! gather rhs and initial guess values to populate petsc vectors DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx DO k=1,Nr color = INT(cg3d_petsc_color(i,j,k,bi,bj)) rank = cg3d_color_rank (color) dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj)) & - local_offset if (dof_index.ge.0 .and. rank.eq.mpimywid) THEN rhs_values(dof_index+1) = cg3d_b(i,j,k,bi,bj) solution_values(dof_index+1) = cg3d_x(i,j,k,bi,bj) endif ENDDO ENDDO ENDDO ENDDO ENDDO call VecSetValues(rhs, local_dofs, indices, rhs_values, & INSERT_VALUES, ierr) call VecAssemblyBegin(rhs, ierr) call VecAssemblyEnd(rhs, ierr) call VecSetValues(solution, local_dofs, indices, & solution_values, INSERT_VALUES, ierr) call VecAssemblyBegin(solution, ierr) call VecAssemblyEnd(solution, ierr) if (cg3d_need2create_mat) then CALL TIMER_START ('CG3D_PETSC_MATCREATE',myThid) ! IF USING v3.0 THEN ! call MatCreateMPIAIJ (PETSC_COMM_WORLD, call MatCreateAIJ (PETSC_COMM_WORLD, & local_dofs, local_dofs, & global_dofs, global_dofs, & 7, PETSC_NULL_INTEGER, & 7, PETSC_NULL_INTEGER, & matrix_petsc_cg3d, ierr) ! populate petsc matrix DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx DO k=1,Nr dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj)) ! & - local_offset color = INT(cg3d_petsc_color(i,j,k,bi,bj)) rank = cg3d_color_rank (color) IF (dof_index .ge. 0 .and. rank.eq.mpimywid) THEN DO col_iter=1,7 indices_col(col_iter) = 0 mat_values(col_iter,1) = 0. _d 0 ENDDO col_iter = 0 ! ASSUME THE STENCIL ! ac3d (i,j,k) --> 1 ! aw3d (i,j,k) --> 2 ! as3d (i,j,k) --> 3 ! av3d (i,j,k) --> 4 ! aw3d (i+1,j,k) --> 5 ! as3d (i,j+1,k) --> 6 ! av3d (i,j,k+1) --> 7 dof_index_col = cg3d_petsc_dofs(i,j,k,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = ac3d(i,j,k,bi,bj) indices_col(col_iter) = dof_index_col endif dof_index_col = cg3d_petsc_dofs(i-1,j,k,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = aw3d(i,j,k,bi,bj) indices_col(col_iter) = dof_index_col endif dof_index_col = cg3d_petsc_dofs(i,j-1,k,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = as3d(i,j,k,bi,bj) indices_col(col_iter) = dof_index_col endif if (k.gt.1) then dof_index_col = cg3d_petsc_dofs(i,j,k-1,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = av3d(i,j,k,bi,bj) indices_col(col_iter) = dof_index_col endif endif dof_index_col = cg3d_petsc_dofs(i+1,j,k,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = aw3d(i+1,j,k,bi,bj) indices_col(col_iter) = dof_index_col endif dof_index_col = cg3d_petsc_dofs(i,j+1,k,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = as3d(i,j+1,k,bi,bj) indices_col(col_iter) = dof_index_col endif if (k.lt.Nr) then dof_index_col = cg3d_petsc_dofs(i,j,k+1,bi,bj) if (dof_index_col.ge.0) then col_iter = col_iter + 1 mat_values(col_iter,1) = av3d(i,j,k+1,bi,bj) indices_col(col_iter) = dof_index_col endif endif call matSetValues (matrix_petsc_cg3d, 1, dof_index, col_iter, & indices_col, & mat_values,INSERT_VALUES,ierr) ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO call MatAssemblyBegin(matrix_petsc_cg3d,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(matrix_petsc_cg3d,MAT_FINAL_ASSEMBLY,ierr) call KSPCreate(PETSC_COMM_WORLD, ksp_cg3d, ierr) call KSPSetOperators(ksp_cg3d, & matrix_petsc_cg3d, matrix_petsc_cg3d, & DIFFERENT_NONZERO_PATTERN, ierr) IF (CG3D_PETSC_PRECOND_TYPE.eq.'MUMPS') then call KSPSetType(ksp_cg3d,KSPPREONLY,ierr) ELSE SELECT CASE (CG3D_PETSC_SOLVER_TYPE) CASE ('CG') PRINT *, "PETSC SOLVER: SELECTED CG" call KSPSetType(ksp_cg3d, KSPCG, ierr) CASE ('GMRES') PRINT *, "PETSC SOLVER: SELECTED GMRES" call KSPSetType(ksp_cg3d, KSPGMRES, ierr) CASE ('BICG') PRINT *, "PETSC SOLVER: SELECTED BICG" call KSPSetType(ksp_cg3d, KSPBICG, ierr) CASE DEFAULT PRINT *, "PETSC SOLVER: SELECTED DEFAULT" call KSPSetType(ksp_cg3d, KSPCG, ierr) END SELECT ENDIF call KSPGetPC(ksp_cg3d, pc_cg3d, ierr) call KSPSetTolerances(ksp_cg3d,tolerance, & PETSC_DEFAULT_DOUBLE_PRECISION, & PETSC_DEFAULT_DOUBLE_PRECISION, & maxiter,ierr) SELECT CASE (CG3D_PETSC_PRECOND_TYPE) CASE ('BLOCKJACOBI') PRINT *, "PETSC PRECOND: SELECTED BJACOBI" call PCSetType(pc_cg3d, PCBJACOBI, ierr) call kspsetup (ksp_cg3d, ierr) call PCBJacobiGetSubKSP (pc_cg3d,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER, & subksp,ierr); call KSPGetPC (subksp, subpc, ierr) call PCSetType (subpc, PCILU, ierr) ! call PCFactorSetLevels(subpc,0,ierr) CASE ('JACOBI') PRINT *, "PETSC PRECOND: SELECTED JACOBI" call PCSetType(pc_cg3d, PCJACOBI, ierr) CASE ('ILU') PRINT *, "PETSC PRECOND: SELECTED ILU" call PCSetType(pc_cg3d, PCILU, ierr) CASE ('GAMG') PRINT *, "PETSC PRECOND: SELECTED GAMG" call PCSetType(pc_cg3d, PCGAMG, ierr) call PCGAMGSetCoarseEqLim(pc_cg3d,10000,ierr) call PCGAMGSetSymGraph(pc_cg3d, PETSC_TRUE,ierr) call PCGAMGSetNSmooths(pc_cg3d, 0,ierr) call PCGAMGSetThreshold(pc_cg3d, .001,ierr) ! call PCGAMGSetReuseProl(pc,PETSC_FALSE,ierr) call kspsetup (ksp_cg3d, ierr) CASE ('MUMPS') PRINT *, "PETSC PRECOND: SELECTED MUMPS" call PCSetType(pc_cg3d,PCLU,ierr) call PCFactorSetMatSolverPackage(pc_cg3d,MATSOLVERMUMPS,ierr) call PCFactorSetUpMatSolverPackage(pc_cg3d,ierr) call PCFactorGetMatrix(pc_cg3d,mumpsFac_cg3d,ierr) call MatMumpsSetIcntl(mumpsfac_cg3d,7,2,ierr) call MatMumpsSetIcntl(mumpsfac_cg3d,24,1,ierr) ! call MatMumpsSetCntl(mumpsfac,3,2.e-6,ierr) call kspsetup (ksp_cg3d, ierr) CASE DEFAULT PRINT *, "PETSC PRECOND: SELECTED DEFAULT" call PCSetType(pc_cg3d, PCBJACOBI, ierr) END SELECT CALL TIMER_STOP ('CG3D_PETSC_MATCREATE',myThid) endif ! need3create_mat CALL TIMER_START ('CG3D_PETSC_SOLVE',myThid) call KSPSolve(ksp_cg3d, rhs, solution, ierr) CALL TIMER_STOP ('CG3D_PETSC_SOLVE',myThid) call KSPGetIterationNumber(ksp_cg3d,iters,ierr) maxIter = iters call VecGetValues(solution,local_dofs,indices, & solution_values,ierr) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx DO k=1,Nr dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj)) & - local_offset color = INT(cg3d_petsc_color(i,j,k,bi,bj)) rank = cg3d_color_rank (color) if (dof_index.ge.0.and.rank.eq.mpimywid) THEN cg3d_x(i,j,k,bi,bj) = solution_values(dof_index+1) endif ENDDO ENDDO ENDDO ENDDO ENDDO if (cg3d_need2destroy_mat) then call KSPDestroy (ksp_cg3d, ierr) call MatDestroy (matrix_petsc_cg3d, ierr) endif call VecDestroy (rhs, ierr) call VecDestroy (solution, ierr) #endif RETURN END