diff --git a/zsolver_inv.f90 b/zsolver_inv.f90
index 3af2003ec0164c1185df5521cfbeceb0efa853e6..27e8b34be802160337cb183cad6af1fab621d5a6 100644
--- a/zsolver_inv.f90
+++ b/zsolver_inv.f90
@@ -144,6 +144,7 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
USE MODI_FFT55
USE MODI_GET_HALO
USE MODI_FLAT_INV
+ USE MODI_DOTPROD
!
IMPLICIT NONE
!
@@ -241,11 +242,23 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZBAND_YT ! array in Y slices distribution transpose
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZBAND_YRT ! array in Y slices distribution transpose
!
- REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZF_1_Y,ZCORREC
+ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZF_1_Y,ZCORREC,ZP
+ !
+ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRESIDUE,ZDELTA,ZQ,ZKSI
+ !
+ REAL :: ZDOT_DELTA,ZLAMBDA,ZALPHA
+ !
+ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZY_RES_FLAT,ZF_1_Y_FLAT
+ !
+ INTEGER :: JM , KITR
+ !
+ LOGICAL :: GRICHARDSON
+ !
+ LOGICAL,SAVE :: GFIRST_CALL_ZSOL = .TRUE.
!
INTEGER :: IIU,IJU
REAL :: ZMEAN
- INTEGER :: IT , NT
+ INTEGER :: IT , NT
REAL :: ZOMEGA
!-------------------------------------------------------------------------------
!
@@ -279,13 +292,21 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
END IF
ALLOCATE(ZF_1_Y(IIU,IJU,IKU))
ALLOCATE(ZCORREC(IIU,IJU,IKU))
+ ALLOCATE(ZP(IIU,IJU,IKU))
+ ALLOCATE(ZRESIDUE(IIU,IJU,IKU))
+ ALLOCATE(ZDELTA(IIU,IJU,IKU))
+ ALLOCATE(ZQ(IIU,IJU,IKU))
+ ALLOCATE(ZKSI(IIU,IJU,IKU))
+ ALLOCATE(ZY_RES_FLAT(IIU,IJU,IKU))
+ ALLOCATE(ZF_1_Y_FLAT(IIU,IJU,IKU))
!
+ ZDXM2 = PDXHATM*PDXHATM
+ ZDYM2 = PDYHATM*PDYHATM
!
-!!$ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &! -1
-!!$ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,ZF_1_Y) ! F (Y - Q (PHI)))
-!!$
-!!$ CALL PF(ZY,ZF_1_Y)
-!!$ ZY = ZY - PY
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,ZF_1_Y_FLAT) ! F (Y - Q (PHI)))
+ CALL PF(ZY_RES_FLAT,ZF_1_Y_FLAT)
+ ZY_RES_FLAT = ZY_RES_FLAT - PY
!
!-------------------------------------------------------------------------------
!
@@ -295,7 +316,7 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
!
!* 3.1 copy the RHS in a local array REMAP functions will shift the indices for the FFT
!
- PF_1_Y = 0.
+ !PF_1_Y = 0.
ZY = PY
!
!* 3.2 form homogeneous boundary condition used by the FFT for non-periodic
@@ -347,10 +368,6 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
!* 5. MATRIX INVERSION FOR THE FLAT OPERATOR
! --------------------------------------
!
-
- ZDXM2 = PDXHATM*PDXHATM
- ZDYM2 = PDYHATM*PDYHATM
-
IF (L2D) THEN
CALL FAST_SUBSTITUTION_2D(ZBAND_YR,ZBETX,PBFB,ZGAM,PCF,ZAF &
,ZBAND_Y,IIY,IJY,IKU)
@@ -358,38 +375,128 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
CALL FAST_SUBSTITUTION_3D(ZF_1_Y,ZBETX,PBFB,ZGAM,PCF,PAF &
,ZY,IIU,IJU,IKU)
-
CALL PF_1_Y_BOUND(ZF_1_Y)
+ CALL GET_HALO(ZF_1_Y)
- NT = 0 ; ZOMEGA = 1.0 ! FLAT_INV 4IT
- NT = 11 ; ZOMEGA = 1.0 ! 12IT
- !NT = 101 ; ZOMEGA = 1.0 ! 6IT
- NT = 301 ; ZOMEGA = 1.0 ! 5IT
-
- PF_1_Y = ZF_1_Y ! when no first guess is available, we take the solution
- ! for the flat problem
+ IF (GFIRST_CALL_ZSOL) THEN
+ !GFIRST_CALL_ZSOL = .FALSE.
+ PF_1_Y = ZF_1_Y ! when no first guess is available, we take the solution
+ ! for the flat problem
+ END IF
!
- CALL PF_1_Y_BOUND(PF_1_Y)
!
+GRICHARDSON = .TRUE. ! .FALSE. ! .TRUE.
+ ! FLAT_INV :: 4IT 3D_128X128 / 5IT REU 128x128
+ KITR = 11 ! CONRES :: 16IT 3D 128X128 / 7IT REU 128x128
+ NT = 11 ; ZOMEGA = 1.0 ! RICH :: 36IT 3D 128X128 / 26IT REU 128x128
+IF ( GRICHARDSON ) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!! RICHARDSON
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !NT = 1000 ; ZOMEGA = 1.0
+ !NT = 11 ; ZOMEGA = 1.0 ! 12IT
+ !NT = 101 ; ZOMEGA = 1.0 ! 6IT
+ !NT = 401 ; ZOMEGA = 1.0 ! 1.0 ! 5IT
+
DO IT = 1,NT
!
- CALL PF(ZY,PF_1_Y) ! Q (PF_1_Y)
+ CALL PF(ZP,PF_1_Y) ! Q (PF_1_Y)
!
ZCORREC = 0.
!
CALL FAST_SUBSTITUTION_3D(ZCORREC,ZBETX,PBFB,ZGAM,PCF,PAF &
- ,ZY,IIU,IJU,IKU) ! zcorrec = F-1 * Q (PF_1_Y)
+ ,ZP,IIU,IJU,IKU) ! zcorrec = F-1 * Q (PF_1_Y)
CALL PF_1_Y_BOUND(ZCORREC)
! -1 -1
! update the iterative solution PHI = PHI + relax* (F (Y) - F * Q (PHI))
!
- ZY = ZF_1_Y - ZCORREC ! for totalview
+ ZP = ZF_1_Y - ZCORREC ! for totalview
PF_1_Y = PF_1_Y + ZOMEGA * (ZF_1_Y - ZCORREC)
+ CALL PF_1_Y_BOUND(PF_1_Y)
CALL GET_HALO(PF_1_Y)
- !CALL PF_1_Y_BOUND(PF_1_Y)
!
END DO
+ELSE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!! CONJUGE RESIDUAL
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!
+!* 1.1 compute the vector: r^(0) = Q(PHI) - Y
+!
+CALL PF(ZRESIDUE,PF_1_Y)
+ZRESIDUE = ZRESIDUE - PY
+!
+!* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y )
+!
+CALL FAST_SUBSTITUTION_3D(ZP,ZBETX,PBFB,ZGAM,PCF,PAF &
+ ,ZRESIDUE,IIU,IJU,IKU)
+CALL PF_1_Y_BOUND(ZP)
+CALL GET_HALO(ZP)
+!
+!* 1.3 compute the vector: delta^(0) = Q ( p^(0) )
+!
+CALL PF(ZDELTA,ZP)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. ITERATIVE LOOP
+! --------------
+!
+DO JM = 1,KITR
+ !
+ !* 2.1 compute the step LAMBDA
+ !
+ ZDOT_DELTA = DOTPROD(ZDELTA, ZDELTA,HLBCX,HLBCY) ! norm of DELTA
+ ZLAMBDA = - DOTPROD(ZRESIDUE,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA
+ !
+ !* 2.2 update the pressure function PHI
+ !
+ PF_1_Y = PF_1_Y + ZLAMBDA * ZP
+ CALL PF_1_Y_BOUND(PF_1_Y)
+ CALL GET_HALO(PF_1_Y)
+ !
+ !
+ IF( JM == KITR ) EXIT
+ !
+ !
+ !* 2.3 update the residual error: r
+ !
+ ZRESIDUE = ZRESIDUE + ZLAMBDA * ZDELTA
+ !
+ !* 2.4 compute the vector: q = F^(-1)*( Q(PHI) - Y )
+ !
+ !CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ ! PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZQ, &
+ ! PBFB,&
+ ! PBF_SXP2_YP1_Z) !JUAN Z_SPLITTING
+ CALL FAST_SUBSTITUTION_3D(ZQ,ZBETX,PBFB,ZGAM,PCF,PAF &
+ ,ZRESIDUE,IIU,IJU,IKU)
+ CALL PF_1_Y_BOUND(ZQ)
+ CALL GET_HALO(ZQ)
+ !
+ !* 2.5 compute the auxiliary field: ksi = Q ( q )
+ !
+ ! ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ)
+ CALL PF(ZKSI,ZQ)
+ ! -1
+ !* 2.6 compute the step ALPHA
+ !
+ ZALPHA = - DOTPROD(ZKSI,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA ! lambda
+ !
+ !* 2.7 update p and DELTA
+ !
+ ZP = ZQ + ZALPHA * ZP
+ ZDELTA = ZKSI + ZALPHA * ZDELTA
+ !
+END DO ! end of the loop for the iterative solver
+ENDIF
!!$ZMEAN = SUM ( PF_1_Y(IIB:IIE,IJB:IJE,IKU) ) / FLOAT( (IIE-IIB+1)*(IJE-IJB+1) )
!!$PF_1_Y(IIB:IIE,IJB:IJE,:) = PF_1_Y(IIB:IIE,IJB:IJE,:) - ZMEAN
!CALL PF_1_Y_BOUND(PF_1_Y)
@@ -419,6 +526,9 @@ END DO
!
!PF_1_Y = ZF_1_Y
!
+ CALL PF(ZRESIDUE,PF_1_Y) ! for totalview
+ ZRESIDUE = ZRESIDUE - ZY
+ !
CALL GET_HALO(PF_1_Y)
!-------------------------------------------------------------------------------
!
@@ -442,11 +552,13 @@ DO JK=IKB,IKE
END DO
END DO
END DO
-ZY(IIB:IIE,IJB:IJE,IKB-1) = PBFB(IIB:IIE,IJB:IJE,IKB-1) * PF_1_Y(IIB:IIE,IJB:IJE,IKB-1) &
+PY(IIB:IIE,IJB:IJE,IKB-1) = PBFB(IIB:IIE,IJB:IJE,IKB-1) * PF_1_Y(IIB:IIE,IJB:IJE,IKB-1) &
+ PCF(IKB-1) * PF_1_Y(IIB:IIE,IJB:IJE,IKB)
-ZY(IIB:IIE,IJB:IJE,IKE+1) = PAF(IKE+1) * PF_1_Y(IIB:IIE,IJB:IJE,IKE) &
- + PBFB(IIB:IIE,IJB:IJE,IKE+1) * PF_1_Y(IIB:IIE,IJB:IJE,IKE+1)
-
+PY(IIB:IIE,IJB:IJE,IKE+1) = PAF(IKE+1) * PF_1_Y(IIB:IIE,IJB:IJE,IKE) &
+ + PBFB(IIB:IIE,IJB:IJE,IKE+1) * PF_1_Y(IIB:IIE,IJB:IJE,IKE+1)
+!
+CALL GET_HALO(PY)
+!
END SUBROUTINE PF
SUBROUTINE PF_1_Y_BOUND(PF_1_Y)
@@ -561,19 +673,19 @@ REAL , DIMENSION (:,:,:) :: PF_1_Y
IF (.NOT. L2D .AND. HLBCX(1)/='CYCL' .AND. HLBCY(1)/='CYCL') THEN
! the following verticals are not used
IF ( (LWEST_ll(HSPLITTING='B')).AND.(LSOUTH_ll(HSPLITTING='B')) ) THEN
- PF_1_Y(IIB-1,IJB-1,:)=PF_1_Y(IIB,IJB,:) ! 0.0
+ PF_1_Y(IIB-1,IJB-1,:)= PF_1_Y(IIB,IJB,:) ! 0.0
END IF
!
IF ( (LWEST_ll(HSPLITTING='B')).AND.(LNORTH_ll(HSPLITTING='B')) ) THEN
- PF_1_Y(IIB-1,IJE+1,:)=PF_1_Y(IIB,IJE,:) ! 0.0
+ PF_1_Y(IIB-1,IJE+1,:)= PF_1_Y(IIB,IJE,:) ! 0.0
END IF
!
IF ( (LEAST_ll(HSPLITTING='B')).AND.(LSOUTH_ll(HSPLITTING='B')) ) THEN
- PF_1_Y(IIE+1,IJB-1,:)=PF_1_Y(IIE,IJB,:) ! 0.0
+ PF_1_Y(IIE+1,IJB-1,:)= PF_1_Y(IIE,IJB,:) ! 0.0
END IF
!
IF ( (LEAST_ll(HSPLITTING='B')).AND.(LNORTH_ll(HSPLITTING='B')) ) THEN
- PF_1_Y(IIE+1,IJE+1,:)=PF_1_Y(IIE,IJE,:) ! 0.0
+ PF_1_Y(IIE+1,IJE+1,:)= PF_1_Y(IIE,IJE,:) ! 0.0
END IF
END IF
END SUBROUTINE PF_1_Y_BOUND