1034 lines
30 KiB
Fortran
1034 lines
30 KiB
Fortran
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 :: 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)
|
|
|
|
count = 2 * nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*2*nx*ny*nz*8
|
|
call MPI_FILE_READ_AT_ALL(fh, offset, y1, 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)
|
|
|
|
count = 2 * nx * ny * nz
|
|
offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 3*nx*ny*nz_all*8 + myid*2*nx*ny*nz*8
|
|
call MPI_FILE_WRITE_AT(fh, offset, y1, count, MPI_REAL8, mpi_status, mpi_err)
|
|
|
|
|
|
call MPI_FILE_CLOSE(fh, mpi_err)
|
|
call MPI_INFO_FREE(mpi_info, 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
|
|
|
|
! 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
|
|
|
|
|
|
f_(:,:,:,1) = 0.0 ! continuity
|
|
|
|
! 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) ! 1:rho
|
|
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 )
|
|
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
|
|
f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + diff*(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 )
|
|
! + D* (rho* d2(Yr)/dz2 + d(rho)/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) ! rho
|
|
CALL pdfp (fbuf1, fzzl, fzzu, fbuf4, hy, xx*yy, zz, yy, 3)
|
|
|
|
fbuf2 = fbuf2 * fbuf4
|
|
fbuf1 = fbuf1 * fbuf3 + fbuf2
|
|
|
|
f_(:,:,:,2) = f_(:,:,:,2) + diff*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 )
|
|
! + D* (rho* d2(Yr)/dx2 + d(rho)/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 tp2 (yxbuf1, r1_(:,:,k,1), yy, xx) ! rho
|
|
CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho)
|
|
|
|
yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho)
|
|
yxbuf1 = -diff*(yxbuf2 + yxbuf1 * yxbuf3) ! -D( ... + rho * 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 )
|
|
! + D* (rho* d2(Yr)/dyy2 + d(rho)/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) ! rho
|
|
CALL dfp(yy,hy,xybuf1,xybuf4,xx,2) ! d/dy(rho)
|
|
|
|
xybuf2 = xybuf2 * xybuf4 ! d/dy(Y) * d/dy(rho)
|
|
xybuf1 = xybuf2 + xybuf1 * xybuf3 ! ... + rho * d2/dy2(Y)
|
|
|
|
f_(:,:,k,2) = f_(:,:,k,2) + diff*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
|
|
|
|
|