diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 83980b649bfbfa066bb1974b4635c90e46627abc..f5d38df5761efb79d8c2b6dd133d9a3b26c15697 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -93,7 +93,7 @@ subroutine geomcomplex_ellip(epsi,epsi_x,epsi_y,epsi_z,nxi,nxf,ny,nyi,nyf,nzi,nz
epsi_x(i,j,k)=pointVelocity(1)
epsi_y(i,j,k)=pointVelocity(2)
epsi_z(i,j,k)=pointVelocity(3)
- ! write(*,*) remp
+ ! write(*,*) 'Set ep1_x(',i,j,k,') = ',pointVelocity(1)
enddo
enddo
enddo
diff --git a/src/genepsi3d.f90 b/src/genepsi3d.f90
index 78e4b2e2ef2f835a341012bceb22a588ea6c3935..0876c090fff86888081359c1ae406ed19795e350 100644
--- a/src/genepsi3d.f90
+++ b/src/genepsi3d.f90
@@ -175,9 +175,9 @@ contains
real(mytype),dimension(nobjmax,xsize(2),xsize(3)) :: xi,xf
real(mytype),dimension(nobjmax,ysize(1),ysize(3)) :: yi,yf
real(mytype),dimension(nobjmax,zsize(1),zsize(2)) :: zi,zf
- real(mytype),dimension(nxraf,xsize(2),xsize(3)) :: xepsi
- real(mytype),dimension(ysize(1),nyraf,ysize(3)) :: yepsi
- real(mytype),dimension(zsize(1),zsize(2),nzraf) :: zepsi
+ real(mytype),dimension(nxraf,xsize(2),xsize(3)) :: xepsi,xep1,xep2,xep3
+ real(mytype),dimension(ysize(1),nyraf,ysize(3)) :: yepsi,yep1,yep2,yep3
+ real(mytype),dimension(zsize(1),zsize(2),nzraf) :: zepsi,zep1,zep2,zep3
real(mytype),dimension(ny) :: yp
real(mytype),dimension(nyraf) :: ypraf
real(mytype) :: dxraf,dyraf,dzraf
@@ -198,7 +198,7 @@ contains
integer :: numvis
integer :: mpi_aux_i, code
- ! write(*,*)'Inside gene_epsi_3D'
+ write(*,*)'Inside gene_epsi_3D'
!x-pencil
ep1=zero
call geomcomplex(ep1,ep1_ux,ep1_uy,ep1_uz,xstart(1),xend(1),ny,xstart(2),xend(2),xstart(3),xend(3),dx,yp,dz,one)
@@ -209,7 +209,8 @@ contains
dxraf =xlx/real(nxraf-1, mytype)
endif
xepsi=zero
- call geomcomplex(xepsi,ep1_ux,ep1_uy,ep1_uz,1,nxraf,ny,xstart(2),xend(2),xstart(3),xend(3),dxraf,yp,dz,one)
+ write(*,*) 'Doing x pencil'
+ call geomcomplex(xepsi,xep1,xep2,xep3,1,nxraf,ny,xstart(2),xend(2),xstart(3),xend(3),dxraf,yp,dz,one)
! if (nrank==0) print*,' step 2'
!y-pencil
if(ncly)then
@@ -224,7 +225,7 @@ contains
enddo
if(.not.ncly)ypraf(nyraf)=yp(ny)
yepsi=zero
- call geomcomplex(yepsi,ep1_ux,ep1_uy,ep1_uz,ystart(1),yend(1),nyraf,1,nyraf,ystart(3),yend(3),dx,ypraf,dz,one)
+ call geomcomplex(yepsi,yep1,yep2,yep3,ystart(1),yend(1),nyraf,1,nyraf,ystart(3),yend(3),dx,ypraf,dz,one)
! if (nrank==0) print*,' step 3'
!z-pencil
if(nclz)then
@@ -233,7 +234,7 @@ contains
dzraf=zlz/real(nzraf-1, mytype)
endif
zepsi=zero
- call geomcomplex(zepsi,ep1_ux,ep1_uy,ep1_uz,zstart(1),zend(1),ny,zstart(2),zend(2),1,nzraf,dx,yp,dzraf,one)
+ call geomcomplex(zepsi,zep1,zep2,zep3,zstart(1),zend(1),ny,zstart(2),zend(2),1,nzraf,dx,yp,dzraf,one)
! if (nrank==0) print*,' step 4'
!x-pencil
diff --git a/src/navier.f90 b/src/navier.f90
index ce849661ca20873e0c90ca149c8d044487b881e8..0c405e50d84979b67354ac497470af1aeec208e5 100644
--- a/src/navier.f90
+++ b/src/navier.f90
@@ -20,7 +20,7 @@ contains
!! DESCRIPTION: Takes the intermediate momentum field as input,
!! computes div and solves pressure-Poisson equation.
!############################################################################
- SUBROUTINE solve_poisson(pp3, px1, py1, pz1, rho1, ux1, uy1, uz1, ep1, drho1, divu3)
+ SUBROUTINE solve_poisson(pp3, px1, py1, pz1, rho1, ux1, uy1, uz1, ep1, ep1x,ep1y,ep1z, drho1, divu3)
USE decomp_2d, ONLY : mytype, xsize, zsize, ph1, nrank, real_type
USE decomp_2d_poisson, ONLY : poisson
@@ -34,7 +34,7 @@ contains
!! Inputs
REAL(mytype), DIMENSION(xsize(1), xsize(2), xsize(3)), INTENT(IN) :: ux1, uy1, uz1
- REAL(mytype), DIMENSION(xsize(1), xsize(2), xsize(3)), INTENT(IN) :: ep1
+ REAL(mytype), DIMENSION(xsize(1), xsize(2), xsize(3)), INTENT(IN) :: ep1,ep1x,ep1y,ep1z
REAL(mytype), DIMENSION(xsize(1), xsize(2), xsize(3), nrhotime), INTENT(IN) :: rho1
REAL(mytype), DIMENSION(xsize(1), xsize(2), xsize(3), ntime), INTENT(IN) :: drho1
REAL(mytype), DIMENSION(zsize(1), zsize(2), zsize(3)), INTENT(IN) :: divu3
@@ -70,7 +70,7 @@ contains
CALL momentum_to_velocity(rho1, ux1, uy1, uz1)
ENDIF
- CALL divergence(pp3(:,:,:,1),rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock)
+ CALL divergence(pp3(:,:,:,1),rho1,ux1,uy1,uz1,ep1,ep1x,ep1y,ep1z,drho1,divu3,nlock)
IF (ilmn.AND.ivarcoeff) THEN
dv3(:,:,:) = pp3(:,:,:,1)
ENDIF
@@ -82,7 +82,7 @@ contains
IF (.NOT.converged) THEN
!! Evaluate additional RHS terms
- CALL calc_varcoeff_rhs(pp3(:,:,:,1), rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rho0, &
+ CALL calc_varcoeff_rhs(pp3(:,:,:,1), rho1, px1, py1, pz1, dv3, drho1, ep1,ep1x,ep1y,ep1z, divu3, rho0, &
poissiter)
ENDIF
@@ -254,7 +254,7 @@ contains
! output : pp3 (on pressure mesh)
!written by SL 2018
!############################################################################
- subroutine divergence (pp3,rho1,ux1,uy1,uz1,ep1,drho1,divu3,nlock)
+ subroutine divergence (pp3,rho1,ux1,uy1,uz1,ep1,ep1x,ep1y,ep1z,drho1,divu3,nlock)
USE param
USE decomp_2d
@@ -270,7 +270,7 @@ contains
! TYPE(DECOMP_INFO) :: ph1,ph3,ph4
!X PENCILS NX NY NZ -->NXM NY NZ
- real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1,uy1,uz1,ep1
+ real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(in) :: ux1,uy1,uz1,ep1,ep1x,ep1y,ep1z
real(mytype),dimension(xsize(1),xsize(2),xsize(3),ntime),intent(in) :: drho1
real(mytype),dimension(xsize(1),xsize(2),xsize(3),nrhotime),intent(in) :: rho1
!Z PENCILS NXM NYM NZ -->NXM NYM NZM
@@ -289,9 +289,9 @@ contains
tc1(:,:,:) = uz1(:,:,:)
elseif (itype.eq.itype_ellip) then
write(*,*) 'Computing flow field in navier correctly.'
- ta1(:,:,:) = (one - ep1(:,:,:)) * ux1(:,:,:) + ep1(:,:,:)*ep1_ux(:,:,:)
- tb1(:,:,:) = (one - ep1(:,:,:)) * uy1(:,:,:) + ep1(:,:,:)*ep1_uy(:,:,:)
- tc1(:,:,:) = (one - ep1(:,:,:)) * uz1(:,:,:) + ep1(:,:,:)*ep1_uz(:,:,:)
+ ta1(:,:,:) = (one - ep1(:,:,:)) * ux1(:,:,:) + ep1(:,:,:)*ep1x(:,:,:)
+ tb1(:,:,:) = (one - ep1(:,:,:)) * uy1(:,:,:) + ep1(:,:,:)*ep1y(:,:,:)
+ tc1(:,:,:) = (one - ep1(:,:,:)) * uz1(:,:,:) + ep1(:,:,:)*ep1z(:,:,:)
else
ta1(:,:,:) = (one - ep1(:,:,:)) * ux1(:,:,:) + ep1(:,:,:)*lvx
tb1(:,:,:) = (one - ep1(:,:,:)) * uy1(:,:,:) + ep1(:,:,:)*lvy
@@ -1152,7 +1152,7 @@ contains
!! DESCRIPTION: Computes RHS of the variable-coefficient Poisson solver
!!
!############################################################################
- SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, divu3, rho0, poissiter)
+ SUBROUTINE calc_varcoeff_rhs(pp3, rho1, px1, py1, pz1, dv3, drho1, ep1, ep1x,ep1y,ep1z, divu3, rho0, poissiter)
USE MPI
@@ -1171,7 +1171,7 @@ contains
REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3)) :: px1, py1, pz1
REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3), nrhotime) :: rho1
REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3), ntime) :: drho1
- REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3)) :: ep1
+ REAL(mytype), INTENT(IN), DIMENSION(xsize(1), xsize(2), xsize(3)) :: ep1,ep1x,ep1y,ep1z
REAL(mytype), INTENT(IN), DIMENSION(zsize(1), zsize(2), zsize(3)) :: divu3
REAL(mytype), INTENT(IN), DIMENSION(ph1%zst(1):ph1%zen(1), ph1%zst(2):ph1%zen(2), nzmsize) :: dv3
real(mytype) :: rho0
@@ -1195,7 +1195,7 @@ contains
tc1(:,:,:) = (one - rho0 / rho1(:,:,:,1)) * pz1(:,:,:)
nlock = -1 !! Don't do any funny business with LMN
- CALL divergence(pp3,rho1,ta1,tb1,tc1,ep1,drho1,divu3,nlock)
+ CALL divergence(pp3,rho1,ta1,tb1,tc1,ep1,ep1x,ep1y,ep1z,drho1,divu3,nlock)
!! lapl(p) = div((1 - rho0/rho) grad(p)) + rho0(div(u*) - div(u))
!! dv3 contains div(u*) - div(u)
diff --git a/src/variables.f90 b/src/variables.f90
index bb3917716d1b767c97dceb664ee275588ecf2df8..437cfb955d0d23f0ce14af02654d8f8095f59044 100644
--- a/src/variables.f90
+++ b/src/variables.f90
@@ -22,6 +22,7 @@ module var
real(mytype), save, allocatable, dimension(:,:,:,:) :: phi1, phi2, phi3
real(mytype), save, allocatable, dimension(:,:,:) :: px1, py1, pz1
real(mytype), save, allocatable, dimension(:,:,:) :: ep1, diss1, pre1
+ real(mytype), save, allocatable, dimension(:,:,:) :: ep1_ux, ep1_uy, ep1_uz
real(mytype), save, allocatable, dimension(:,:,:,:) :: dux1,duy1,duz1 ! Output of convdiff
real(mytype), save, allocatable, dimension(:,:,:,:,:) :: dphi1
real(mytype), save, allocatable, dimension(:,:,:) :: mu1,mu2,mu3
@@ -156,6 +157,12 @@ contains
di1 = zero
call alloc_x(ep1)
ep1 = zero
+ call alloc_x(ep1_ux)
+ ep1_ux = zero
+ call alloc_x(ep1_uy)
+ ep1_uy = zero
+ call alloc_x(ep1_uz)
+ ep1_uz = zero
if (ilmn) then
call alloc_x(mu1, opt_global=.true.)
mu1 = one
diff --git a/src/xcompact3d.f90 b/src/xcompact3d.f90
index ded7abbdf634aa7ff29e554833b15fc45e0484f6..139a51342bd7ea091c12bb8a1660e2267d7a1f02 100644
--- a/src/xcompact3d.f90
+++ b/src/xcompact3d.f90
@@ -64,7 +64,7 @@ program xcompact3d
call pre_correc(ux1,uy1,uz1,ep1)
call calc_divu_constraint(divu3,rho1,phi1)
- call solve_poisson(pp3,px1,py1,pz1,rho1,ux1,uy1,uz1,ep1,drho1,divu3)
+ call solve_poisson(pp3,px1,py1,pz1,rho1,ux1,uy1,uz1,ep1,ep1_ux,ep1_uy,ep1_uz,drho1,divu3)
call cor_vel(ux1,uy1,uz1,px1,py1,pz1)
if (ilmn) then