module m_fdm_calc use m_parameters use m_io use m_compact implicit none !variables real*8, dimension(:,:,:), allocatable :: u_,v_,w_ real*8, dimension(:,:,:), allocatable :: dm real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf real*8, dimension(:,:), allocatable :: fzzl, fzzu real*8, dimension(:,:,:), allocatable :: fbuf1, fbuf2, fbuf3, fbuf4 real*8, dimension(:,:), allocatable :: yxbuf1, yxbuf2, yxbuf3, yxbuf4 real*8, dimension(:,:), allocatable :: xybuf1, xybuf2, xybuf3, xybuf4 real*8 :: in_yr,out_yr,refwr,minf integer :: fullsavenum !,svfx,svfy integer :: fdmcyc,fdmsavecount real*8 :: fdmtime,fdmdt real*8 :: sumc,oldsumc,time_int,sum_wrate real*8 :: sumc1,sum_wrate1 real*8 :: visdt, convdt integer :: fdmstep real*8 :: coe,tt1,tt2,tt3,tt4 real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon real*8, dimension(3) :: velmax, velmin, velmax1, velmin1 integer, parameter :: neq = 2 logical :: fors !=========================================================================== !=========================================================================== contains subroutine fdm_restart_read integer(kind=MPI_INTEGER_KIND) :: fh integer(kind=MPI_OFFSET_KIND) :: offset write(fname,"('fort.',i4)") restartnum if (myid.eq.master) then write (6,*) '********************************************************' write (6,*) ' FDM field is being initialized for restarting.' write (6,457) restartnum 457 format(' Restart file number : ',i6) OPEN (restartnum,file=fname,form='unformatted',status='unknown') read (restartnum) fdmtime,tt1,tt2,tt3,oldsumc,time_int read (restartnum) fdmcyc,fdmdt,dummyu_ read (restartnum) tt1,tt2 !dt_fdmsave,dt_fullsave read (restartnum) t_fdmsave,t_fullsave read (restartnum) in_yr,out_yr CLOSE (restartnum) if(fdmtime.ge.t_fdmsave) t_fdmsave=t_fdmsave+dt_fdmsave if(fdmtime.ge.t_fullsave) t_fullsave=t_fullsave+dt_fullsave write (6,456) fdmtime,fdmcyc 456 format(' Restart time : ',f10.5,' / FDM cycle : ',i6) write(*,*) write(*,*) ' Save Times : FDM FULL ' write(*,454) t_fdmsave,t_fullsave 454 format(' ',f7.3,' ',f7.3) write(*,*)' Save Intervals : FDM FULL','from input' write(*,453) dt_fdmsave,dt_fullsave 453 format(' ',f7.3,' ',f7.3) write(*,*) if (ABS(fdmtime-TIME).le.1.0e-10) then write(*,*) ' Spectral and FDM times are consistent.' else write(*,*)' !! Warning : Spectral and FDM times are different !!' write(*,455) TIME,fdmtime 455 format(' !! Spectral Time : ',f10.5,' / FDM Time : ',f10.5) endif end if count = 1 call MPI_BCAST(fdmtime, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(tt1, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(tt2, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(tt3, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(oldsumc, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(time_int, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(fdmcyc, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(fdmdt, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(dummyu_, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) !call MPI_BCAST(dt_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) !call MPI_BCAST(dt_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) !call MPI_BCAST(t_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) !call MPI_BCAST(t_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(in_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) call MPI_BCAST(out_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) ! opening the file call MPI_INFO_CREATE(mpi_info, mpi_err) call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_RDONLY,mpi_info,fh,mpi_err) count = nx * ny * nz offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_READ_AT_ALL(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_READ_AT_ALL(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_READ_AT_ALL(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_READ_AT_ALL(fh, offset, y1(1,1,1,1), count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 4*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_READ_AT_ALL(fh, offset, y1(1,1,1,2), count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_CLOSE(fh, mpi_err) call MPI_INFO_FREE(mpi_info, mpi_err) end subroutine fdm_restart_read subroutine fdm_restart_write integer(kind=MPI_INTEGER_KIND) :: fh integer(kind=MPI_OFFSET_KIND) :: offset ! --------------- writing process ------------------ if (myid.eq.0) write(*,*) '=======================================================' if (myid.eq.0) write(*,*) 'Full results are being written',fullsavenum !fname = 'fort.'//savenum write(fname,"('fort.',i4)") fullsavenum ! opening the file call MPI_INFO_CREATE(mpi_info, mpi_err) call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_WRONLY+MPI_MODE_CREATE,mpi_info,fh,mpi_err) ! the master node writes the header with parameters if (myid.eq.0) then count = 1 call MPI_FILE_WRITE(fh, 8*6, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, fdmtime, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, nx, count, MPI_INTEGER8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, ny, count, MPI_INTEGER8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, nz_all, count, MPI_INTEGER8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, oldsumc, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, time_int, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*6, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*3, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, fdmcyc, count, MPI_INTEGER8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, DT, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, dummyu_, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*3, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, dt_fdmsave, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, dt_fullsave, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, t_fdmsave, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, t_fullsave, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, in_yr, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, out_yr, count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*2, count, MPI_INTEGER4, mpi_status, mpi_err) call MPI_FILE_WRITE(fh, 8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 5*nx*ny*nz_all*8 call MPI_FILE_WRITE_AT(fh, offset, 8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err) end if count = nx * ny * nz offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_WRITE_AT(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_WRITE_AT(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_WRITE_AT(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_WRITE_AT(fh, offset, y1(1,1,1,1), count, MPI_REAL8, mpi_status, mpi_err) offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 4*nx*ny*nz_all*8 + myid*nx*ny*nz*8 call MPI_FILE_WRITE_AT(fh, offset, y1(1,1,1,2), count, MPI_REAL8, mpi_status, mpi_err) call MPI_FILE_CLOSE(fh, mpi_err) call MPI_INFO_FREE(mpi_info, mpi_err) call MPI_Barrier(MPI_COMM_TASK, mpi_err) if (myid.eq.0) write(*,*) '=======================================================' fullsavenum=fullsavenum+1 end subroutine fdm_restart_write subroutine prepare_fdm implicit none integer :: i,j,ii,k real*8 :: fl_location,wrate,yr,c real*8 :: tar_sumc real*8 :: yy(nx) if (fdm_sw .eq. 0) then return endif fors = .true. allocate(u_(nx,ny,nz)) allocate(v_(nx,ny,nz)) allocate(w_(nx,ny,nz)) u_=0.0 v_=0.0 w_=0.0 allocate(dm(nx,ny,nz)) dm = diff ! DQ initializing fdmsavecount=1 !FDM save count sum_wrate=0. sum_wrate1=0. sumc=0. !for adjusting mean velocity sumc1=0. !for adjusting mean velocity oldsumc=0. time_int=0. visdt=9999. convdt=9999. fdmcyc=0 fdmtime=0. fullsavenum=1000 !full save file allocate(y1(nx,ny,nz,neq)) allocate(y2(nx,ny,nz,neq)) allocate(yf(nx,ny,nz,neq)) allocate(fzzu(nx*ny,2)) allocate(fzzl(nx*ny,2)) allocate(fbuf1(nx,ny,nz)) allocate(fbuf2(nx,ny,nz)) allocate(fbuf3(nx,ny,nz)) allocate(fbuf4(nx,ny,nz)) allocate(xybuf1(nx,ny)) allocate(xybuf2(nx,ny)) allocate(xybuf3(nx,ny)) allocate(xybuf4(nx,ny)) allocate(yxbuf1(ny,nx)) allocate(yxbuf2(ny,nx)) allocate(yxbuf3(ny,nx)) allocate(yxbuf4(ny,nx)) y1=0.0 y2=0.0 yf=0.0 CALL ludcmp(nx,ny,nz_all,1,0,0) refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon !FDM normal start================================== if (restartnum==0) then if (myid.eq.master) then OPEN(305,FILE='sfield.bin',form='unformatted',status='unknown') DO i=1,nx READ (305) yy(i) ! Yr ENDDO CLOSE (305) end if call MPI_BCAST(yy, nx, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err) in_yr=yy(1) ! inlet_Yr out_yr=yy(nx) ! outlet_Yr do j=1,nz do i=1,ny do ii=1,nx y1(ii,i,j,1)=1. ! rho initializing y1(ii,i,j,2)=yy(ii) ! Yr initializing enddo enddo enddo ! flame location setup sum_wrate=0.; sumc=0. DO ii=1,nz DO j=1,ny DO i=1,nx yr=y1(i,j,ii,2)/y1(i,j,ii,1) c=1.-yr IF (c.lt.0.) c=0. wrate=pre*yr*exp(-ac/(1.+bc*c)) ! cold boundary difficulty treatment IF(c.le.c_ref) THEN wrate=min_wr IF(c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ & min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) ENDIF sum_wrate=sum_wrate+wrate*(hx*hy*hy) sumc=sumc+(1.-yr) ENDDO ENDDO ENDDO sum_wrate1=sum_wrate/(hy*hy*ny*nz_all) CALL MPI_ALLREDUCE(sum_wrate1,sum_wrate,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) if (myid.eq.0) write(*,633) sum_wrate CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) fl_location=(hx*(nx-1.))*(1.-(sumc1/(nx*ny*nz_all))) sumc=sumc1*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.)) !CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) oldsumc=sumc if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100. 633 format (' ** Consumption Speed, Sc = ',f7.4) 634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.') ! Restart setup ================================================================ else CALL fdm_restart_read fullsavenum=restartnum+1 endif ! End of Restart setup ======================================================== return end subroutine prepare_fdm subroutine fdm_exe implicit none integer :: i,j,k,ii real*8 :: fl_location,delu,wrate,yr real*8 :: c real*8 :: uvel,vvel,wvel ! J. Kwon fdmcyc=itime ! Mean velocty setup do k=1,nz do i=1,nx do j=1,ny u_(i,j,k)=u_(i,j,k)+dummyu_ enddo enddo enddo if (myid.eq.0) write(*,*) if (myid.eq.0) write(*,*) '==========================================================' if (myid.eq.0) WRITE(*,932) TIME,DT 932 format(' Spectral results at time = ',f10.5,', dT = ',f7.5) if (myid.eq.0) write(*,933) dummyu_ 933 format(' ** Mean U = ',f7.4) !---------------------------------------------------------------- ! FDM DNS BGN !---------------------------------------------------------------- velmax1(1)=maxval(u_) velmin1(1)=minval(u_) velmax1(2)=maxval(v_) velmin1(2)=minval(v_) velmax1(3)=maxval(w_) velmin1(3)=minval(w_) ! MPI_REDUCE [uvw](max|min) CALL MPI_ALLREDUCE(velmax1,velmax,3,MPI_REAL8,MPI_MAX,MPI_COMM_TASK,mpi_err) CALL MPI_ALLREDUCE(velmin1,velmin,3,MPI_REAL8,MPI_MIN,MPI_COMM_TASK,mpi_err) umax=velmax(1) umin=velmin(1) vmax=velmax(2) vmin=velmin(2) wmax=velmax(3) wmin=velmin(3) uvel = max(abs(umax), abs(umin)) vvel = max(abs(vmax), abs(vmin)) wvel = max(abs(wmax), abs(wmin)) umax = umax - dummyu_ umin = umin - dummyu_ visdt=max(0.,0.3*fdmcfl*hx**2./nu) convdt=max(0.,fdmcfl*hx/(uvel+vvel+wvel)) fdmdt=min(DT,visdt,convdt) DT=fdmdt call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf) fdmtime=time+fdmdt if (myid.eq.0) write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt sumc=0. sum_wrate=0. do k=1,nz do j=1,ny do i=1,nx yr=y1(i,j,k,2)/y1(i,j,k,1) wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) IF((1.-yr).le.c_ref) THEN wrate=min_wr IF((1.-yr).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-yr-c_ref))+ & min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) ENDIF !get sum_wrate sum_wrate=sum_wrate+wrate*(hx*hy*hy) sumc=sumc+(1.-yr) enddo enddo enddo ! Real time results for Sc and flame location. sum_wrate1=sum_wrate/(hy*hy*REAL(ny*nz_all)) CALL MPI_ALLREDUCE(sum_wrate1,sum_wrate,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) if (myid.eq.0) write(*,633) sum_wrate 633 format (' ** Consumption Speed, Sc = ',f7.4) CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) fl_location=(hx*REAL(nx-1))*(1.-(sumc1/(REAL(nx*ny*nz_all)))) sumc1=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.)) CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err) if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100. 634 format (' ** Flame Location = ',f7.3, ' % point of x-domain.') if (myid.eq.0) write(St_data,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx-1)*hx)*100. ! Control the inflow mean velocity, dummyu_ if (myid.eq.0) write(*,'(a20,2f12.8)') ' ** oldsumc, sumc = ',oldsumc,sumc if(swadtv.ne.0.and.mod((fdmcyc+1),swadtv).eq.0) then time_int=fdmtime-time_int delu=(sumc-oldsumc)/time_int if ((fdmtime).gt.startad.and.oldsumc.ne.0.) then dummyu_=dummyu_+delu endif time_int=fdmtime endif !---------------------------------------------------------------------------- ! FDM DNS END !---------------------------------------------------------------------------- if (myid.eq.0) write(*,*) 'Number of whole calculation = ',fdmcyc if (myid.eq.0) write(*,*) '==========================================================' if (myid.eq.0) write(*,*) if (myid.eq.0) call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin ) return end subroutine fdm_exe SUBROUTINE fns(r1_,f_,xx,yy,zz,uu_,vv_,ww_) implicit none integer, parameter :: nb = 4 integer :: i,j,k,xx,yy,zz,ii,jj,kk integer :: n integer :: idx1, idx2 real*8 :: wrate,yr,yp real*8 :: r1_(xx,yy,zz,neq),f_(xx,yy,zz,neq) real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: uz(4,zz),duz(4,zz),d2uz(zz) real*8 :: uuz(zz) real*8 :: y real*8 :: T0, T1, Ts, As f_(:,:,:,1) = 0.0 ! continuity ! diffusivity T0 = 1.0 T1 = (1.0 + bc) Ts = (dm_ratio*T1 - (T1**(3./2.))) / (T1**(3./2.) - dm_ratio) As = (T0 + Ts) DO k=1,zz DO j=1,yy DO i=1,xx y=(1.0 - r1_(i,j,k,2)/r1_(i,j,k,1)) * bc + 1.0 dm(i,j,k) = diff * As * sqrt(y) / (1. + Ts/y) ENDDO ENDDO ENDDO if (myid.eq.0) write(*,*) 'min(dm)',minval(dm),'max(dm)',maxval(dm) ! reaction source term DO k=1,zz DO j=1,yy DO i=1,xx y=r1_(i,j,k,2)/r1_(i,j,k,1) ! 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 f_(i,j,k,2) = - wrate ! species conservation ENDDO ENDDO ENDDO !! z-direction ! sequential run if (numprocs.eq.1) then DO j=1,yy DO i=1,xx DO k=1,zz uz(1,k)=r1_(i,j,k,1)*dm(i,j,k) ! 1:rho*D uz(2,k)=r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y 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) ENDDO CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3) CALL d2fp(yy,hy,uuz(:),d2uz(:),1,3) DO k=1,zz ! -( d(rho*w)/dz ) f_(i,j,k,1) = f_(i,j,k,1) - duz(3,k) ! continuity ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz ) f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + (uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv. ENDDO ENDDO ENDDO ! parallel run else ! -( d(rho*w)/dz ) fbuf1(:,:,:) = r1_(:,:,:,1)*ww_(:,:,:) ! rho*w CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3) f_(:,:,:,1) = f_(:,:,:,1) - fbuf2(:,:,:) ! continuity ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz ! = -( d(rho*w*Yr)/dz ) ! + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz ) fbuf1(:,:,:) = r1_(:,:,:,2)*ww_(:,:,:) ! rho*w*Y CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3) f_(:,:,:,2) = f_(:,:,:,2) - fbuf2(:,:,:) ! species conserv. fbuf1(:,:,:) = r1_(:,:,:,2)/r1_(:,:,:,1) ! Y CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3) CALL pd2fp(fbuf1, fzzl, fzzu, fbuf3, hy, xx*yy, zz, yy, 3) fbuf1(:,:,:) = r1_(:,:,:,1) * dm(:,:,:) ! rho*D CALL pdfp (fbuf1, fzzl, fzzu, fbuf4, hy, xx*yy, zz, yy, 3) fbuf2 = fbuf2 * fbuf4 ! dY/dz * d(rho*D)/dz fbuf1 = fbuf1 * fbuf3 + fbuf2 ! rho*D * d2(rho*D)/dz2 + ... f_(:,:,:,2) = f_(:,:,:,2) + fbuf1(:,:,:) ! species conserv. endif ! x-direction DO k=1,zz ! -( d(rho*u)/dx ) CALL tp2mul (yxbuf1, r1_(:,:,k,1), uu_(:,:,k), yy, xx) ! rho*u CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u) CALL tp2subasgn (f_(:,:,k,1), yxbuf2, xx, yy) ! continuity ! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx ! = -( d(rho*u*Yr)/dx ) ! + (rho*D * d2(Yr)/dx2 + d(rho*D)/dx * d(Yr)/dx ) CALL tp2mul (yxbuf1, r1_(:,:,k,2), uu_(:,:,k), yy, xx) ! rho*u*Y CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u*Y) CALL tp2subasgn (f_(:,:,k,2), yxbuf2, xx, yy) ! species conservation CALL tp2div (yxbuf1, r1_(:,:,k,2), r1_(:,:,k,1), yy, xx)! Y CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(Y) CALL d2fnonp(xx,hx,yxbuf1,yxbuf3,yy,1) ! d2/dx2(Y) CALL tp2mul (yxbuf1, r1_(:,:,k,1), dm(:,:,k), yy, xx) ! rho*D CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho*D) yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho*D) yxbuf1 = -(yxbuf2 + yxbuf1 * yxbuf3) ! -( ... + rho*D * d2/dx2(Y)) CALL tp2subasgn (f_(:,:,k,2), yxbuf1, xx, yy) ! species conservation ENDDO !! y-direction DO k=1,zz ! -( d(rho*v)/dy ) xybuf1(:,:)=r1_(:,:,k,1)*vv_(:,:,k) ! rho*v CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(rho*v) f_(:,:,k,1)=f_(:,:,k,1)-xybuf2(:,:) ! continuity ! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy ! = -( d(rho*v*Yr)/dy ) ! + (rho*D* d2(Yr)/dyy2 + d(rho*D)/dy * d(Yr)/dy ) xybuf1(:,:) = r1_(:,:,k,2)*vv_(:,:,k) ! rho*v*Y CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(rho*v*Y) f_(:,:,k,2) = f_(:,:,k,2) - xybuf2(:,:) ! species conserv. xybuf1(:,:) = r1_(:,:,k,2)/r1_(:,:,k,1) ! Y CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(Y) CALL d2fp(yy,hy,xybuf1,xybuf3,xx,2) ! d2/dy2(Y) xybuf1(:,:) = r1_(:,:,k,1) * dm(:,:,k) ! rho*D CALL dfp(yy,hy,xybuf1,xybuf4,xx,2) ! d/dy(rho*D) xybuf2 = xybuf2 * xybuf4 ! d/dy(Y) * d/dy(rho*D) xybuf1 = xybuf2 + xybuf1 * xybuf3 ! ... + rho*D * d2/dy2(Y) f_(:,:,k,2) = f_(:,:,k,2) + xybuf1(:,:) ! species conserv. ENDDO ! Boundary condition DO k=1,zz DO j=1,yy DO i=1,yrsw f_(i,j,k,2)=r1_(i,j,k,1)*0.+f_(i,j,k,1)*in_yr ENDDO IF (uu_(xx,j,k).lt.0.) f_(xx,j,k,2)=f_(xx,j,k,1)*r1_(xx,j,k,2)/r1_(xx,j,k,1) ENDDO ENDDO return END SUBROUTINE fns subroutine EE1(xx,yy,zz,uu_,vv_,ww_,yy1,rhs1) implicit none integer :: xx,yy,zz real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: yy1(xx,yy,zz,neq),rhs1(xx,yy,zz,neq) integer :: i, j, k CALL fns(yy1,rhs1,xx,yy,zz,uu_,vv_,ww_) yy1 = yy1 + fdmdt * rhs1 !==========rho=1 treatment yy1(:,:,:,2) = yy1(:,:,:,2)/yy1(:,:,:,1) yy1(:,:,:,1) = 1. DO k = 1,zz DO j = 1,yy DO i = 1,xx !==========Max Yr=1 treatment yy1(i,j,k,2)=MIN(in_yr,yy1(i,j,k,2)) !==========Min Yr=0 treatment ! yy1(i,j,k,2)=MAX(out_yr,yy1(i,j,k,2)) ENDDO ENDDO ENDDO return END SUBROUTINE EE1 subroutine AB2(xx,yy,zz,uu_,vv_,ww_,yy1,rhs0,rhs1) implicit none integer :: xx,yy,zz real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: yy1(xx,yy,zz,neq),rhs0(xx,yy,zz,neq),rhs1(xx,yy,zz,neq) integer :: i, j, k CALL fns(yy1,rhs1,xx,yy,zz,uu_,vv_,ww_) yy1 = yy1 + fdmdt * ( 1.5d0 * rhs1 - 0.5d0 * rhs0 ) rhs0 = rhs1 !==========rho=1 treatment yy1(:,:,:,2) = yy1(:,:,:,2)/yy1(:,:,:,1) yy1(:,:,:,1) = 1. DO k = 1,zz DO j = 1,yy DO i = 1,xx !==========Max Yr=1 treatment yy1(i,j,k,2)=MIN(in_yr,yy1(i,j,k,2)) !==========Min Yr=0 treatment ! yy1(i,j,k,2)=MAX(out_yr,yy1(i,j,k,2)) ENDDO ENDDO ENDDO return END SUBROUTINE AB2 subroutine RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf) implicit none integer :: istage,xx,yy,zz,i real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq) CALL rotarysubstep(1, xx, yy, zz, uu_, vv_, ww_, yy1, yy1, yy2) CALL rotarysubstep(2, xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf) CALL rotarysubstep(3, xx, yy, zz, uu_, vv_, ww_, yy2, yyf, yy1) CALL rotarysubstep(4, xx, yy, zz, uu_, vv_, ww_, yyf, yy1, yy2) CALL lastsubstep (xx, yy, zz, uu_, vv_, ww_, yy1, yy2, yyf) return END SUBROUTINE RK4 SUBROUTINE solve(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf) IMPLICIT NONE integer :: i,j,k,xx,yy,zz real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq) ! advance reacting scalars - either Euler or Adams-Bashforth if (fors) then call EE1(xx,yy,zz,uu_,vv_,ww_,yy1,yy2) fors = .false. else call AB2(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf) yy2 = yyf end if return END SUBROUTINE solve subroutine rotarysubstep(istage,xx,yy,zz,uu_,vv_,ww_,r0,r1,r2) implicit none integer :: xx,yy,zz,istage real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz) real*8 :: r0(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq) real*8 :: a(5),b(5) real*8 :: at, bt integer :: i,j,k,l a(1)= 970286171893./4311952581923. a(2)= 6584761158862./12103376702013. a(3)= 2251764453980./15575788980749. a(4)= 26877169314380./34165994151039. a(5)=0. b(1)= 1153189308089./22510343858157. b(2)= 1772645290293./4653164025191. b(3)= -1672844663538./4480602732383. b(4)= 2114624349019./3568978502595. b(5)= 5198255086312./14908931495163. CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_) at=a(istage)*fdmdt bt=(b(istage) - a(istage))*fdmdt DO l = 1,neq DO k = 1,zz DO j = 1,yy DO i = 1,xx r1(i,j,k,l) = r1(i,j,k,l) + at*r2(i,j,k,l) r2(i,j,k,l) = r1(i,j,k,l) + bt*r2(i,j,k,l) ENDDO ENDDO ENDDO ENDDO return END SUBROUTINE rotarysubstep subroutine lastsubstep(xx,yy,zz,uu_,vv_,ww_,r0,r1,r2) implicit none integer :: xx, yy, zz real*8 :: bt real*8 :: r0(xx,yy,zz,neq), r1(xx,yy,zz,neq), r2(xx,yy,zz,neq) real*8 :: uu_(xx,yy,zz), vv_(xx,yy,zz), ww_(xx,yy,zz) real*8 :: b integer :: i, j, k b = 5198255086312./14908931495163. CALL fns(r0,r2,xx,yy,zz,uu_,vv_,ww_) bt = b*fdmdt r0 = r1 + bt*r2 !==========rho=1 treatment r0(:,:,:,2) = r0(:,:,:,2)/r0(:,:,:,1) r0(:,:,:,1) = 1. DO k = 1,zz DO j = 1,yy DO i = 1,xx !==========Max Yr=1 treatment r0(i,j,k,2)=MIN(in_yr,r0(i,j,k,2)) !==========Min Yr=0 treatment ! r0(i,j,k,2)=MAX(out_yr,r0(i,j,k,2)) ENDDO ENDDO ENDDO return END SUBROUTINE lastsubstep subroutine substep(ri,r1,r2,f,xx,yy,zz,istage,uu_,vv_,ww_) implicit none integer :: i,j,k,xx,yy,zz,istage real*8 :: at,bt , wrate , yr real*8 :: ri(xx,yy,zz,neq),r1(xx,yy,zz,neq),r2(xx,yy,zz,neq),f(xx,yy,zz,neq) 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 a(1)= 970286171893./4311952581923. a(2)= 6584761158862./12103376702013. a(3)= 2251764453980./15575788980749. a(4)= 26877169314380./34165994151039. a(5)=0. b(1)= 1153189308089./22510343858157. b(2)= 1772645290293./4653164025191. b(3)= -1672844663538./4480602732383. b(4)= 2114624349019./3568978502595. b(5)= 5198255086312./14908931495163. CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_) IF(istage<5) THEN at=a(istage)*fdmdt bt=b(istage)*fdmdt r2 = r1 + bt*f r1 = r1 + at*f ELSE bt=b(istage)*fdmdt r1 = r1 + bt*f !==========rho=1 treatment r1(:,:,:,2) = r1(:,:,:,2)/r1(:,:,:,1) r1(:,:,:,1) = 1. DO k = 1,zz DO j = 1,yy DO i = 1,xx !==========Max Yr=1 treatment r1(i,j,k,2)=MIN(in_yr,r1(i,j,k,2)) !==========Min Yr=0 treatment ! r1(i,j,k,2)=MAX(out_yr,r1(i,j,k,2)) ENDDO ENDDO ENDDO ENDIF return END SUBROUTINE substep subroutine tp2 (a, b, n1, n2) ! a = transpose(b) implicit none integer,intent(in) :: n1, n2 real*8,intent(out) :: a(n1,n2) real*8,intent(in) :: b(n2,n1) integer :: i,j,ii,jj integer,parameter :: nb = 16 DO jj=1,n2,nb DO ii=1,n1,nb DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 a(i,j) = b(j,i) ENDDO ENDDO ENDDO ENDDO end subroutine tp2 subroutine tp2div (a, b, c, n1, n2) ! a = transpose(b/c) implicit none integer,intent(in) :: n1, n2 real*8,intent(in) :: b(n2,n1), c(n2,n1) real*8,intent(out) :: a(n1,n2) integer :: i,j,ii,jj integer,parameter :: nb = 16 DO jj=1,n2,nb DO ii=1,n1,nb DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 a(i,j) = b(j,i) / c(j,i) ENDDO ENDDO ENDDO ENDDO end subroutine tp2div subroutine tp2mul (a, b, c, n1, n2) ! a = transpose(b*c) implicit none integer,intent(in) :: n1, n2 real*8,intent(in) :: b(n2,n1), c(n2,n1) real*8,intent(out) :: a(n1,n2) integer :: i,j,ii,jj integer,parameter :: nb = 16 DO jj=1,n2,nb DO ii=1,n1,nb DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 a(i,j) = b(j,i) * c(j,i) ENDDO ENDDO ENDDO ENDDO end subroutine tp2mul subroutine tp2subasgn (a, b, n1, n2) ! a = a - transpose(b) implicit none integer,intent(in) :: n1, n2 real*8,intent(inout) :: a(n1,n2) real*8,intent(in) :: b(n2,n1) integer :: i,j,ii,jj integer,parameter :: nb = 16 DO jj=1,n2,nb DO ii=1,n1,nb DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 a(i,j) = a(i,j) - b(j,i) ENDDO ENDDO ENDDO ENDDO end subroutine tp2subasgn end module m_fdm_calc