From fe6444366cab9954a1152c5136d6803052f6849a Mon Sep 17 00:00:00 2001 From: park0d Date: Tue, 21 Mar 2017 01:03:39 +0900 Subject: [PATCH] add scalar equation, disable reaction --- Makefile | 2 +- m_fdm_calc.f90 | 154 ++++++++++++++++++++++++++++++------------------- 2 files changed, 95 insertions(+), 61 deletions(-) diff --git a/Makefile b/Makefile index d1635eb..439ca93 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,7 @@ MPIF90 = blah - FFTW_HOME = ../libs/fftw-3.2.2 + FFTW_HOME = ../libs/fftw MPIF90 = mpif90 FCFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include LDFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 375ca44..ae8dd50 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -8,6 +8,11 @@ module m_fdm_calc !variables + integer, parameter :: nscalar = 2 + integer, parameter :: neq = 1 + nscalar + integer, parameter :: nd1 = 2 + 2 * nscalar + integer, parameter :: nd2 = nscalar + real*8, dimension(:,:,:), allocatable :: u_,v_,w_ real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf real*8, dimension(:,:), allocatable :: fz, dfz, fzz, dfzz @@ -234,20 +239,20 @@ module m_fdm_calc fdmtime=0. fullsavenum=1000 !full save file - allocate(y1(2,nx,ny,nz)) - allocate(y2(2,nx,ny,nz)) - allocate(yf(2,nx,ny,nz)) + allocate(y1(neq,nx,ny,nz)) + allocate(y2(neq,nx,ny,nz)) + allocate(yf(neq,nx,ny,nz)) - allocate(fz(4*nx*ny,nz)) - allocate(dfz(4*nx*ny,nz)) - allocate(fzz(nx*ny,nz)) - allocate(dfzz(nx*ny,nz)) + allocate( fz(nd1*nx*ny,nz)) + allocate(dfz(nd1*nx*ny,nz)) + allocate( fzz(nd2*nx*ny,nz)) + allocate(dfzz(nd2*nx*ny,nz)) - allocate(fzu(4*nx*ny,2)) - allocate(fzl(4*nx*ny,2)) + allocate(fzu(nd1*nx*ny,2)) + allocate(fzl(nd1*nx*ny,2)) - allocate(fzzu(nx*ny,2)) - allocate(fzzl(nx*ny,2)) + allocate(fzzu(nd2*nx*ny,2)) + allocate(fzzl(nd2*nx*ny,2)) y1=0.0 y2=0.0 @@ -280,6 +285,7 @@ module m_fdm_calc do ii=1,nx y1(1,ii,i,j)=1. ! rho initializing y1(2,ii,i,j)=yy(ii) ! Yr initializing + y1(3,ii,i,j)=0.0 ! Yr initializing enddo enddo enddo @@ -467,16 +473,16 @@ module m_fdm_calc integer :: i,j,k,xx,yy,zz,ii integer :: n - integer :: idx1, idx2 + integer :: idx1, idx2, idx3 real*8 :: wrate,yr,yp - real*8 :: r1_(2,xx,yy,zz),f_(2,xx,yy,zz) + real*8 :: r1_(neq,xx,yy,zz),f_(neq,xx,yy,zz) real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) - real*8 :: ux(4,xx),dux(4,xx),d2ux(xx) - real*8 :: uy(4,yy),duy(4,yy),d2uy(yy) - real*8 :: uz(4,zz),duz(4,zz),d2uz(zz) - real*8 :: uux(xx) - real*8 :: uuy(yy) - real*8 :: uuz(zz) + real*8 :: ux(nd1,xx),dux(nd1,xx),d2ux(nd2,xx) + real*8 :: uy(nd1,yy),duy(nd1,yy),d2uy(nd2,yy) + real*8 :: uz(nd1,zz),duz(nd1,zz),d2uz(nd2,zz) + real*8 :: uux(nd2,xx) + real*8 :: uuy(nd2,yy) + real*8 :: uuz(nd2,zz) real*8 :: y @@ -485,17 +491,19 @@ module m_fdm_calc DO j=1,yy DO i=1,xx - y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y +! y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y - wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate - IF ((1.-y).le.c_ref) THEN - wrate=min_wr - IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ & - min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) - ENDIF + wrate=0.0 +! wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate +! IF ((1.-y).le.c_ref) THEN +! wrate=min_wr +! IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ & +! min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) +! ENDIF f_(1,i,j,k) = 0.0 ! continuity f_(2,i,j,k) = - wrate ! species conservation + f_(3,i,j,k) = - wrate ! species conservation ENDDO ENDDO @@ -509,14 +517,18 @@ module m_fdm_calc DO i=1,xx DO k=1,zz uz(1,k)=r1_(1,i,j,k) ! 1:rho - uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y + uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF uz(3,k)=uz(1,k)*ww_(i,j,k) ! 3:rho*w - uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y - uuz (k)=uz(2,k) + uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*YF + uz(5,k)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO + uz(6,k)=uz(3,k)*uz(5,k) ! 6:rho*w*YO + + uuz (1,k)=uz(2,k) + uuz (2,k)=uz(5,k) ENDDO - CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3) - CALL d2fp(yy,hy,uuz(:),d2uz(:),1,3) + CALL dfp (yy,hy, uz(1:nd1,:), duz(1:nd1,:),nd1,3) + CALL d2fp(yy,hy,uuz(1:nd2,:),d2uz(1:nd2,:),nd2,3) DO k=1,zz ! -( d(rho*w)/dz ) @@ -525,7 +537,8 @@ module m_fdm_calc ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz ) - f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv. + f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(1,k)+duz(1,k)*duz(2,k)) ! species conserv. + f_(3,i,j,k) = f_(3,i,j,k) - duz(6,k) + diff*(uz(1,k)*d2uz(2,k)+duz(1,k)*duz(5,k)) ! species conserv. ENDDO ENDDO ENDDO @@ -535,32 +548,39 @@ module m_fdm_calc DO k=1,zz DO j=1,yy DO i=1,xx - idx2 = xx*(j-1)+i - idx1 = (idx2-1)*4 + idx3 = xx*(j-1)+i + idx1 = (idx3-1)*nd1 + idx2 = (idx3-1)*nd2 fz(idx1+1,k) = r1_(1,i,j,k) ! 1:rho - fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y + fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF fz(idx1+3,k) = r1_(1,i,j,k)*ww_(i,j,k) ! 3:rho*w - fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! 4:rho*w*Y - fzz(idx2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) + fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! 4:rho*w*YF + fz(idx1+5,k) = r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO + fz(idx1+6,k) = r1_(3,i,j,k)*ww_(i,j,k) ! 6:rho*w*YO + + fzz(idx2+1,k) = r1_(2,i,j,k)/r1_(1,i,j,k) + fzz(idx2+2,k) = r1_(3,i,j,k)/r1_(1,i,j,k) ENDDO ENDDO ENDDO - CALL par_dfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3) - CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3) + CALL par_dfp (fz, fzl, fzu, dfz, hy, nd1*xx*yy, zz, yy, 3) + CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, nd2*xx*yy, zz, yy, 3) DO k=1,zz DO j=1,yy DO i=1,xx - idx2 = xx*(j-1)+i - idx1 = (idx2-1)*4 + idx3 = xx*(j-1)+i + idx1 = (idx3-1)*nd1 + idx2 = (idx3-1)*nd2 ! -( d(rho*w)/dz ) f_(1,i,j,k) = f_(1,i,j,k) - dfz(idx1+3,k) ! continuity ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz ) - f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv. + f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2+1,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv. + f_(3,i,j,k) = f_(3,i,j,k) - dfz(idx1+6,k) + diff*(fz(idx1+1,k)*dfzz(idx2+2,k)+dfz(idx1+1,k)*dfz(idx1+5,k)) ! species conserv. ENDDO ENDDO @@ -573,14 +593,18 @@ module m_fdm_calc DO j=1,yy DO i=1,xx ux(1,i)=r1_(1,i,j,k) ! 1:rho - ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y + ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF ux(3,i)=ux(1,i)*uu_(i,j,k) ! 3:rho*u - ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*Y - uux (i)=ux(2,i) + ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*YF + ux(5,i)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO + ux(6,i)=ux(3,i)*ux(5,i) ! 6:rho*u*YO + + uux (1,i)=ux(2,i) + uux (2,i)=ux(5,i) ENDDO - CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1) - CALL d2fnonp(xx,hx,uux(:),d2ux(:),1,1) + CALL dfnonp (xx,hx, ux(1:nd1,:), dux(1:nd1,:),nd1,1) + CALL d2fnonp(xx,hx,uux(1:nd2,:),d2ux(1:nd2,:),nd2,1) DO i=1,xx @@ -590,7 +614,8 @@ module m_fdm_calc ! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx ! = -( d(rho*u*Yr)/dx ) ! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx ) - f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) ! species conservation + f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(1,i)+dux(1,i)*dux(2,i)) ! species conservation + f_(3,i,j,k) = f_(3,i,j,k) - dux(6,i) + diff*(ux(1,i)*d2ux(2,i)+dux(1,i)*dux(5,i)) ! species conservation ENDDO ENDDO @@ -601,14 +626,18 @@ module m_fdm_calc DO i=1,xx DO j=1,yy uy(1,j)=r1_(1,i,j,k) ! 1:rho - uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y - uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v - uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*Y - uuy (j)=uy(2,j) + uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF + uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v + uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*YF + uy(5,j)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO + uy(6,j)=uy(3,j)*uy(5,j) ! 6:rho*v*YO + + uuy (1,j)=uy(2,j) + uuy (2,j)=uy(5,j) ENDDO - CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2) - CALL d2fp(yy,hy,uuy(:),d2uy(:),1,2) + CALL dfp (yy,hy, uy(1:nd1,:), duy(1:nd1,:),nd1,2) + CALL d2fp(yy,hy,uuy(1:nd2,:),d2uy(1:nd2,:),nd2,2) DO j=1,yy ! -( d(rho*v)/dy ) @@ -617,7 +646,8 @@ module m_fdm_calc ! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy ! = -( d(rho*v*Yr)/dy ) ! + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy ) - f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(j)+duy(1,j)*duy(2,j)) ! species conserv. + f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(1,j)+duy(1,j)*duy(2,j)) ! species conserv. + f_(3,i,j,k)=f_(3,i,j,k)-duy(6,j) + diff*(uy(1,j)*d2uy(2,j)+duy(1,j)*duy(5,j)) ! species conserv. ENDDO ENDDO ENDDO @@ -628,8 +658,10 @@ module m_fdm_calc DO j=1,yy DO i=1,yrsw f_(2,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr + f_(3,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr ENDDO IF (uu_(xx,j,k).lt.0.) f_(2,xx,j,k)=f_(1,xx,j,k)*r1_(2,xx,j,k)/r1_(1,xx,j,k) + IF (uu_(xx,j,k).lt.0.) f_(3,xx,j,k)=f_(1,xx,j,k)*r1_(3,xx,j,k)/r1_(1,xx,j,k) ENDDO ENDDO @@ -642,7 +674,7 @@ module m_fdm_calc integer :: istage,xx,yy,zz,i real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) - real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz) + real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz) istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_) @@ -660,7 +692,7 @@ module m_fdm_calc integer :: i,j,k,xx,yy,zz real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) - real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz) + real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz) CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf) @@ -674,7 +706,7 @@ module m_fdm_calc integer :: i,j,k,xx,yy,zz,istage real*8 :: at,bt , wrate , yr - real*8 :: ri(2,xx,yy,zz),r1(2,xx,yy,zz),r2(2,xx,yy,zz),f(2,xx,yy,zz) + real*8 :: ri(neq,xx,yy,zz),r1(neq,xx,yy,zz),r2(neq,xx,yy,zz),f(neq,xx,yy,zz) real*8 :: a(5),b(5) real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv @@ -699,7 +731,7 @@ module m_fdm_calc DO k=1,zz DO j=1,yy DO i=1,xx - DO nv=1,2 + DO nv=1,neq r1(nv,i,j,k)=r1(nv,i,j,k)+at*f(nv,i,j,k) r2(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k) ENDDO @@ -711,15 +743,17 @@ module m_fdm_calc DO k=1,zz DO j=1,yy DO i=1,xx - DO nv=1,2 + DO nv=1,neq r1(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k) ENDDO !==========rho=1 treatment r1(2,i,j,k)=r1(2,i,j,k)/r1(1,i,j,k) + r1(3,i,j,k)=r1(3,i,j,k)/r1(1,i,j,k) r1(1,i,j,k)=1. !==========Max Yr=1 treatment r1(2,i,j,k)=MIN(in_yr,r1(2,i,j,k)) + r1(3,i,j,k)=MIN(in_yr,r1(3,i,j,k)) !==========Min Yr=0 treatment ! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))