C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_cg_solve.F,v 1.9 2015/02/21 19:09:53 dgoldberg Exp $ C $Name: $ #include "STREAMICE_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE STREAMICE_CG_SOLVE( U cg_Uin, ! x-velocities U cg_Vin, ! y-velocities I cg_Bu, ! force in x dir I cg_Bv, ! force in y dir I A_uu, ! section of matrix that multiplies u and projects on u I A_uv, ! section of matrix that multiplies v and projects on u I A_vu, ! section of matrix that multiplies u and projects on v I A_vv, ! section of matrix that multiplies v and projects on v I tolerance, O iters, I maxIter, I myThid ) C /============================================================\ C | SUBROUTINE | C | o | C |============================================================| C | | C \============================================================/ IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "STREAMICE.h" #include "STREAMICE_CG.h" !#ifdef ALLOW_PETSC !#include "finclude/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 iters INTEGER maxIter _RL tolerance _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1) C LOCAL VARIABLES INTEGER i, j, bi, bj, cg_halo, conv_flag INTEGER iter, is, js, ie, je, colx, coly, 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 !#ifdef ALLOW_PETSC ! INTEGER indices(2*(snx*nsx*sny*nsy)) ! INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1) ! _RL rhs_values(2*(snx*nsx*sny*nsy)) ! _RL solution_values(2*(snx*nsx*sny*nsy)) ! _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy)) ! _RL mat_values (18,1), mat_val_return(1) ! INTEGER indices_col(18) ! INTEGER local_dofs, global_dofs, dof_index, dof_index_col ! INTEGER local_offset ! Mat matrix ! KSP ksp ! PC pc ! Vec rhs ! Vec solution ! PetscErrorCode ierr !#ifdef ALLOW_USE_MPI ! integer mpiRC, mpiMyWid !#endif !#endif #ifdef ALLOW_STREAMICE CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid) #ifndef STREAMICE_SERIAL_TRISOLVE #ifdef ALLOW_PETSC if (streamice_use_petsc) then CALL STREAMICE_CG_SOLVE_PETSC( U cg_Uin, ! x-velocities U cg_Vin, ! y-velocities I cg_Bu, ! force in x dir I cg_Bv, ! force in y dir I A_uu, ! section of matrix that multiplies u and projects on u I A_uv, ! section of matrix that multiplies v and projects on u I A_vu, ! section of matrix that multiplies u and projects on v I A_vv, ! section of matrix that multiplies v and projects on v I tolerance, I iters, O maxiter, I myThid ) else #endif /* ALLOW_PETSC */ iters = maxIter conv_flag = 0 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx Zu_SI (i,j,bi,bj) = 0. _d 0 Zv_SI (i,j,bi,bj) = 0. _d 0 Ru_SI (i,j,bi,bj) = 0. _d 0 Rv_SI (i,j,bi,bj) = 0. _d 0 Au_SI (i,j,bi,bj) = 0. _d 0 Av_SI (i,j,bi,bj) = 0. _d 0 Du_SI (i,j,bi,bj) = 0. _d 0 Dv_SI (i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO C FIND INITIAL RESIDUAL, and initialize r ! #ifdef STREAMICE_CONSTRUCT_MATRIX DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx DO colx=-1,1 DO coly=-1,1 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + & A_uu(i,j,bi,bj,colx,coly)* & cg_Uin(i+colx,j+coly,bi,bj)+ & A_uv(i,j,bi,bj,colx,coly)* & cg_Vin(i+colx,j+coly,bi,bj) Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + & A_vu(i,j,bi,bj,colx,coly)* & cg_Uin(i+colx,j+coly,bi,bj)+ & A_vv(i,j,bi,bj,colx,coly)* & cg_Vin(i+colx,j+coly,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO _EXCH_XY_RL( Au_SI, myThid ) _EXCH_XY_RL( Av_SI, myThid ) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)- & Au_SI(i,j,bi,bj) Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)- & Av_SI(i,j,bi,bj) ENDDO ENDDO dot_p1_tile(bi,bj) = 0. _d 0 dot_p2_tile(bi,bj) = 0. _d 0 ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) resid_0 = sqrt(dot_p1) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) WRITE(msgBuf,'(A,I1,I1,E14.7)') 'CONJ GRAD INIT RESID LOCAL, ', & bi,bj, dot_p1_tile(bi,bj) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) enddo enddo WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ', & resid_0 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) C CCCCCCCCCCCCCCCCCCCC DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj) IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO cg_halo = min(OLx-1,OLy-1) conv_flag = 0 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj) Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO resid = resid_0 iters = 0 c !!!!!!!!!!!!!!!!!! c !! !! c !! MAIN CG LOOP !! c !! !! c !!!!!!!!!!!!!!!!!! c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!) WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid) do iter = 1, maxIter if (resid .gt. tolerance*resid_0) then c to avoid using "exit" iters = iters + 1 is = 1 - cg_halo ie = sNx + cg_halo js = 1 - cg_halo je = sNy + cg_halo DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Au_SI(i,j,bi,bj) = 0. _d 0 Av_SI(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO ! IF (STREAMICE_construct_matrix) THEN ! #ifdef STREAMICE_CONSTRUCT_MATRIX DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=js,je DO i=is,ie DO colx=-1,1 DO coly=-1,1 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + & A_uu(i,j,bi,bj,colx,coly)* & Du_SI(i+colx,j+coly,bi,bj)+ & A_uv(i,j,bi,bj,colx,coly)* & Dv_SI(i+colx,j+coly,bi,bj) Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + & A_vu(i,j,bi,bj,colx,coly)* & Du_SI(i+colx,j+coly,bi,bj)+ & A_vv(i,j,bi,bj,colx,coly)* & Dv_SI(i+colx,j+coly,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO ! else ! #else ! ! CALL STREAMICE_CG_ACTION( myThid, ! O Au_SI, ! O Av_SI, ! I Du_SI, ! I Dv_SI, ! I is,ie,js,je) ! ! ! ENDIF ! ! #endif DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) dot_p1_tile(bi,bj) = 0. _d 0 dot_p2_tile(bi,bj) = 0. _d 0 ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* & Ru_SI(i,j,bi,bj) dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)* & Au_SI(i,j,bi,bj) ENDIF IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* & Rv_SI(i,j,bi,bj) dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)* & Av_SI(i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) alpha_k = dot_p1/dot_p2 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+ & alpha_k*Du_SI(i,j,bi,bj) Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj) Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)- & alpha_k*Au_SI(i,j,bi,bj) Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) / & DIAGu_SI(i,j,bi,bj) ENDIF IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+ & alpha_k*Dv_SI(i,j,bi,bj) Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj) Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)- & alpha_k*Av_SI(i,j,bi,bj) Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) / & DIAGv_SI(i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) dot_p1_tile(bi,bj) = 0. _d 0 dot_p2_tile(bi,bj) = 0. _d 0 ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* & Ru_SI(i,j,bi,bj) dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)* & Ru_old_SI(i,j,bi,bj) ENDIF IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* & Rv_SI(i,j,bi,bj) dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)* & Rv_old_SI(i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) beta_k = dot_p1/dot_p2 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+ & Zu_SI(i,j,bi,bj) IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+ & Zv_SI(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) dot_p1_tile(bi,bj) = 0. _d 0 ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) resid = sqrt(dot_p1) ! IF (iter .eq. 1) then ! print *, alpha_k, beta_k, resid ! ENDIF cg_halo = cg_halo - 1 if (cg_halo .eq. 0) then cg_halo = min(OLx-1,OLy-1) _EXCH_XY_RL( Du_SI, myThid ) _EXCH_XY_RL( Dv_SI, myThid ) _EXCH_XY_RL( Ru_SI, myThid ) _EXCH_XY_RL( Rv_SI, myThid ) _EXCH_XY_RL( cg_Uin, myThid ) _EXCH_XY_RL( cg_Vin, myThid ) endif endif enddo ! end of CG loop c to avoid using "exit" c if iters has reached max_iters there is no convergence IF (iters .lt. maxIter) THEN conv_flag = 1 ENDIF PRINT *, "GOT HERE CG ITERATIONS", iters ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1-OLy,sNy+OLy ! DO i=1-OLy,sNx+OLy ! IF (STREAMICE_umask(i,j,bi,bj).eq.3.0) ! & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj) ! IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0) ! & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! _EXCH_XY_RL( cg_Uin, myThid ) ! _EXCH_XY_RL( cg_Vin, myThid ) #ifdef ALLOW_PETSC endif !if (streamice_use_petsc) #endif #else /* STREAMICE_SERIAL_TRISOLVE */ iters = 0 CALL STREAMICE_TRIDIAG_SOLVE( U cg_Uin, ! x-velocities U cg_Vin, U cg_Bu, ! force in x dir I A_uu, ! section of matrix that multiplies u and projects on u I STREAMICE_umask, I myThid ) #endif CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid) #endif RETURN END