research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Compare commits

base: research:mpi_fdm
Branches Tags
research:mpi_fdm
research:two-step
research:more-iterations
research:diffusivity
research:test_branch
research:pipeline-compact
research:ppc
..
compare: research:more-iterations
Branches Tags
research:two-step
research:more-iterations
research:diffusivity
research:test_branch
research:pipeline-compact
research:ppc
research:mpi_fdm

34 commits
mpi_fdm ... more-itera

Author SHA1 Message Date
ignis
95d915b4c6 file_ext extended from 6-digit to 8-digit 2019-04-06 23:19:33 +09:00
ignis
0849b31e7b diffusivity calculation using sutherland formula 2019-04-04 22:16:37 +09:00
ignis
d798fe343e added mpi barrier after restat write to avoid send/recv time out 2019-04-04 17:08:52 +09:00
ignis
8431f9cac8 fix bug rho => rho*D 2019-03-25 04:41:23 +09:00
ignis
2934fb0e44 T^0.76 diffusivity 2019-03-24 13:59:40 +09:00
ignis
f38b995827 governing equation changed to consider variable diffusivity 2019-03-24 10:48:55 +09:00
Dohyun Kim
0776f3ef6e fdm_restart write/read update according to solution array layout change 2018-03-20 21:29:06 +09:00
ignis
e16ead959d Adams-Bashforth method for reacting scalar time advancement 2018-01-21 06:24:41 +09:00
ignis
f5b0550436 pipeline td solve, distinct tag for each rhs 2017-12-24 20:42:30 +09:00
ignis
4a8a67d241 fixed bug in ialias array construction 2017-12-24 20:19:22 +09:00
ignis
db7732025a remove unnecessary lines in gather_4 and added comments in main 2017-12-24 20:18:40 +09:00
ignis
a56c4d978c remove loops in parallel z-dir fns 2017-11-11 00:43:27 +09:00
ignis
e659052e26 z-dir rhs exchange mpi subroutine change, rk substep advance opt 2017-11-10 23:39:52 +09:00
ignis
b426be89ad transpose subroutine convention change and minor mods 2017-11-10 16:37:35 +09:00
ignis
749e77de32 corse grained and term-wise derivative computation 2017-11-09 03:42:30 +09:00
ignis
6933297130 remove unnecessary ifs and optimize parallel td solve 2017-11-09 03:11:09 +09:00
ignis
75f0feaa02 Makefile update for optimization flags 2017-11-09 03:05:23 +09:00
ignis
27bf7ab6dd Makefile added option for profiling and solution array axis roll 2017-11-07 03:52:47 +09:00
ignis
3ee83a9fff Temperton circulant tridiagonal algorithm 2017-11-04 01:36:43 +09:00
ignis
56fc50e9f0 par_ptdslv performance tune 2017-11-03 22:29:37 +09:00
ignis
b39e4cbef4 Makefile fftw cluster edition 2017-11-03 03:35:07 +09:00
ignis
dac299c060 Makefile update to use gfortran compiler 2017-08-02 19:39:19 +09:00
ignis
3225cd757b fully mpi working 2017-02-06 12:49:15 +09:00
ignis
6c75483f1e pipeline parallel compact scheme 2017-01-24 15:50:20 +09:00
ignis
383b2c5aa7 init_velocity seperate n_shell range test and e_spec value test 2017-01-24 15:50:20 +09:00
ignis
6f0c1f44fa fdm restart for mpi 2017-01-24 15:50:20 +09:00
ignis
4f60757a9a mpi communication established 2017-01-24 15:50:20 +09:00
ignis
016bd9fb80 gather4 only performs inverse transform 2017-01-24 15:50:20 +09:00
ignis
804698f717 every process call prepare_fdm and fdm_exe, do not work correctly 2017-01-24 15:50:20 +09:00
ignis
4f7bb89461 remove omp 2017-01-24 15:50:20 +09:00
ignis
50bc5073c1 remove old HIT3D save after writing fdm save 2017-01-24 15:50:20 +09:00
ignis
9c48452b0f disable io 2017-01-24 15:50:20 +09:00
ignis
6d7b31c555 write tp_field flush after each x loop 2017-01-24 15:03:34 +09:00
CombustionLab.POSTECH
cfe19ffb17 job_runlimit default to one month 2017-01-24 14:09:18 +09:00
16 changed files with 1401 additions and 303 deletions
Show stats Expand all files Collapse all files

86
Makefile
View file

@ -4,77 +4,37 @@
MPIF90 = blah
ifeq ($(OMPI), 1)
FFTW_HOME = ../libs/fftw-3.2.2/lib
OPENMPI_HOME = /home/comb/ignis/openmpi/
MPIF90 = $(OPENMPI_HOME)/bin/mpif90
FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include -I$(OPENMPI_HOME)/include
LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -L$(OPENMPI_HOME)/lib -lfftw3 -lm -shared-intel
FCFLAGS_F77 = -Mextend
FFTW_HOME = $(FFTW_DIR)
MPIF90 = mpif90
FCFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 $(MPIPP_LD_FLAGS) $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 =
MODULES += m_openmpi.o
else
FFTW_HOME = ../libs/fftw-3.2.2/lib
# compile for mpiP MPI Profiling
ifeq ($(PROFILER), mpiP)
FFTW_HOME = $(FFTW_DIR)
MPIF90 = mpif90
FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 = -Mextend
MODULES += m_mpi.o
endif
# Franklin cluster, NERSC
ifeq ($(NERSC_HOST), franklin)
MPIF90 = ftn
FCFLAGS = -target=linux -i8 -r8 -O4 -c
LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
FCFLAGS_F77 = -Mextend
endif
# Hopper cluster, NERSC
ifeq ($(NERSC_HOST), hopper)
MPIF90 = ftn
FCFLAGS = -target=linux -i8 -r8 -O4 -c
LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
FCFLAGS_F77 = -Mextend
endif
# Yellowrail cluster, LANL
ifeq ($(HOSTNAME), yr-fe1.lanl.gov)
MPIF90 = mpif90
FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
endif
# Coyote cluster, LANL
ifeq ($(HOSTNAME), cy-c2.lanl.gov)
MPIF90 = mpif90
FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
endif
# Linux Cluster, CMU
ifeq ($(HOSTNAME), karman.me.cmu.edu)
FFTW_HOME = ../libs/fftw-3.2.2/lib
MPIF90 = mpif90
FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 = -Mextend
FCFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
MPIP_LD_FLAGS = -L$(MPIP_DIR)/lib -lmpiP -lm -lbfd -liberty -lunwind
LDFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 $(MPIPP_LD_FLAGS) $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 =
endif
# sparrow.stanford.edu, Mac OS X box
ifeq ($(HOSTNAME), sparrow.stanford.edu)
MPIF90 = mpif90
FCFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -c
FCFLAGS_F77 = -ffixed-form -ffixed-line-length-none
FCFLAGS_F90 = -ffree-form -ffree-line-length-none
LDFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -lmpi -lmpi -lfftw3 -lm
# tau source instrumentation
ifeq ($(PROFILER), tau_source)
FFTW_HOME = $(FFTW_DIR)
MPIF90 = tau_f90.sh
FCFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -Ofast -march=native -Wall -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
FCFLAGS_F77 =
endif
# Program name
PROG = hit3d.x

2
begin_new.f90
View file

@ -15,7 +15,7 @@ subroutine begin_new
! defining the iteration number
ITIME = 0
file_ext = '000000'
file_ext = '00000000'

42
gather_4.f90
View file

@ -23,8 +23,6 @@ subroutine gather_4
! number of variables to write out
n_out = 3
if (int_scalars) n_out = n_out + n_scalars
if (les .and. n_les>0) n_out = n_out + n_les
! putting all variables in wrk array
do k = 1,nz
@ -44,47 +42,13 @@ subroutine gather_4
! velocities
call xFFT3d(-1,1)
fname = 'u.'//file_ext
u_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
call xFFT3d(-1,2)
fname = 'v.'//file_ext
v_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
call xFFT3d(-1,3)
fname = 'w.'//file_ext
u_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
v_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
w_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
! scalars
if (int_scalars) then
do n = 1,n_scalars
call xFFT3d(-1,3+n)
write(fname,"('sc',i2.2,'.',a6)") n,file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n)
call write_tmp4
end do
end if
! LES quantities
if (les) then
! turbulent viscosity
if (allocated(turb_visc)) then
write(fname,"('nu_t.',a6)") file_ext
tmp4 = turb_visc
call write_tmp4
end if
if (n_les > 0) then
do n = 1, n_les
call xFFT3d(-1,3+n_scalars+n)
write(fname,"('les',i1,'.',a6)") n,file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n_scalars+n)
call write_tmp4
end do
end if
end if
return
end subroutine gather_4

2
get_file_ext.f90
View file

@ -5,7 +5,7 @@
use m_io, only : file_ext
implicit none
write(file_ext,"(i6.6)") itime
write(file_ext,"(i8.8)") itime
return
end subroutine get_file_ext

4
init_velocity.f90
View file

@ -266,8 +266,10 @@ subroutine init_velocity
do i = 1,nx+2
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kmax .and. e_spec(n_shell) .gt. zip) then
if (n_shell .gt. 0 .and. n_shell .le. kmax ) then
if (e_spec(n_shell) .gt. zip) then
fields(i,j,k,1:3) = fields(i,j,k,1:3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
endif
else
fields(i,j,k,1:3) = zip
end if

6
io_write_4.f90
View file

@ -75,7 +75,7 @@ subroutine io_write_4
if (int_scalars) then
do n = 1,n_scalars
call xFFT3d(-1,3+n)
write(fname,"('sc',i2.2,'.',a6)") n,file_ext
write(fname,"('sc',i2.2,'.',a8)") n,file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n)
call write_tmp4
@ -86,7 +86,7 @@ subroutine io_write_4
if (les) then
! turbulent viscosity
if (allocated(turb_visc)) then
write(fname,"('nu_t.',a6)") file_ext
write(fname,"('nu_t.',a8)") file_ext
tmp4 = turb_visc
call write_tmp4
end if
@ -94,7 +94,7 @@ subroutine io_write_4
if (n_les > 0) then
do n = 1, n_les
call xFFT3d(-1,3+n_scalars+n)
write(fname,"('les',i1,'.',a6)") n,file_ext
write(fname,"('les',i1,'.',a8)") n,file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3+n_scalars+n)
call write_tmp4
end do

852
m_compact.f90
View file

@ -1,14 +1,20 @@
MODULE m_compact
use m_openmpi
IMPLICIT NONE
PRIVATE
REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs
! lyzf,lyzs,wyzf,wyzs
REAL*8, DIMENSION(:), ALLOCATABLE :: zz1, zz2, vz1, vz2
INTEGER :: nxc,nyc,nzc
REAL*8, PARAMETER :: ezero = 1.0e-14
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp,pdfp,pd2fp
CONTAINS
@ -65,6 +71,22 @@
call nonp_lud(3,nzc)
ENDIF
ALLOCATE(vz1(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(zz1(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(vz2(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(zz2(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
call ctdlu (3.0, nzc, vz1, zz1)
call ctdlu (11.0/2.0, nzc, vz2, zz2)
! CALL x_lud
! CALL yz_lud
@ -141,6 +163,280 @@
l(n)=1./d
END SUBROUTINE ptdlu
SUBROUTINE ctdlu(lambda,n,v,z)
INTEGER, INTENT(IN) :: n
REAL*8, INTENT(IN) :: lambda
REAL*8, INTENT(OUT) :: v(n)
REAL*8, INTENT(OUT) :: z(n)
INTEGER :: i
REAL*8 :: sigma, alpha
v(1)=0.0
v(2)=lambda
DO i=3,n
v(i)=lambda - 1./v(i-1)
ENDDO
if (lambda > 0.0) then
alpha = (-lambda + sqrt(lambda**2 - 4.0))/2.0
else
alpha = (-lambda - sqrt(lambda**2 - 4.0))/2.0
endif
sigma = (1.0 + alpha**2) / (lambda * (1.0 - alpha**2) * (1.0 - alpha**n))
DO i=1,n
z(i) = sigma * (alpha**(i-1) + alpha**(n-i+1))
ENDDO
END SUBROUTINE ctdlu
SUBROUTINE pdfp(x,xu,xl,dx,h,nd,n,nall,dir)
INTEGER,INTENT(IN) :: nall,n,nd,dir
REAL*8,INTENT(IN) :: h
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL*8,INTENT(OUT),DIMENSION(nd,2) :: xu, xl
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL*8 :: r1,r2,h1
integer (kind=MPI_INTEGER_KIND) :: idx
integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
logical :: recvLow, recvUpp
integer (kind=MPI_INTEGER_KIND), parameter :: rhstag = 0
h1=1./h
r1=7./3.
r2=1./12.
!if (myid.eq.0) write(*,*) "parallel dfp commuication"
if (myid.eq.master) then
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, numprocs-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, numprocs-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, 0, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
else
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
endif
call MPI_STARTALL (4, requests, mpi_err)
DO i=3,n-2
DO j=1,nd
dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))*h1
ENDDO
ENDDO
recvLow = .true.
recvUpp = .true.
do while (recvLow.or.recvUpp)
call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)
if (idx.eq.1) then
DO j=1,nd
dx(j,1) =(r1*( x(j,2)-xl(j,2)) +r2*( x(j,3)-xl(j,1))) *h1
ENDDO
DO j=1,nd
dx(j,2) =(r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,2))) *h1
ENDDO
recvLow = .false.
endif
if (idx.eq.2) then
DO j=1,nd
dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
ENDDO
DO j=1,nd
dx(j,n) =(r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
ENDDO
recvUpp = .false.
endif
enddo
IF (dir.eq.1) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
IF (dir.eq.2) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
IF (dir.eq.3) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
END SUBROUTINE pdfp
SUBROUTINE pd2fp(x,xu,xl,dx,h,nd,n,nall,dir)
INTEGER,INTENT(IN) :: nall,n,nd,dir
REAL*8,INTENT(IN) :: h
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL*8 :: h2,r1,r2,t1,t2
REAL*8,DIMENSION(nd,2) :: xu, xl
integer (kind=MPI_INTEGER_KIND) :: idx
integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
logical :: recvLow, recvUpp
integer (kind=MPI_INTEGER_KIND), parameter :: rhstag = 0
h2=1./(h*h)
r1=6.
r2=3./8.
!if (myid.eq.0) write(*,*) "parallel d2fp commuication"
if (myid.eq.master) then
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, numprocs-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, numprocs-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, 0, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
else
call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(1), mpi_err)
call MPI_SEND_INIT (x(1,1), 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(2), mpi_err)
call MPI_RECV_INIT (xl, 2*nd, MPI_REAL8, myid-1, rhstag, &
MPI_COMM_TASK, requests(3), mpi_err)
call MPI_RECV_INIT (xu, 2*nd, MPI_REAL8, myid+1, rhstag, &
MPI_COMM_TASK, requests(4), mpi_err)
endif
call MPI_STARTALL (4, requests, mpi_err)
DO i=3,n-2
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
dx(j,i)=(r1*t1+r2*t2)*h2
ENDDO
ENDDO
recvLow = .true.
recvUpp = .true.
do while (recvLow.or.recvUpp)
call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)
if (idx.eq.1) then
DO j=1,nd
t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
dx(j,1)=(r1*t1+r2*t2)*h2
ENDDO
DO j=1,nd
t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
dx(j,2)=(r1*t1+r2*t2)*h2
ENDDO
recvLow = .false.
endif
if (idx.eq.2) then
DO j=1,nd
t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
dx(j,n-1)=(r1*t1+r2*t2)*h2
ENDDO
DO j=1,nd
t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
dx(j,n)=(r1*t1+r2*t2)*h2
ENDDO
recvUpp = .false.
endif
enddo
IF (dir.eq.1) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! x-direction
IF (dir.eq.2) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! y-direction
IF (dir.eq.3) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! z-direction
END SUBROUTINE pd2fp
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL*8,INTENT(IN) :: h
@ -219,6 +515,85 @@
END SUBROUTINE dfp
SUBROUTINE par_dfp(x,xu,xl,dx,h,nd,n,nall,dir)
INTEGER,INTENT(IN) :: nall,n,nd,dir
REAL*8,INTENT(IN) :: h
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL*8 :: r1,r2,h1
REAL*8,DIMENSION(nd,2) :: xu, xl
h1=1./h
r1=7./3.
r2=1./12.
!if (myid.eq.0) write(*,*) "parallel dfp commuication"
if (myid.eq.master) then
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
else
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
endif
DO j=1,nd
dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
dx(j,n)= (r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
dx(j,1)= (r1*( x(j,2)-xl(j,2)) +r2*( x(j,3)-xl(j,1)))*h1
dx(j,2)= (r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,2)))*h1
ENDDO
DO i=3,n-2
DO j=1,nd
dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))*h1
ENDDO
ENDDO
IF (dir.eq.1) CALL par_ptdslv(dx,lxf,wxf,nd,n,nall) ! x-direction
IF (dir.eq.2) CALL par_ptdslv(dx,lyf,wyf,nd,n,nall) ! x-direction
IF (dir.eq.3) CALL par_ptdslv(dx,lzf,wzf,nd,n,nall) ! x-direction
END SUBROUTINE par_dfp
SUBROUTINE ptdslv(r,n,l,w,nd)
INTEGER,INTENT(IN) :: n,nd
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
@ -247,6 +622,362 @@
ENDDO
END SUBROUTINE ptdslv
SUBROUTINE ctdslv(r,va,za,nd,n,nall)
INTEGER,PARAMETER :: nb = 256
INTEGER,PARAMETER :: nb2 = 64
INTEGER,INTENT(IN) :: n,nd,nall
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL*8,INTENT(IN),DIMENSION(nall) :: va
REAL*8,INTENT(IN),DIMENSION(nall) :: za
REAL*8,DIMENSION(n) :: v,z
INTEGER i,j
INTEGER ii,jj
REAL*8, DIMENSION(nd) :: x1, x0
INTEGER npart, nbase, nlow, nupp
INTEGER :: pid, np
REAL*8, DIMENSION(nb) :: r0, r1
pid = myid
np = numprocs
npart = nall / np
nbase = pid * npart
nlow = pid * npart + 1
nupp = (pid + 1) * npart
v = va(nlow:nupp)
z = za(nlow:nupp)
if (npart.lt.4) then
! assertion fail
endif
x1 = 0.0
x0 = 0.0
do i=1,n
x0 = x0 + z(i)*r(:,i)
enddo
call MPI_REDUCE(x0, x1, nd, MPI_REAL8, MPI_SUM, master, MPI_COMM_TASK, mpi_err)
! first process
if (myid.eq.0) then
CALL MPI_SSEND(x1, nd, MPI_REAL8, np-1, 100, MPI_COMM_TASK, mpi_err)
r(:,1) = x1
r(:,2) = r(:,2) - x1
! last process
elseif (pid.eq.(np-1)) then
CALL MPI_RECV(x1, nd, MPI_REAL8, master, 100, MPI_COMM_TASK, mpi_status, mpi_err)
r(:,n) = r(:,n) - x1
! intermediate process
else
endif
! first process
if (myid.eq.0) then
DO jj=1,nd,nb
DO j=jj,jj+nb-1
r(j,2)=r(j,2)/v(2)
ENDDO
DO i=3,n
DO j=jj,jj+nb-1
r(j,i)=(r(j,i)-r(j,i-1))/v(i)
ENDDO
ENDDO
r0 = r(jj:jj+nb-1,n)
CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, jj, MPI_COMM_TASK, mpi_err)
ENDDO
DO jj=1,nd,nb
CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, nd+jj, MPI_COMM_TASK, mpi_status, mpi_err)
i=n
DO j=jj,jj+nb-1
r(j,i) = r(j,i) - r1(j-jj+1)/v(i)
ENDDO
DO i=n-1,2,-1
DO j=jj,jj+nb-1
r(j,i) = r(j,i) - r(j,i+1)/v(i)
ENDDO
ENDDO
ENDDO
! last process
elseif (pid.eq.(np-1)) then
DO jj=1,nd,nb
CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, jj, MPI_COMM_TASK, mpi_status, mpi_err)
i=1
DO j=jj,jj+nb-1
r(j,i)=(r(j,i) - r1(j-jj+1))/v(i)
ENDDO
DO i=2,n
DO j=jj,jj+nb-1
r(j,i)=(r(j,i) - r(j,i-1))/v(i)
ENDDO
ENDDO
ENDDO
DO jj=1,nd,nb
DO i=n-1,1,-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r(j,i+1)/v(i)
ENDDO
ENDDO
r0 = r(jj:jj+nb-1,1)
CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, nd+jj, MPI_COMM_TASK, mpi_err)
ENDDO
! intermediate process
else
DO jj=1,nd,nb
CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, jj, MPI_COMM_TASK, mpi_status, mpi_err)
i=1
DO j=jj,jj+nb-1
r(j,i)=(r(j,i)-r1(j-jj+1))/v(i)
ENDDO
DO i=2,n
DO j=jj,jj+nb-1
r(j,i)=(r(j,i)-r(j,i-1))/v(i)
ENDDO
ENDDO
r0 = r(jj:jj+nb-1,n)
CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, jj, MPI_COMM_TASK, mpi_err)
ENDDO
DO jj=1,nd,nb
CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, nd+jj, MPI_COMM_TASK, mpi_status, mpi_err)
i=n
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r1(j-jj+1)/v(i)
ENDDO
DO i=n-1,1,-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r(j,i+1)/v(i)
ENDDO
ENDDO
r0 = r(jj:jj+nb-1,1)
CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, nd+jj, MPI_COMM_TASK, mpi_err)
ENDDO
endif
END SUBROUTINE ctdslv
SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall)
INTEGER,PARAMETER :: nb = 64
INTEGER,INTENT(IN) :: n,nd,nall
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL*8,INTENT(IN),DIMENSION(nall) :: la,wa
REAL*8,DIMENSION(n) :: l,w
INTEGER i,j
INTEGER ii,jj
REAL*8, DIMENSION(nd) :: sum
INTEGER npart, nbase, nlow, nupp
INTEGER :: pid, np
REAL*8, DIMENSION(nb) :: r0, r1, sum0
REAL*8, DIMENSION(nb,2) :: buf
pid = myid
np = numprocs
npart = nall / np
nbase = pid * npart
nlow = pid * npart + 1
nupp = (pid + 1) * npart
l = la(nlow:nupp)
w = wa(nlow:nupp)
if (npart.lt.4) then
! assertion fail
endif
call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
! first process
if (myid.eq.0) then
DO jj=1,nd,nb
DO j=jj,jj+nb-1
sum(j)=w(1)*r(j,1)
r(j,1)=r(j,1)*l(1)
ENDDO
DO i=2,n
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r(j,i-1)
sum(j)=sum(j)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
ENDDO
buf(:,1) = r(jj:jj+nb-1,n)
buf(:,2) = sum(jj:jj+nb-1)
CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
ENDDO
DO jj=1,nd,nb
CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
r0 = buf(:,1)
r1 = buf(:,2)
i=n
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
ENDDO
DO i=n-1,1,-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
ENDDO
ENDDO
ENDDO
! last process
elseif (pid.eq.(np-1)) then
DO jj=1,nd,nb
CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
r0 = buf(:,1)
sum0 = buf(:,2)
i=1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r0(j-jj+1)
sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
DO i=2,n-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r(j,i-1)
sum(j)=sum(j)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
ENDDO
DO j=jj,jj+nb-1
r(j,n)=l(n)*(r(j,n)-sum(j))
ENDDO
ENDDO
DO jj=1,nd,nb
DO j=jj,jj+nb-1
r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
ENDDO
DO i=n-2,1,-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
ENDDO
ENDDO
buf(:,1) = r(jj:jj+nb-1,1)
buf(:,2) = r(jj:jj+nb-1,n)
CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
ENDDO
! intermediate process
else
DO jj=1,nd,nb
CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
r0 = buf(:,1)
sum0 = buf(:,2)
i=1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r0(j-jj+1)
sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
DO i=2,n
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-r(j,i-1)
sum(j)=sum(j)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
ENDDO
buf(:,1) = r(jj:jj+nb-1,n)
buf(:,2) = sum(jj:jj+nb-1)
CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
ENDDO
DO jj=1,nd,nb
CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
r0 = buf(:,1)
r1 = buf(:,2)
i=n
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
ENDDO
DO i=n-1,1,-1
DO j=jj,jj+nb-1
r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
ENDDO
ENDDO
buf(:,1) = r(jj:jj+nb-1,1)
!buf(:,2) = r1
CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
ENDDO
endif
END SUBROUTINE par_ptdslv
SUBROUTINE d2fp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL*8,INTENT(IN) :: h
@ -260,30 +991,22 @@
DO j=1,nd
t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2))
t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3))
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
dx(j,n-1)=(r1*t1+r2*t2)
t1 = (x(j,1)-2.*x(j,n)+x(j,n-1))
t2 = (x(j,2)-2.*x(j,n)+x(j,n-2))
IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
dx(j,n)=(r1*t1+r2*t2)
t1 = (x(j,2)-2.*x(j,1)+x(j,n))
t2 = (x(j,3)-2.*x(j,1)+x(j,n-1))
IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0.
IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0.
! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
dx(j,1)=(r1*t1+r2*t2)
t1 = (x(j,3)-2.*x(j,2)+x(j,1))
t2 = (x(j,4)-2.*x(j,2)+x(j,n))
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0.
! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
! +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
dx(j,2)=(r1*t1+r2*t2)
@ -297,8 +1020,6 @@
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
dx(j,i)=(r1*t1+r2*t2)
@ -310,6 +1031,109 @@
IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
END SUBROUTINE d2fp
SUBROUTINE par_d2fp(x,xu,xl,dx,h,nd,n,nall,dir)
INTEGER,INTENT(IN) :: nall,n,nd,dir
REAL*8,INTENT(IN) :: h
REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL*8 :: h2,r1,r2,t1,t2
REAL*8,DIMENSION(nd,2) :: xu, xl
h2=1./(h*h)
r1=6.
r2=3./8.
!if (myid.eq.0) write(*,*) "parallel d2fp commuication"
if (myid.eq.master) then
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
elseif (myid.eq.numprocs-1) then
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
else
call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, &
MPI_COMM_TASK, mpi_status, mpi_err)
call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
MPI_COMM_TASK, mpi_request, mpi_err)
call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, &
MPI_COMM_TASK, mpi_status, mpi_err)
endif
DO j=1,nd
t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
dx(j,n-1)=(r1*t1+r2*t2)*h2
t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
dx(j,n)=(r1*t1+r2*t2)*h2
t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
dx(j,1)=(r1*t1+r2*t2)*h2
t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
! +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
dx(j,2)=(r1*t1+r2*t2)*h2
ENDDO
DO i=3,n-2
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
dx(j,i)=(r1*t1+r2*t2)*h2
ENDDO
ENDDO
IF (dir.eq.1) CALL par_ptdslv(dx,lxs,wxs,nd,n,nall) ! x-direction
IF (dir.eq.2) CALL par_ptdslv(dx,lys,wys,nd,n,nall) ! y-direction
IF (dir.eq.3) CALL par_ptdslv(dx,lzs,wzs,nd,n,nall) ! z-direction
END SUBROUTINE par_d2fp
SUBROUTINE tdslv(r,n,l,nd)
INTEGER,INTENT(IN) :: n,nd
REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
@ -354,10 +1178,6 @@
dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1))
dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2))
dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3))
IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0.
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0.
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0.
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0.
dx(j,1)=dx(j,1)*h2
dx(j,2)=dx(j,2)*h2*r3
dx(j,n-1)=dx(j,n-1)*h2*r3
@ -367,8 +1187,6 @@
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))

683
m_fdm_calc.f90
View file

@ -9,7 +9,13 @@ module m_fdm_calc
!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
@ -23,10 +29,9 @@ module m_fdm_calc
real*8 :: umax,umin,vmax,vmin,wmax,wmin ! J. Kwon
real*8, dimension(3) :: velmax, velmin, velmax1, velmin1
integer :: vtype1, vtype2
integer, dimension(128) :: scnt1, sdisp1, stype1, rcnt1, rdisp1, rtype1
integer, dimension(128) :: scnt2, sdisp2, stype2, rcnt2, rdisp2, rtype2
integer, parameter :: neq = 2
logical :: fors
!===========================================================================
!===========================================================================
@ -115,9 +120,11 @@ module m_fdm_calc
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)
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)
@ -190,14 +197,18 @@ module m_fdm_calc
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)
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
@ -215,6 +226,8 @@ module m_fdm_calc
return
endif
fors = .true.
allocate(u_(nx,ny,nz))
allocate(v_(nx,ny,nz))
allocate(w_(nx,ny,nz))
@ -223,6 +236,10 @@ module m_fdm_calc
v_=0.0
w_=0.0
allocate(dm(nx,ny,nz))
dm = diff
! DQ initializing
fdmsavecount=1 !FDM save count
sum_wrate=0.
@ -237,9 +254,27 @@ 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(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
@ -247,34 +282,6 @@ module m_fdm_calc
CALL ludcmp(nx,ny,nz_all,1,0,0)
CALL MPI_TYPE_VECTOR (nz, nx*nz, nx*ny, MPI_REAL8, vtype1, mpi_err)
CALL MPI_TYPE_COMMIT (vtype1, mpi_err)
CALL MPI_TYPE_VECTOR (nz, 2*nx*nz, 2*nx*ny, MPI_REAL8, vtype2, mpi_err)
CALL MPI_TYPE_COMMIT (vtype2, mpi_err)
do i = 1, (ny/nz)
scnt1(i) = 1
rcnt1(i) = nx * nz * nz
sdisp1(i) = (i-1) * nx * nz * 8
rdisp1(i) = (i-1) * nx * nz * nz * 8
stype1(i) = vtype1
rtype1(i) = MPI_REAL8
scnt2(i) = 1
rcnt2(i) = 2 * nx * nz * nz
sdisp2(i) = (i-1) * 2 * nx * nz * 8
rdisp2(i) = (i-1) * 2 * nx * nz * nz * 8
stype2(i) = vtype2
rtype2(i) = MPI_REAL8
end do
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! Kwon
@ -295,11 +302,11 @@ module m_fdm_calc
in_yr=yy(1) ! inlet_Yr
out_yr=yy(nx) ! outlet_Yr
do i=1,ny
do j=1,nz
do j=1,nz
do i=1,ny
do ii=1,nx
y1(1,ii,i,j)=1. ! rho initializing
y1(2,ii,i,j)=yy(ii) ! Yr initializing
y1(ii,i,j,1)=1. ! rho initializing
y1(ii,i,j,2)=yy(ii) ! Yr initializing
enddo
enddo
enddo
@ -310,7 +317,7 @@ module m_fdm_calc
DO ii=1,nz
DO j=1,ny
DO i=1,nx
yr=y1(2,i,j,ii)/y1(1,i,j,ii)
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))
@ -330,9 +337,10 @@ module m_fdm_calc
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
fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz_all)))
sumc1=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.))
CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
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.
@ -423,7 +431,7 @@ module m_fdm_calc
do k=1,nz
do j=1,ny
do i=1,nx
yr=y1(2,i,j,k)/y1(1,i,j,k)
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
@ -444,7 +452,8 @@ module m_fdm_calc
if (myid.eq.0) write(*,633) sum_wrate
633 format (' ** Consumption Speed, Sc = ',f7.4)
fl_location=(hx*REAL(nx-1))*(1.-(sumc/(REAL(nx*ny*nz_all))))
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.
@ -483,126 +492,190 @@ module m_fdm_calc
implicit none
integer :: i,j,k,xx,yy,zz,ii
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_(2,xx,yy,zz),f_(2,xx,yy,zz)
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 :: ux(4,xx),dux(4,xx),d2ux(xx)
real*8 :: uy(4,yy),duy(4,yy),d2uy(yy)
real*8 :: uz(4,yy),duz(4,yy),d2uz(yy)
real*8 :: Ly(yy,zz),Dy
real*8 :: r1_t(2,xx,zz,yy),f_t(2,xx,zz,yy)
real*8 :: ww_t(xx,zz,yy)
real*8 :: uz(4,zz),duz(4,zz),d2uz(zz)
real*8 :: uuz(zz)
real*8 :: y
real*8 :: T0, T1, Ts, As
!CALL MPI_ALLTOALLW (ww_, scnt1, sdisp1, stype1, ww_t, rcnt1, rdisp1, rtype1, MPI_COMM_TASK, mpi_err)
f_(:,:,:,1) = 0.0 ! continuity
!CALL MPI_ALLTOALLW (r1_, scnt2, sdisp2, stype2, r1_t, rcnt2, rdisp2, rtype2, MPI_COMM_TASK, mpi_err)
! diffusivity
T0 = 1.0
T1 = (1.0 + bc)
do i = 1,numprocs
CALL MPI_ISEND(ww_(1,zz*(i-1)+1,1), scnt1(i), stype1(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
end do
do i = 1,numprocs
CALL MPI_RECV(ww_t(1,1,zz*(i-1)+1), rcnt1(i), rtype1(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
end do
Ts = (dm_ratio*T1 - (T1**(3./2.))) / (T1**(3./2.) - dm_ratio)
As = (T0 + Ts)
do i = 1,numprocs
CALL MPI_ISEND(r1_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
end do
do i = 1,numprocs
CALL MPI_RECV(r1_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
end do
DO k=1,zz
DO j=1,yy
DO i=1,xx
!! z-direction
DO j=1,zz
DO i=1,xx
DO k=1,yy
uz(1,k)=r1_t(1,i,j,k) ! 1:rho
uz(2,k)=r1_t(2,i,j,k)/r1_t(1,i,j,k) ! 2:Y
uz(3,k)=uz(1,k)*ww_t(i,j,k) ! 3:rho*w
uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y
ENDDO
y=(1.0 - r1_(i,j,k,2)/r1_(i,j,k,1)) * bc + 1.0
CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3)
CALL d2fp(yy,hy,uz(2,:),d2uz(:),1,3)
dm(i,j,k) = diff * As * sqrt(y) / (1. + Ts/y)
DO k=1,yy
! -( d(rho*w)/dz )
f_t(1,i,j,k)=-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_t(2,i,j,k)=-duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv.
ENDDO
ENDDO
ENDDO
!CALL MPI_ALLTOALLW (f_t, rcnt2, rdisp2, rtype2, f_, scnt2, sdisp2, stype2, MPI_COMM_TASK, mpi_err)
if (myid.eq.0) write(*,*) 'min(dm)',minval(dm),'max(dm)',maxval(dm)
do i = 1,numprocs
CALL MPI_ISEND(f_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
end do
do i = 1,numprocs
CALL MPI_RECV(f_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
end do
! x-direction
! reaction source term
DO k=1,zz
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(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
ENDDO
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
CALL d2fnonp(xx,hx,ux(2,:),d2ux(:),1,1)
y=r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y
DO i=1,xx
wrate=pre*ux(2,i)*exp(-ac/(1.+bc*(1.-ux(2,i)))) !wrate
IF ((1.-ux(2,i)).le.c_ref) THEN
wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate
IF ((1.-y).le.c_ref) THEN
wrate=min_wr
IF ((1.-ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(2,i)-c_ref))+ &
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 )
f_(1,i,j,k) = f_(1,i,j,k) - dux(3,i) ! continuity
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 )
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)) - wrate ! species conservation
! + (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
ENDDO
ENDDO
!! y-direction
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
ENDDO
DO k=1,zz
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
CALL d2fp(yy,hy,uy(2,:),d2uy(:),1,2)
! -( d(rho*v)/dy )
DO j=1,yy
! -( d(rho*v)/dy )
f_(1,i,j,k)=f_(1,i,j,k)-duy(3,j) ! continuity
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 )
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.
ENDDO
ENDDO
! + (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
@ -610,29 +683,103 @@ module m_fdm_calc
DO k=1,zz
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_(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_(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_(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(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
real*8 :: yy1(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=2; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=3; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=4; CALL substep(yy1,yy2,yy1,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
istage=5; CALL substep(yy2,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
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
@ -643,21 +790,117 @@ 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(xx,yy,zz,neq),yy2(xx,yy,zz,neq),yyf(xx,yy,zz,neq)
CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
! 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(2,xx,yy,zz),r1(2,xx,yy,zz),r2(2,xx,yy,zz),f(2,xx,yy,zz)
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
@ -668,6 +911,7 @@ module m_fdm_calc
a(3)= 2251764453980./15575788980749.
a(4)= 26877169314380./34165994151039.
a(5)=0.
b(1)= 1153189308089./22510343858157.
b(2)= 1772645290293./4653164025191.
b(3)= -1672844663538./4480602732383.
@ -677,43 +921,146 @@ module m_fdm_calc
CALL fns(ri,f,xx,yy,zz,uu_,vv_,ww_)
IF(istage<5) THEN
at=a(istage)*fdmdt
bt=(b(istage)-a(istage))*fdmdt
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
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
ENDDO
ENDDO
ENDDO
ELSE
bt=b(istage)*fdmdt
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
r1(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
ENDDO
r2 = r1 + bt*f
r1 = r1 + at*f
ELSE
bt=b(istage)*fdmdt
r1 = r1 + bt*f
!==========rho=1 treatment
r1(2,i,j,k)=r1(2,i,j,k)/r1(1,i,j,k)
r1(1,i,j,k)=1.
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(2,i,j,k)=MIN(in_yr,r1(2,i,j,k))
r1(i,j,k,2)=MIN(in_yr,r1(i,j,k,2))
!==========Min Yr=0 treatment
! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))
! 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

2
m_io.f90
View file

@ -9,7 +9,7 @@ module m_io
use m_openmpi
implicit none
character*6 :: file_ext, ex_file_ext
character*8 :: file_ext, ex_file_ext
character*80 :: fname
! output file handle
integer :: in=10, out=11

3
m_mpi.f90
View file

@ -64,7 +64,8 @@ contains
! initializing MPI environment
call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
!call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
call MPI_INIT(mpi_err)
call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err)
call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err)

3
m_openmpi.f90
View file

@ -64,7 +64,8 @@ contains
! initializing MPI environment
call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
!call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
call MPI_INIT(mpi_err)
call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err)
call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err)

4
m_parameters.f90
View file

@ -114,6 +114,7 @@ module m_parameters
real*8 :: ac,bc,c_cut,c_ref,min_wr,prof_wr,hx,hy,diff,pre
real*8 :: dt_fullsave,dt_fdmsave,t_fullsave,t_fdmsave
real*8 :: fdmcfl,startad,dummyu_,schmidt
real*8 :: dm_ratio
real*8 :: bgn_save
! files for statistics
@ -581,6 +582,9 @@ contains
read(in,*,ERR=9000,END=9000) fdmcfl
write(out,'(a30,3x,e20.10)') 'FDM_CFL =',fdmcfl
read(in,*,ERR=9000,END=9000) dm_ratio
write(out,'(a30,3x,e20.10)') 'Dm_Ratio =',dm_ratio
! initialize ----------------------
t_fdmsave=dt_fdmsave
t_fullsave=dt_fullsave

4
m_timing.f90
View file

@ -19,8 +19,8 @@ contains
! make the default job runlimit to be 6 months
job_runlimit = 6 * 30 * 24 * 60
! make the default job runlimit to be 12 hours
job_runlimit = 12 * 60
! make the default job runlimit to be 12 hours -> one month
job_runlimit = 12 * 60 * 30 * 2
write(out,'(a30,3x,i20)') 'job_runlimit (default):',job_runlimit

5
main.f90
View file

@ -129,9 +129,10 @@ program x_code
end if
end if
if (fdm_sw.ne.0) then ! Reacting
! reacting scalars
if (fdm_sw.ne.0) then
call gather_4
call fdm_exe
call fdm_exe
endif
! RHS for scalars

4
write_stats.f90
View file

@ -31,6 +31,7 @@ subroutine write_tp_field
do i=1,nx,spx
write(tp_field,'(6e18.9)') dx*real(i),dy*real(j),dz*real(k),u_(i,j,k),v_(i,j,k),w_(i,j,k) ! u',v',w'
enddo
call flush(tp_field)
enddo
enddo
@ -49,13 +50,12 @@ subroutine write_tp_field
(1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
v_(i,j,ii),w_(i,j,ii)
enddo
call flush(tp_field)
enddo
enddo
endif
call flush(tp_field)
897 format('ZONE T="time= ',f10.5,'" I= ',i4,' J= ',i4,' K= ',i4,' F=POINT')
end subroutine write_tp_field

2
x_fftw.f90
View file

@ -298,7 +298,7 @@ CONTAINS
! Definition of the array ialias.
! The array ialias is just the number of wavenumbers at (i,j,k) that have
! their magnitude higher than nx/3. This is needed in dealiasing procedures.
rnx3 = real(nx/3, 8)
rnx3 = real(2*nx, 8) / real(3*lx, 8)
do k = 1,nz
if (abs(aky(k)) .gt. rnx3) ialias(:,:,k) = 1
do j = 1,ny