Skip to content
Snippets Groups Projects
Commit 11437b1e authored by Juan Escobar's avatar Juan Escobar
Browse files

Juan 18/12/2018: version with CONRES & RICH workinf with IT=11, 128x128 3D+REU... 

Juan 18/12/2018: version with CONRES & RICH workinf with IT=11, 128x128 3D+REU , 4P , 4* solver than ZRESI
parent 52a2a20b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 145 additions and 33 deletions
zsolver_inv.f90
+ 145
− 33
View file @ 11437b1e
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme