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

Original HIT Source

Browse source
This commit is contained in:
ignis 2014-04-16 15:17:00 +09:00
commit 4de8a61a0e
34 changed files with 11618 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
.gitignore vendored Normal file
View file

@ -0,0 +1,5 @@
*.[oa]
*.mod
*.x
scripts/*
tags

123
Makefile Normal file
View file

@ -0,0 +1,123 @@
# Makefile with logic in it. For different hosts runs different things
# Define the correct flags and compilers
MPIF90 = blah
ifeq ($(HOSTNAME), ace2-1.postech.ac.kr)
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
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
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
endif
# Program name
PROG = hit3d.x
# Modules
MODULES = m_openmpi.o\
m_io.o\
m_parameters.o\
m_work.o\
m_fields.o\
m_timing.o\
x_fftw.o\
m_filter_xfftw.o\
m_particles.o\
m_stats.o\
m_force.o\
m_rand_knuth.o\
m_les.o\
RANDu.o
# Objects
OBJ = main.o\
begin_new.o\
begin_restart.o\
dealias_all.o\
get_file_ext.o\
init_velocity.o\
init_scalars.o\
io_write_4.o\
my_dt.o\
my_exit.o\
pressure.o\
restart_io.o\
rhs_velocity.o\
rhs_scalars.o\
velocity_rescale.o\
write_tmp4.o
# -------------------------------------------------------
# link
$(PROG): $(MODULES) $(OBJ)
$(MPIF90) $(MODULES) $(OBJ) -o $(PROG) $(LDFLAGS)
# -------------------------------------------------------
# compile
$(OBJ): $(MODULES)
%.o: %.f
$(MPIF90) $(FCFLAGS) $(FCFLAGS_F77) $<
%.o: %.f90
$(MPIF90) $(FCFLAGS) $<
clean:
rm *.o *.mod $(PROG)

115
Makefile_Orig Normal file
View file

@ -0,0 +1,115 @@
# Makefile with logic in it. For different hosts runs different things
# Define the correct flags and compilers
MPIF90 = blah
# 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
# WCR cluster, Stanford
ifeq ($(HOSTNAME), glacial.stanford.edu)
FFTW_HOME = /home/chumakov/packages/fftw-3.2/pgi
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
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
endif
# Program name
PROG = hit3d.x
# Modules
MODULES = m_openmpi.o\
m_io.o\
m_parameters.o\
m_work.o\
m_fields.o\
m_timing.o\
x_fftw.o\
m_filter_xfftw.o\
m_particles.o\
m_stats.o\
m_force.o\
m_rand_knuth.o\
m_les.o\
RANDu.o
# Objects
OBJ = main.o\
begin_new.o\
begin_restart.o\
dealias_all.o\
get_file_ext.o\
init_velocity.o\
init_scalars.o\
io_write_4.o\
my_dt.o\
my_exit.o\
pressure.o\
restart_io.o\
rhs_velocity.o\
rhs_scalars.o\
velocity_rescale.o\
write_tmp4.o
# -------------------------------------------------------
# link
$(PROG): $(MODULES) $(OBJ)
$(MPIF90) $(MODULES) $(OBJ) -o $(PROG) $(LDFLAGS)
# -------------------------------------------------------
# compile
$(OBJ): $(MODULES)
%.o: %.f
$(MPIF90) $(FCFLAGS) $(FCFLAGS_F77) $<
%.o: %.f90
$(MPIF90) $(FCFLAGS) $<
clean:
rm *.o *.mod $(PROG)

54
RANDu.f Normal file
View file

@ -0,0 +1,54 @@
module RANDu
c..
implicit none
real*8 :: rseed
c..
c double precision random
c..
contains
c----------------------------------------------------
function random(idum3)
c..
c.. initialize with idum<0; Then keep idum>0 unchanged for the same
c.. sequence
c..
integer idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
double precision random,am,eps,rnmx,idum3
parameter (im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1,
+ ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,
+ ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.0-eps)
integer idum2,j,k,iv(ntab),iy
c..
save iv,iy,idum2
data idum2/123456789/, iv/ntab*0/, iy/0/
c..
idum=nint(idum3/987)
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do 11 j=ntab+8,1,-1
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if(idum.lt.0) idum=idum+im1
if(j.le.ntab) iv(j)=idum
11 continue
endif
c.. start here when not initializing
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
if (idum.lt.0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
if (idum2.lt.0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1) iy=iy+imm1
random=min(am*iy,rnmx)
c..
return
end function random
c--------------------------------------------------
end module RANDu

165
README Normal file
View file

@ -0,0 +1,165 @@
--------------------------------------------------------------------------------
GENERAL
--------------------------------------------------------------------------------
HIT3DP is a pseudospectral DNS code, that is, it performs direct numerical
simulation of incompressible isotripic homogeneous turbulence with or without
forcing. The code has capability of carrying passive scalars, Lagrangian
particles and Large Eddy Simulation
--------------------------------------------------------------------------------
EXTRA PACKAGES
--------------------------------------------------------------------------------
The code is written in Fortran 90 and uses the two open libraries:
- Open-MPI (www.open-mpi.org)
- FFTW3 (www.fftw.org)
--------------------------------------------------------------------------------
LICENSING
--------------------------------------------------------------------------------
The code is distributed under the terms of GNU GPLv3 license. You can read
the full text of the license at http://www.gnu.org/licenses/gpl.html
Copyright (C) 2006-2010 Sergei Chumakov, Natalia Vladimirova, Misha Stepanov
--------------------------------------------------------------------------------
COMPILING THE CODE
--------------------------------------------------------------------------------
First, edit the Makefile:
- add a section that corresponds to the name of your machine. Ideally it should
be a wrapper from your MPI implementation.
- define the name of the F90 compiler
- define FCFLAGS and LDFLAGS. They should include the include directories, the
flags that link FFTW3 and MPI implementation.
Run "gmake".
--------------------------------------------------------------------------------
RUNNING THE CODE
--------------------------------------------------------------------------------
The directory "scripts" provides some examples of the batch job submission files.
The directory "scripts" contains the following files:
00_example.in a sample input file
snapshot.gp a Gnuplot instruction file that creates two plots that
can get attached to the notification emails
coyote.sub Running script for the Coyote cluster at LANL
wcr.sub Example script for WCR cluster at Center for Turbuience Research
at Stanford University
--------------------------------------------------------------------------------
THE INPUT FILE
--------------------------------------------------------------------------------
NX,NY,NZ Number of grid points in one dimension. The grid will be NX x NY x NZ.
The physical dimensions will be 2*pi x 2*pi x 2*pi
ITMIN The timestep number of the restart file. The restart files have names
such as "test__this.64.123456". Here, "test__this" is the run name,
"64" signifies that the file is written with double precision and
"123456" is the timestep number. If the ITMIN is set to 0, the
subroutine that defines the initial conditionis for the flow is called.
ITMAX The maximum number of timesteps in the simulation.
IPRNT1 How often to generate the statistics.
IPRNT2 How often to write restart files
IWRITE4 How often to write the real*4 files that are used for post-processing.
TMAX The runtime of the simulation (not the wallclocok time)
TRESCALE The time at which to rescale the velocity. This is used in decaying
simulations when we want to establish some correlations first and
then rescale the velocity field so it has higher kinetic energy.
TSCALAR When to start moving the passive scalars.
flow_type Parameter that switches the flow type
0 - decaying turbulence
1 - forced turbulence
RE The local Reynolds number (1/nu, where nu is viscosity)
DT The timestep.
If DT is negative, then the timestep is fixed to be (-DT)
If DT is positive, the timestep is found from the stability
criteria for the time-stepping scheme that is used.
ISPCV1 Initial spectrum type (see init_velocity.f90)
mv1 initial infrared exponent in the spectrum
wm0v1 initial peak wavenumber in the spectrum
force_type The type of the forcing that is applied for the case of
forced turbulence.
1 - forcing from Michaels PRL paper (PRL #79(18) p.3411)
So far no other forcing has been implemented
KFMAX The upper bound for the forcing band in the Fourier space.
FAMP The magnitude of the forcing (usually set to 0.5)
det_rand The parameter that switches the random generation for the
random seeds for the code.
DEFUNCTIONAL. In the current version of the code, the seeds for the
random number generator are fixed and are taken from the input file.
The fixed seeds have the effect of producing the initial data that
looks similar for different resolutions (the large features of
initial flow in 32^3 simulation will look similar to the large features
of a 1024^3 simulation if the seeds are the same).
RN1, RN2, RN3 - random number seeds
DEALIAS The parameter that switches between the dealiasing algorithms.
0 - the standard 3/2-rule (or 2/3 rule). Faster computations, but
fewer modes.
1 - the phase shift combined with truncation. This retains much more
modes than the 2/3-rule, while increasing the computations 60% or so.
The most economical mode for DNS in terms of flops per the number of
Fourier modes in the resulting data.
np The number of Lagrangian particles
* particle tracking mechanism:
0 - trilinear interpolations
1 - 4-point cubic interpolation
time_p time in the simulation when to release the particles in the flow
particle_filter_size
The particles can be advected by fully resolved field or by locally averaged
field. The filter size determines the size of the filter that is applied
to the velocity field before computing the particles' velocities.
les_model The LES model. See m_les.f90 for list of the current models.
NUMS The number of passive scalars to carry around
The last section contains the parameters of the passive scalars. Each scalar
must have the type, Schmidt number, infrared exponent, peak wavenumber and
reaction rate.
TYPE:
0 The scalar that is forced by the mean gradient.
1-9 The initial conditions for the scalar are generated using Fourier space.
1: Exponential spectrum
2: von Karman spectrum
3: double-delta PDF
>10 The initial conditions for the scalar are generated in the real space.
11: single slab of the scalar.
12: two slabs of the scalar
13: 1-D sinusoidal wave of the scalar
The reaction rate parameter is defunctional in this version of the code.
--------------------------------------------------------------------------------
ANY QUESTIONS?
--------------------------------------------------------------------------------
email Sergei Chumakov at chumakov@stanford.edu

34
begin_new.f90 Normal file
View file

@ -0,0 +1,34 @@
subroutine begin_new
use m_openmpi
use m_parameters
use m_fields
use m_stats
implicit none
! defining time
TIME = zip
! deciding if we advance scalars or not
if (TSCALAR.le.zip .and. n_scalars.gt.0) int_scalars = .true.
! defining the iteration number
ITIME = 0
file_ext = '000000'
if (task.eq.'hydro') then
call init_velocity
if (n_scalars.gt.0) call init_scalars
call io_write_4
end if
if (task_split) then
call fields_to_stats
if (task.eq.'stats') call stat_main
end if
return
end subroutine begin_new

88
begin_restart.f90 Normal file
View file

@ -0,0 +1,88 @@
subroutine begin_restart
use m_parameters
use m_particles
use m_fields
use m_timing
implicit none
real*8, allocatable :: zhopa(:,:,:,:)
integer :: i
write(out,"(70('='))")
write(out,*) ' RESTART '
write(out,"(70('='))")
call flush(out)
ITIME = ITMIN
call get_file_ext
! reading the restart file
if (task.eq.'hydro') call restart_read_parallel
! broadcasting the current simulation time from the restart file
call MPI_BCAST(time,1,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
! redefining the timestep.
! the code uses Adams-Bashforth time-stepping scheme,
! which is 2nd order accurate in time. After the restart, thefirst
! timestep is done using a simple Euler scheme (1st order in time).
! To help maintain any sensible accuracy, we need to start with
! the timestep which is smaller than the last.
if (variable_dt) then
dt = half * dt
call MPI_BCAST(dt,1,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
write(out,*) "Making the timestep smaller: ",dt
call flush(out)
end if
!!$!================================================================================
!!$! Checking the parallel read speed (reading 100 times)
!!$! Currently the speedup factor from using parallel read is about 2.5
!!$!================================================================================
!!$ if(task.eq.'hydro') then
!!$ call m_timing_check
!!$ write(out,*) 'Start! ',cpu_min,cpu_sec
!!$ do i = 1,200
!!$ call restart_read! _parallel
!!$ end do
!!$ call m_timing_check
!!$ write(out,*) 'Finish!',cpu_min,cpu_sec
!!$
!!$ end if
!!$ call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
!!$ stop 'checking the restart read speed'
!!$
!!$!================================================================================
!!$! Checking the accuracy of the parallel read
!!$! (reading parallel first, then serial and comparing)
!!$!================================================================================
! if (task.eq.'hydro') then
!
! call restart_read_parallel
!
! allocate(zhopa(nx+2,ny,nz,3+n_scalars), stat=ierr)
! if (ierr.ne.0) stop 'cannot allocate zhopa'
! zhopa = zip
! zhopa = fields
!
! call restart_read
!
! print "(10e15.6)",maxval(abs(zhopa-fields))
! end if
! call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
! stop 'checking the restart read'
!!$!================================================================================
! deciding whether we advance scalars or not
if (n_scalars.gt.0 .and. time.gt.TSCALAR) then
int_scalars = .true.
write(out,"('Advancing ',i3,' scalars.')") n_scalars
end if
if (task.eq.'parts') call particles_init
return
end subroutine begin_restart

35
dealias_all.f90 Normal file
View file

@ -0,0 +1,35 @@
subroutine dealias_all
use m_parameters
use m_fields
use x_fftw
implicit none
integer :: i,j,k
real*8 :: akmax
! "2/3-rule". For good explanation, see the following paper:
! R.S.Rogallo, "Numerical Exoeriments in Homogeneous Turbuelnce"
! NASA Technical Memorandum, NASA 1981.
! Or contact authors of the code, they should have a PDF copy somewhere.
!
! truncating all the modes in which at least one component of the k-vector
! has magnitude that is larger than N/3
akmax = real(kmax,8)
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
if (abs(akx(i)) .gt. akmax .or. &
abs(aky(k)) .gt. akmax .or. &
abs(akz(j)) .gt. akmax) then
fields(i,j,k,1:3+n_scalars) = zip
end if
end do
end do
end do
return
end subroutine dealias_all

11
get_file_ext.f90 Normal file
View file

@ -0,0 +1,11 @@
!================================================================================
!================================================================================
subroutine get_file_ext
use m_parameters, only : ITIME
use m_io, only : file_ext
implicit none
write(file_ext,"(i6.6)") itime
return
end subroutine get_file_ext

383
init_scalars.f90 Normal file
View file

@ -0,0 +1,383 @@
!================================================================================
!================================================================================
! Initialization of passive scalars
!================================================================================
subroutine init_scalars
use m_parameters
use m_io
use m_fields
use x_fftw, only : ialias
implicit none
integer :: n_scalar, i, j, k
write(out,*) 'Generating scalars'
call flush(out)
do n_scalar = 1,n_scalars
call init_scalar(n_scalar)
end do
write(out,*) "Generated the scalars."
call flush(out)
! now making sure that the scalars do not have any high
! Fourier harmonics by zeroing out everything that has a wavenumber
! that potentially can produce aliasing
do k = 1, nz
do j = 1, ny
do i = 1, nx+2
if (ialias(i,j,k).gt.0) fields(i,j,k,4:3+n_scalars) = zip
end do
end do
end do
return
end subroutine init_scalars
!================================================================================
!================================================================================
subroutine init_scalar(n_scalar)
use m_parameters
use m_io
use m_fields
implicit none
integer :: n_scalar, ic_type, sc_type
write(out,*) 'Generating scalar #', n_scalar
call flush(out)
sc_type = scalar_type(n_scalar)
ic_type = sc_type - (sc_type/100)*100
if (ic_type.eq.0) then
! gradient source - no need for initial conditions
! thus making the initial scalar field zero
write(out,*) "The scalar type = 0, nothing to generate"
call flush(out)
fields(:,:,:,n_scalar+3) = zip
elseif (ic_type.lt.10) then
! if the last two digits of the scalar type are less than 10,
! the scalar initial conditions are generated based on the
! particular spectrum of the scalar
call init_scalar_spectrum(n_scalar)
else
! if the last two digits are bigger than 10, the scalar is generated
! in physical space and then transformed in the Fourier space
call init_scalar_space(n_scalar)
end if
return
end subroutine init_scalar
!================================================================================
!================================================================================
subroutine init_scalar_spectrum(n_scalar)
!================================================================================
use m_openmpi
use m_parameters
use m_io
use m_fields
use m_work
use x_fftw
use m_rand_knuth
use RANDu
implicit none
integer :: i, j, k, n, n_scalar
integer *8 :: i8
real*8, allocatable :: e_spec(:), e_spec1(:), rr(:)
integer *8, allocatable :: hits(:), hits1(:)
integer :: n_shell
real*8 :: sc_rad1, sc_rad2
real*8 :: wmag, wmag2, ratio, fac, fac2
!--------------------------------------------------------------------------------
write(out,*) " Generating scalar # ",n_scalar
call flush(out)
! Initializing the random sequence with the seed RN2
fac = random(-RN2)
! allocate work arrays
allocate( e_spec(kmax), e_spec1(kmax), hits(kmax), hits1(kmax), &
rr(nx+2), stat=ierr)
! bringing the processors to their own places in the random sequence
! ("2" is there because we're generating two random number fields
! for each scalar field
! using i8 because it's int*8
do i8 = 1,myid*(nx+2)*ny*nz*2
fac = random(RN2)
end do
! now filling the arrays wrk1, wrk2
do n = 1,2
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
wrk(i,j,k,n) = random(RN2)
end do
end do
end do
end do
! bringing the random numbers to the same
! point in the sequence again
do i8 = 1,(numprocs-myid-1)*(nx+2)*ny*nz*2
fac = random(RN2)
end do
! making random array with Gaussian PDF
! out of the two arrays that we generated
wrk(:,:,:,3) = sqrt(-two*log(wrk(:,:,:,1))) * sin(TWO_PI*wrk(:,:,:,2))
! go to Fourier space
call xFFT3d(1,3)
!-------------------------------------------------------------------------------
! Calculating the scalar spectrum
!-------------------------------------------------------------------------------
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
e_spec1 = zip
e_spec = zip
hits = 0
hits1 = 0
! assembling the scalar energy in each shell and number of hits in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
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) then
fac2 = fac * wrk(i,j,k,3)**2
if (akx(i).eq.0.d0) fac2 = fac2 * 0.5d0
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! broadcasting the spectrum
count = kmax
call MPI_BCAST(e_spec,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
!-------------------------------------------------------------------------------
! Now make the spectrum to be as desired
!-------------------------------------------------------------------------------
! first, define the desired spectrum
do k = 1,kmax
wmag = real(k, 8)
ratio = wmag / peak_wavenum_sc(n_scalar)
if (scalar_type(n_scalar).eq.0) then
! Plain Kolmogorov spectrum
e_spec1(k) = wmag**(-5.d0/3.d0)
else if (scalar_type(n_scalar).eq.1 .or. scalar_type(n_scalar).eq.3) then
! Exponential spectrum
e_spec1(k) = ratio**3 / peak_wavenum_sc(n_scalar) * exp(-3.0D0*ratio)
else if (scalar_type(n_scalar).eq.2) then
! Von Karman spectrum
fac = two * PI * ratio
e_spec1(k) = fac**4 / (one + fac**2)**3
else
write(out,*) "INIT_SCALARS: WRONG INITIAL SPECTRUM TYPE: ",scalar_type(n_scalar)
call flush(out)
stop
end if
end do
! normalize it so it has the unit total energy
e_spec1 = e_spec1 / sum(e_spec1(1:kmax))
! now go over all Fourier shells and multiply the velocities in a shell by
! the sqrt of ratio of the resired to the current spectrum
fields(:,:,:,3+n_scalar) = zip
do k = 1,nz
do j = 1,ny
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
fields(i,j,k,3+n_scalar) = wrk(i,j,k,3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
else
fields(i,j,k,3+n_scalar) = zip
end if
end do
end do
end do
!-------------------------------------------------------------------------------
! Creating scalars with double-delta PDF
!-------------------------------------------------------------------------------
if (scalar_type(n_scalar).eq.3) then
wrk(:,:,:,0) = fields(:,:,:,3+n_scalar)
call xFFT3d(-1,0)
! making it double-delta (0.9 and -0.9)
wrk(:,:,:,0) = sign(one,wrk(:,:,:,0)) * 0.9d0
call xFFT3d(1,0)
! smoothing it by zeroing out high harmonics
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .eq. 0 .or. n_shell .ge. kmax*2/3) then
wrk(i,j,k,0) = zip
end if
end do
end do
end do
fields(:,:,:,3+n_scalar) = wrk(:,:,:,0)
end if
! deallocate work arrays
deallocate(e_spec, e_spec1, rr, hits, hits1, stat=ierr)
return
end subroutine init_scalar_spectrum
!================================================================================
!================================================================================
!================================================================================
!================================================================================
subroutine init_scalar_space(n_scalar)
!================================================================================
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: i, k, n_scalar, sc_type, ic_type, nfi
real*8 :: zloc, sctmp, h, xx
nfi = 3 + n_scalar
write(out,*) " Generating scalar # ", n_scalar
call flush(out)
sc_type = scalar_type(n_scalar)
ic_type = sc_type - (sc_type/100)*100
select case (ic_type)
!---------------------------------------------------
! single slab of the scalar
!---------------------------------------------------
case(11)
! how much to smear out the interface
! the minimum interface length is set to 2*dz
! the maximum is set to 2*PI/(peak wavenumber for the scalar, taken from
! the .in file)
h = max(2.*dz, two*PI/peak_wavenum_sc(n_scalar))
! creating array of scalar
do k = 1,nz
zloc = dble(myid*nz + k-1) * dz
sctmp = tanh((zloc-PI*0.5)/h) - tanh((zloc-PI*1.5)/h) - one
wrk(:,:,k,0) = sctmp
end do
! FFT of the scalar
call xFFT3d(1,0)
! putting it into the scalar array
fields(:,:,:, 3+n_scalar) = wrk(:,:,:,0)
if (iammaster) fields(1,1,1,3+n_scalar) = zip
!---------------------------------------------------
! two slabs of the scalar
!---------------------------------------------------
case(12)
write(out,*) "-- Double-slab scalar in real space"
call flush(out)
! how much to smear out the interface
h = max(8.*dz, PI/8.d0)
! creating array of scalar
do i = 1,nx
xx = dble(i-1) * dx
sctmp = tanh((xx-PI*0.25)/h) - tanh((xx-PI*0.75)/h) + tanh((xx-PI*1.25)/h) - tanh((xx-PI*1.75)/h)
wrk(i,:,:,0) = sctmp - one
! now it is between -1 and 1
end do
! FFT of the scalar
call xFFT3d(1,0)
! putting it into the scalar array
fields(:,:,:, 3+n_scalar) = wrk(:,:,:,0)
! making sure that the mean is ero
if (iammaster) fields(1,1,1,3+n_scalar) = zip
!---------------------------------------------------
! N slabs of the scalar
! actually more like N sinusoidal waves
!---------------------------------------------------
case(13)
write(out,*) "-- Multi-slab scalar in real space"
call flush(out)
! making sure that we can support the desired numberof waves with our FFT
peak_wavenum_sc(n_scalar) = min(peak_wavenum_sc(n_scalar),real(nx/8))
! definition of initial scalar
fields(:,:,:,nfi) = zip
do i = 1, nx
fields(i,:,:,nfi) = sin(peak_wavenum_sc(n_scalar) * dx * real(i-1))
end do
call xFFT3d_fields(1,nfi)
case default
write(out,*) "INIT_SCALARS: UNEXPECTED SCALAR TYPE: ", scalar_type(n_scalar)
call flush(out)
stop
end select
write(out,*) "Initialized the scalars."
call flush(out)
return
end subroutine init_scalar_space

295
init_velocity.f90 Normal file
View file

@ -0,0 +1,295 @@
subroutine init_velocity
use m_openmpi
use m_parameters
use m_io
use m_fields
use m_work
use x_fftw
use m_rand_knuth
use RANDu
implicit none
integer :: i, j, k, n
integer*8 :: seed1, seed2, i8, j8, k8
integer :: time_array(8)
real, allocatable :: rr(:)
real*8, allocatable :: e_spec(:), e_spec1(:)
integer *8, allocatable :: hits(:), hits1(:)
integer :: n_shell
real*8 :: sc_rad1, sc_rad2
real*8 :: wmag, wmag2, ratio, fac, fac2
!--------------------------------------------------------------------------------
! First, if it's a Taylor-Green vortex, then initialize and quit
!--------------------------------------------------------------------------------
if (isp_type .eq. -1) then
call init_velocity_taylor_green
return
end if
!================================================================================
allocate( e_spec(kmax), e_spec1(kmax), rr(nx+2), hits(kmax), hits1(kmax), stat=ierr)
if (ierr.ne.0) stop "cannot allocate the init_velocity arrays"
write(out,*) 'generating random velocities'
call flush(out)
!-------------------------------------------------------------------------------
! Generate the velocities
!-------------------------------------------------------------------------------
! initialize the random number sequence by the seed from the first processor
if (myid.eq.0) call system_clock(seed1,seed2)
if (myid.eq.0) then
call date_and_time(values=time_array)
seed1 = time_array(8)
end if
count = 1
call MPI_BCAST(seed1,count,MPI_INTEGER8,0,MPI_COMM_TASK,mpi_err)
!!$ seed1 = 23498675
!!$ call rand_knuth_init(seed1)
!!$ write(out,*) 'seed1 = ',seed1
!!$ call flush(out)
seed1 = RN1
write(out,*) "RANDOM SEED FOR VELOCITIES = ", seed1
call flush(out)
rseed = real(seed1,8)
fac = random(-rseed)
! bringing the processors to their own places in the random sequence
! ("6" is there because we're generating six fields
! using seed1 because it's int*8
!!$ write(out,*) "Will scroll down to my initial position", myid*ny*nz*6
!!$ call flush(out)
!!$ do i = 1,myid*ny*nz*6
!!$ call rand_knuth(rr,nx+2)
!!$ end do
do i8 = 1,myid*(nx+2)*ny*nz*6
fac = random(rseed)
end do
!!$ write(out,*) "Scrolled."
!!$ call flush(out)
! now filling the arrays wrk1...wrk6
do n = 1,6
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
wrk(i,j,k,n) = random(rseed)
end do
end do
end do
end do
! just in case, bringing the random numbers to the same
! point in the sequence again
do seed1 = 1,int((numprocs-myid-1)*(nx+2)*ny*nz*6,8)
fac = random(rseed)
end do
! making three random arrays with Gaussian PDF
! out of the six arrays that we generated
wrk(:,:,:,1:3) = sqrt(-two*log(wrk(:,:,:,1:3))) * sin(TWO_PI*wrk(:,:,:,4:6))
! --- Making three arrays that have Gaussian PDF and the incompressibility property
! go to Fourier space
do n = 1,3
call xFFT3d(1,n)
end do
! assemble the arrays in wrk4..6, only the wavenumbers below kmax
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
n_shell = nint(sqrt(akx(i)**2 + aky(k)**2 + akz(j)**2))
if (n_shell .gt. 0 .and. n_shell .le. kmax) then
wrk(i ,j,k,4) = - (aky(k)*wrk(i+1,j,k,3) - akz(j)*wrk(i+1,j,k,2))
wrk(i+1,j,k,4) = aky(k)*wrk(i ,j,k,3) - akz(j)*wrk(i ,j,k,2)
wrk(i ,j,k,5) = - (akz(j)*wrk(i+1,j,k,1) - akx(i+1)*wrk(i+1,j,k,3))
wrk(i+1,j,k,5) = akz(j)*wrk(i ,j,k,1) - akx(i )*wrk(i ,j,k,3)
wrk(i ,j,k,6) = - (akx(i+1)*wrk(i+1,j,k,2) - aky(k)*wrk(i+1,j,k,1))
wrk(i+1,j,k,6) = akx(i )*wrk(i ,j,k,2) - aky(k)*wrk(i ,j,k,1)
else
wrk(i:i+1,j,k,4:6) = zip
end if
end do
end do
end do
fields(:,:,:,1:3) = wrk(:,:,:,4:6)
!-------------------------------------------------------------------------------
! Making the spectrum to be what it should
!-------------------------------------------------------------------------------
! --- first get the energy spectrum (copied from m_stat.f90)
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
e_spec1 = zip
e_spec = zip
hits = 0
hits1 = 0
! assembling the total energy in each shell and number of hits in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
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) then
fac2 = fac * (fields(i,j,k,1)**2 + fields(i,j,k,2)**2 + fields(i,j,k,3)**2)
if (akx(i).eq.0.d0) fac2 = fac2 * 0.5d0
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
count = kmax
call MPI_BCAST(e_spec,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
!-------------------------------------------------------------------------------
! Now make the spectrum to be as desired
!-------------------------------------------------------------------------------
! first, define the desired spectrum
do k = 1,kmax
wmag = real(k, 8)
ratio = wmag / peak_wavenum
if (isp_type.eq.0) then
! Plain Kolmogorov spectrum
e_spec1(k) = wmag**(-5.d0/3.d0)
else if (isp_type.eq.1) then
! Exponential spectrum
e_spec1(k) = ratio**3 / peak_wavenum * exp(-3.0D0*ratio)
else if (isp_type.eq.3) then
! Von Karman spectrum
fac = two * PI * ratio
e_spec1(k) = fac**4 / (one + fac**2)**3
else
write(out,*) "ERROR: WRONG INITIAL SPECTRUM TYPE: ",isp_type
call flush(out)
stop
end if
end do
! normalize it so it has the unit total energy
e_spec1 = e_spec1 / sum(e_spec1(1:kmax))
! now go over all Fourier shells and multiply the velocities in a shell by
! the sqrt of ratio of the resired to the current spectrum
do k = 1,nz
do j = 1,ny
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
fields(i,j,k,1:3) = fields(i,j,k,1:3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
else
fields(i,j,k,1:3) = zip
end if
end do
end do
end do
write(out,*) "Generated the velocities."
call flush(out)
! deallocate work arrays
deallocate(e_spec, e_spec1, rr, hits, hits1, stat=ierr)
return
end subroutine init_velocity
!================================================================================
! Initialize the velocities with Taylor-Green vortex
!================================================================================
subroutine init_velocity_taylor_green
use m_openmpi
use m_parameters
use m_io
use m_fields
use m_work
use x_fftw
implicit none
logical :: verbose = .true.
integer :: i, j, k, n
real*8 :: xx, yy, zz
if (verbose) write(out,*) " --- Initial velocity field is Taylor-Green vortex"
if (verbose) call flush(out)
do k = 1, nz
zz = real(nz*myid + k - 1, 8) * dz
do j = 1, ny
yy = real(j-1, 8) * dy
do i = 1, nx
xx = real(i-1,8) *dx
wrk(i,j,k,1) = sin(xx) * cos(yy) * cos(zz)
wrk(i,j,k,2) = - cos(xx) * sin(yy) * cos(zz)
wrk(i,j,k,3) = zip
end do
end do
end do
call xFFT3D(1, 1)
call xFFT3D(1, 2)
call xFFT3D(1, 3)
fields(:,:,:,1:3) = wrk(:,:,:,1:3)
if (verbose) write(out,*) " --- initialized."
if (verbose) call flush(out)
return
end subroutine init_velocity_taylor_green

92
io_write_4.f90 Normal file
View file

@ -0,0 +1,92 @@
subroutine io_write_4
! Writing out the velocities and scalars in X-space
! to the real*4 file
use m_parameters
use m_fields
use m_work
use x_fftw
use m_les
implicit none
integer :: n_out, n, i, j, k
real*8 :: wmag2, rkmax2
! every variable will undergo a mode truncation for all modes
! that are higher than kmax. This will ensure that the written
! variables are isotropic
rkmax2 = real(kmax,8)**2
! 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
do j = 1,ny
do i = 1,nx+2
wmag2 = akx(i)**2 + aky(k)**2 + akz(j)**2
if (wmag2 .gt. rkmax2) then
wrk(i,j,k,1:n_out) = zip
else
wrk(i,j,k,1:n_out) = fields(i,j,k,1:n_out)
end if
end do
end do
end do
! velocities
call xFFT3d(-1,1)
fname = 'u.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
call write_tmp4
call xFFT3d(-1,2)
fname = 'v.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
call write_tmp4
call xFFT3d(-1,3)
fname = 'w.'//file_ext
tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
call write_tmp4
! 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 io_write_4

414
m_fields.f90 Normal file
View file

@ -0,0 +1,414 @@
module m_fields
implicit none
real*8, allocatable :: fields(:,:,:,:)
!================================================================================
contains
!================================================================================
subroutine m_fields_init
use m_io
use m_parameters
implicit none
integer :: n
n = 3 + n_scalars + n_les
allocate(fields(nx+2,ny,nz,n), stat=ierr)
if (ierr.ne.0) then
write(out,*) "Cannot allocate fields, stopping."
call my_exit(-1)
end if
fields = zip
write(out,"('Allocated ',i3,' fields.')") n
call flush(out)
return
end subroutine m_fields_init
!================================================================================
subroutine m_fields_exit
use m_io
implicit none
if (allocated(fields)) deallocate(fields)
write(out,*) 'fields deallocated.'
call flush(out)
return
end subroutine m_fields_exit
!================================================================================
!--------------------------------------------------------------------------------
! Subroutine that broadcasts the array "fields" from the hydro part to the
! "stats" part of the code
!--------------------------------------------------------------------------------
subroutine fields_to_stats
use m_openmpi
use m_parameters
use m_io
implicit none
integer :: n_field, n_proc, k, ratio, n_scalars_bcast
! broadcasting time from the hydro root process to the whole world
!!$ write(out,*) "Broadcasting fields to stats part."
!!$ call flush(out)
! first send it to the root process of the stats part
count = 1
tag = 0
if (iammaster) then
if (task.eq.'hydro') call MPI_SEND(TIME,count,MPI_REAL8,id_root_stats,tag,MPI_COMM_WORLD,mpi_err)
if (task.eq.'stats') call MPI_RECV(TIME,count,MPI_REAL8,id_root_hydro,tag,MPI_COMM_WORLD,mpi_status,mpi_err)
!!$ write(out,*) "Exchanged information between master processors."
!!$ call flush(out)
end if
! then broadcast it over the "stats" communicator
if (task.eq.'stats') then
call MPI_BCAST(TIME,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
! checking if we need to start advancing scalars
if (.not. int_scalars .and. TIME .gt. TSCALAR) then
int_scalars = .true.
write(out,*) "Starting to move the scalars."
call flush(out)
end if
end if
!!$ write(out,*) "Broadcasted to slave processors."
!!$ call flush(out)
! figuring out how many scalars to broadcast:
! 0 if we do not move scalars
! all if we move scalars
n_scalars_bcast = 0
if (int_scalars) n_scalars_bcast = n_scalars
!!$ write(out,*) "Number of scalars to broadcast:", n_scalars_bcast
!!$ call flush(out)
if (numprocs_hydro .ge. numprocs_stats) then
! using the code structure:
! since the size of the arrays is always 2^n, there is always 2^k slabs
! of a "hydro" array that correspond to q slab of the "stats" array
ratio = numprocs_hydro / numprocs_stats
select case (task)
case ('hydro')
! sending the fields array to the corresponding process in "stats"
do n_field = 1,3+n_scalars_bcast
count = (nx+2)*ny*nz
id_to = numprocs_hydro + floor(real(myid_world) / real(ratio))
tag = myid_world*(3+n_scalars_bcast) + n_field-1
!!$ write(out,*) "Sending :", n_field, count, id_to, tag
!!$ call flush(out)
call MPI_ISEND(fields(1,1,1,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
!!$ write(out,*) "Sent."
!!$ call flush(out)
end do
case ('stats')
! receiving fields from hydro processors
do n_proc = 0,ratio-1
id_from = myid*ratio + n_proc
count = (nx+2)*ny*nz/ratio
k = (nz/ratio) * n_proc + 1
do n_field = 1,3+n_scalars_bcast
tag = id_from*(3+n_scalars_bcast) + n_field-1
!!$ write(out,*) "Receiving :", n_field, count, id_from, tag
!!$ call flush(out)
call MPI_IRECV(fields(1,1,k,n_field),count,MPI_REAL8,&
id_from,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
!!$ write(out,*) "Received."
!!$ call flush(out)
end do
end do
end select
else
! now doing the same for the case when more processors are involved in
! the "stat" part than in "hydro" part
ratio = numprocs_stats / numprocs_hydro
select case (task)
case ('hydro')
do n_proc = 0,ratio-1
count = (nx+2)*ny*nz/ratio
id_to = numprocs_hydro + myid_world*ratio + n_proc
k = (nz/ratio) * n_proc + 1
do n_field = 1,3+n_scalars_bcast
tag = id_to*(3+n_scalars_bcast) + n_field-1
call MPI_ISEND(fields(1,1,k,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
end do
end do
case ('stats')
do n_field = 1,3+n_scalars_bcast
count = (nx+2)*ny*nz
id_from = floor(real(myid) / real(ratio))
tag = myid_world*(3+n_scalars_bcast) + n_field-1
call MPI_RECV(fields(1,1,1,n_field),count,MPI_REAL8,&
id_from,tag,MPI_COMM_WORLD,mpi_status,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
end do
end select
end if
!!$ write(out,*) "broadcasted fields to stats"
!!$ call flush(out)
return
end subroutine fields_to_stats
!================================================================================
!================================================================================
!--------------------------------------------------------------------------------
! Subroutine that broadcasts the arrays that contain velocities in the x-space
! to the "parts" part of the code, for tracking particles
! The velocities are contained in the arrays wrk1...3
! The whole ideology remains similar to the subroutine fields_to_stats, except
! for the fact that the parts part of the code receives the velocities into
! the "fields" array.
!--------------------------------------------------------------------------------
subroutine fields_to_parts
use m_openmpi
use m_parameters
use m_io
use m_work
implicit none
integer :: n_field, n_proc, k, ratio, n_scalars_bcast
! if there are zero particles, return
if (nptot.eq.0) return
! broadcasting time and timestep (dt) from the hydro root process
! first send it to the root process of the "parts" part
count = 1
if (iammaster) then
tag = 0
if (task.eq.'hydro') call MPI_SEND(TIME,count,MPI_REAL8,id_root_parts,tag,MPI_COMM_WORLD,mpi_err)
if (task.eq.'parts') call MPI_RECV(TIME,count,MPI_REAL8,id_root_hydro,tag,MPI_COMM_WORLD,mpi_status,mpi_err)
tag = 1
if (task.eq.'hydro') call MPI_SEND(dt,count,MPI_REAL8,id_root_parts,tag,MPI_COMM_WORLD,mpi_err)
if (task.eq.'parts') call MPI_RECV(dt,count,MPI_REAL8,id_root_hydro,tag,MPI_COMM_WORLD,mpi_status,mpi_err)
end if
! then broadcast them over the "parts" communicator
if (task.eq.'parts') call MPI_BCAST(TIME,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
if (task.eq.'parts') call MPI_BCAST(dt ,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
! if have not yet started moving particles, return
if (TIME.lt.starttime_particles) return
! if this is the first timestep when we need to start moving particles,
! change the int_particle variable
if (TIME.ge.starttime_particles .and. .not. int_particles) then
write(out,*) 'Starting to move particles'
call flush(out)
int_particles = .true.
end if
! figure out if we broadcast scalars (currently not)
n_scalars_bcast = 0
! if (int_scalars) n_scalars_bcast = n_scalars
if (numprocs_hydro .ge. numprocs_parts) then
! using the code structure:
! since the size of the arrays is always 2^n, there is always 2^k slabs
! of a "hydro" array that correspond to q slab of the "parts" array
ratio = numprocs_hydro / numprocs_parts
select case (task)
case ('hydro')
! sending the fields array to the corresponding process in "stats"
do n_field = 1,3+n_scalars_bcast
count = (nx+2)*ny*nz
id_to = id_root_parts + floor(real(myid_world) / real(ratio))
tag = myid_world*(3+n_scalars_bcast) + n_field-1
! if the paricles are advected by fully resolved velocity
! ( that is, particles_filter_size=0) then send the fully resolved
! velocity to the "parts" task
! Else, if the particles are advected by locally averaged velofity, send
! the velocities in the Fourier form. They will be locally averaged and
! processed by the "parts" part of the code
if (particles_filter_size .le. 0.d0) then
call MPI_ISEND(wrk(1,1,1,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
else
call MPI_ISEND(fields(1,1,1,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
!!$
!!$ write(out,*) id_to, n_field, fields(:,1,1,n_field)
!!$ call flush(out)
!!$
end if
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
end do
case ('parts')
! receiving fields from hydro processors
do n_proc = 0,ratio-1
id_from = myid*ratio + n_proc
count = (nx+2)*ny*nz/ratio
k = (nz/ratio) * n_proc + 1
do n_field = 1,3+n_scalars_bcast
tag = id_from*(3+n_scalars_bcast) + n_field-1
call MPI_IRECV(fields(1,1,k,n_field),count,MPI_REAL8,&
id_from,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
!!$
!!$ write(out,*) 'rec',id_from, n_field, fields(:,1,k,n_field)
!!$ call flush(out)
!!$
end do
end do
end select
else
! now doing the same for the case when more processors are involved in
! the "parts" part than in "hydro" part
ratio = numprocs_parts / numprocs_hydro
select case (task)
case ('hydro')
do n_proc = 0,ratio-1
count = (nx+2)*ny*nz/ratio
id_to = id_root_parts + myid_world*ratio + n_proc
k = (nz/ratio) * n_proc + 1
do n_field = 1,3+n_scalars_bcast
tag = id_to*(3+n_scalars_bcast) + n_field-1
if (particles_filter_size .le. 0.d0) then
call MPI_ISEND(wrk(1,1,k,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
else
call MPI_ISEND(fields(1,1,k,n_field),count,MPI_REAL8,id_to,tag,MPI_COMM_WORLD,mpi_request,mpi_err)
end if
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
end do
end do
case ('parts')
do n_field = 1,3+n_scalars_bcast
count = (nx+2)*ny*nz
id_from = floor(real(myid) / real(ratio))
tag = myid_world*(3+n_scalars_bcast) + n_field-1
call MPI_RECV(fields(1,1,1,n_field),count,MPI_REAL8,&
id_from,tag,MPI_COMM_WORLD,mpi_status,mpi_err)
call MPI_WAIT(mpi_request,mpi_status,mpi_err)
end do
end select
end if
return
end subroutine fields_to_parts
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!!$
!!$ subroutine check_bcast
!!$
!!$ use m_parameters
!!$ use m_io
!!$ use m_work
!!$ implicit none
!!$
!!$ integer :: i,j,k,n
!!$ real*8 :: zyu
!!$
!!$ ! defining the fields to be cosines (testing)
!!$ zyu = 1.0
!!$ if (task.eq.'hydro') then
!!$ do n = 1,3
!!$ do k = 1,nz
!!$ do j = 1,ny
!!$ do i = 1,nx
!!$
!!$
!!$ fields(i,j,k,n) = sin(dble(n*i)*dx)
!!$
!!$ end do
!!$ end do
!!$ end do
!!$ end do
!!$ end if
!!$
!!$
!!$ call m_fields_bcast
!!$
!!$
!!$ do n = 1,3
!!$
!!$ if (task.eq.'hydro') write(fname,"('hydro',i1,'.arr')") n
!!$ if (task.eq.'stats') write(fname,"('stats',i1,'.arr')") n
!!$
!!$ tmp4(:,:,:) = fields(1:nx,:,:,n)
!!$ call write_tmp4
!!$ end do
!!$
!!$ return
!!$ end subroutine check_bcast
!!$
end module m_fields

289
m_filter_xfftw.f90 Normal file
View file

@ -0,0 +1,289 @@
!======================================================================
! Module that contains filtering rocedure for the part of the code that
! addvects lagrangian particles
!
! Time-stamp: <2009-05-20 11:50:52 (chumakov)>
!======================================================================
module m_filter_xfftw
use m_parameters
use x_fftw
use m_work
implicit none
! filter type (Gaussian by default)
integer :: filter_type = 2
! filter size
real*8 :: filter_size
! filter function in Fourier space
real*8, allocatable :: filter_g(:,:,:)
!======================================================================
contains
!======================================================================
!======================================================================
! filering subroutine for 2*pi^3-periodic domain
! IN: ss(n,:,:,:)
! OUT: ss(n:,;,:)
!======================================================================
subroutine filter_xfftw(n)
use m_io
implicit none
integer n
integer :: ix, iy, iz
real*8 :: a,b,c,d
!----------------------------------------------------------------------
!!$ call xFFT3d(1,n)
! multiplying wrk(:,:,:,n) by filter_g
! remember both are in the complex form
do iy = 1, nz
do iz = 1, nz_all
do ix = 1, nx + 1, 2
a = wrk(ix,iz,iy,n)
b = wrk(ix+1,iz,iy,n)
c = filter_g(ix,iz,iy)
d = filter_g(ix+1,iz,iy)
wrk(ix , iz, iy, n) = a*c - b*d
wrk(ix + 1, iz, iy, n) = b*c + a*d
end do
end do
end do
!!$ ! perform inverse FFT
!!$ call xFFT3d(-1,n)
return
end subroutine filter_xfftw
!======================================================================
! filering subroutine for 2*pi^3-periodic domain
! IN: ss(n,:,:,:)
! OUT: ss(n:,;,:)
!======================================================================
subroutine filter_xfftw_fields(n)
use m_io
use m_fields
implicit none
integer n
integer :: ix, iy, iz
real*8 :: a,b,c,d
!----------------------------------------------------------------------
!!$ call xFFT3d(1,n)
! multiplying wrk(:,:,:,n) by filter_g
! remember both are in the complex form
do iy = 1, nz
do iz = 1, nz_all
do ix = 1, nx + 1, 2
a = fields(ix,iz,iy,n)
b = fields(ix+1,iz,iy,n)
c = filter_g(ix,iz,iy)
d = filter_g(ix+1,iz,iy)
fields(ix , iz, iy, n) = a*c - b*d
fields(ix + 1, iz, iy, n) = b*c + a*d
end do
end do
end do
!!$ ! perform inverse FFT
!!$ call xFFT3d(-1,n)
return
end subroutine filter_xfftw_fields
!======================================================================
!======================================================================
!======================================================================
! subroutine that initializes the filter function filter_g
! IN: filter_type - type of the filter
! delta - filter's characteristic width
!
! OUT: filter_g - normalized FFT of the filter function
!======================================================================
subroutine filter_xfftw_init
use m_parameters
use m_io
use m_work
implicit none
real*8 delta
integer :: di,dj,dk,i,j,k, m
integer :: idelta,sum1
integer :: ii(8)
real*8 :: rx,ry,rz
real*8 :: const1,const2
filter_size = particles_filter_size
delta = filter_size
if (task.eq.'hydro') delta = dx*four
if (delta .lt. dx*three) then
write(out,*) "filter_xfftw_init: delta is too small:",delta,dx
write(out,*) "Must be at least 3*dx = ",three*dx
call flush(out)
call my_exit(-1)
end if
write(out,"('initializing the filter, filter_size = ',2e15.6)") delta, dx
call flush(out)
! the main idea is as follows:
! 1) create the filter kernel in one of the fields array
! 2) transform it into the Fourier space
! 3) put it into the array filter_g
! number of the slice in the fields array that will be used
m = LBOUND(fields,4)
!---------------------------------------------------------------
! define the filtering function
!---------------------------------------------------------------
case_filter_type: select case (filter_type)
!-------------------------------------------------------------------
! tophat filter
!-------------------------------------------------------------------
case (1)
write(out,*) '-- tophat filter, delta =',delta
write(out,*) "Tophat filter is not working at the moment."
write(out,*) "We're sorry for the inconvenience, stopping."
call my_exit(-1)
fields(:,:,:,m) = zip
idelta = delta / dx / 2
! normalization constant
const1 = real((2*idelta+1)**3,8)
const2 = 1.0d0 / const1
do k = 1,nz
dk = min(myid*nz+k-1,nz*numprocs-(myid*nz+k)+1)
do j = 1,ny
dj = min(j-1,ny-j+1)
do i = 1,nx
di = min(i-1,nx-i+1)
fields(i,j,k,m) = zip
if (di.le.idelta .and. dj.le.idelta .and. dk.le.idelta) then
fields(i,j,k,m) = const2
end if
end do
end do
end do
!-------------------------------------------------------------------
! Gaussian filter
!-------------------------------------------------------------------
case (2)
write(out,*) '-- Gaussian filter, delta =',delta
call flush(out)
fields(:,:,:,m) = zip
const1 = 6.0d0 / delta**2
const2 = sqrt(const1/PI)**3 *dx**3
do k = 1,nz
dk = min(myid*nz+k-1,nz*numprocs-(myid*nz+k)+1)
rz = dx * real(dk,8)
do j = 1,ny
dj = min(j-1,ny-j+1)
ry = dx * real(dj,8)
do i = 1,nx
di = min(i-1,nx-i+1)
rx = dx * real(di,8)
rx = rx*rx+ry*ry+rz*rz
fields(i,j,k,m) = const2*exp(-const1*rx)
if (fields(i,j,k,m) < 1.e-20) fields(i,j,k,m) = zip
end do
end do
end do
!-------------------------------------------------------------------
! linear filter
!-------------------------------------------------------------------
case(3)
write(out,*) '-- linear filter, delta =',delta
! if(myid.eq.0) print*,'-- linear filter, delta =',delta
const1 = 24.0d0 / (PI*delta**3) *dx**3
const2 = delta**2 / 4.0d0
do k = 1,nz
dk = min(myid*nz+k-1,nz*numprocs-(myid*nz+k)+1)
rz = dx * real(dk,8)
do j = 1,ny
dj = min(j-1,ny-j+1)
ry = dx * real(dj,8)
do i = 1,nx
di = min(i-1,nx-i+1)
rx = dx * real(di,8)
rx = rx*rx+ry*ry+rz*rz
fields(i,j,k,m) = zip
if (rx.le.const2) fields(i,j,k,m) = &
const1*(1.0d0-2.0d0*sqrt(rx)/delta)
end do
end do
end do
case default
print *,'FILTER_FFT_INIT: wrong filter_type:',filter_type
stop
end select case_filter_type
! outputting the sum of all elements of G.
! It should equal 1.0.
const1 = sum(fields(1:nx,1:ny,1:nz,m))
call MPI_REDUCE(const1,const2,1,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
if(myid.eq.0) write(out,*) '-- NORM OF G: ',const2
! compute FFT of g
call xFFT3d_fields(1,m)
! putting the FFT of the filter into filter_g.
! allocating the filter array
if (.not.allocated(filter_g)) then
allocate(filter_g(nx+2,ny,nz), stat=ierr)
if (ierr /= 0) stop 'cannot allocate filter_g'
filter_g = zip
write(out,*) 'allocated filter_g'
call flush(out)
end if
filter_g(:,:,:) = fields(:,:,:,m)
!!$! no need to normalize filter_g because our implementation of FFT
!!$! normalizes the result anyways
!!$ const1 = one/real(nx**3,8)
!!$ do k=1,nz
!!$ do j=1,ny
!!$ do i=1,nx+2
!!$ filter_g(i,j,k) = wrk(i,j,k,m)*const1
!!$ end do
!!$ end do
!!$ end do
return
end subroutine filter_xfftw_init
end module m_filter_xfftw

185
m_force.f90 Normal file
View file

@ -0,0 +1,185 @@
module m_force
use m_parameters, only : kfmax, famp, force_type
integer*4 :: n_forced_nodes, n_forced_nodes_total
! coordinated of the forced nodes
integer, allocatable :: ifn(:), jfn(:), kfn(:), k_shell(:)
contains
subroutine m_force_init
use m_parameters
use x_fftw
implicit none
integer :: i, j, k, n, n_shell
! if flow is not forced, return
if (flow_type .ne. 1) return
if (task.ne.'hydro') return
select case (force_type)
case (1)
! Machiels forcing (see article in PRL #79(18) p.3411)
write(out,*) "Forcing #1: Machiels forcing - setting up"
call flush(out)
! find out how many nodes are we forcing and book them
n_forced_nodes = 0
n_forced_nodes_total = 0
do k = 1,nz
do j = 1,ny
do i = 1,nx
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kfmax) then
n_forced_nodes = n_forced_nodes + 1
end if
end do
end do
end do
! reducing to the master process to find out the total number of forced nodes
!!$ write(out,*) 'before reducing'; call flush(out)
count = 1
call MPI_REDUCE(n_forced_nodes, n_forced_nodes_total, count, &
MPI_INTEGER4,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! writing out the # of forced nodes
write(out,*) 'Number of forced nodes for this process:',n_forced_nodes
if (myid.eq.0) write(out,*) ' total number:',n_forced_nodes_total
call flush(out)
! allocating arrays for the coordinates of the forced nodes
allocate(ifn(n_forced_nodes), jfn(n_forced_nodes), kfn(n_forced_nodes), &
k_shell(n_forced_nodes))
! filling up the arrays
n = 0
do k = 1,nz
do j = 1,ny
do i = 1,nx
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kfmax) then
n = n + 1
ifn(n) = i
jfn(n) = j
kfn(n) = k
k_shell(n) = n_shell
end if
end do
end do
end do
!!$ ! writing out the nodes
!!$ do n = 1,n_forced_nodes
!!$ write(out,"(3i4)") ifn(n),jfn(n),kfn(n)
!!$ end do
!!$ call flush(out)
case default
write(out,*) 'WRONG FORCE TYPE:',force_type
write(out,*) 'STOPPING'
call flush(out)
stop
end select
return
end subroutine m_force_init
!================================================================================
subroutine force_velocity
! adding forcing to the arrays wrk(:,:,:,1:3) that already contain the RHS for velocities
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_stats
implicit none
integer :: i, j, k, n_shell, n
real*8 :: fac, fac2
select case (force_type)
case (1)
! Machiels forcing (see article in PRL #79(18) p.3411)
! write(out,*) "Machiels forcing"; call flush(out)
e_spec = zip
e_spec1 = zip
hits = 0
hits1 = 0
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
! assembling the total energy in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
if (n_shell .gt. 0 .and. n_shell .le. kfmax) then
fac2 = fac * (fields(i,j,k,1)**2 + fields(i,j,k,2)**2 + fields(i,j,k,3)**2)
if (akx(i).eq.0.d0) fac2 = fac2 * 0.5d0
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kfmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! getting the total energy in the region [0:kfmax] by integrating the spectrum
if (myid.eq.0) energy = sum(e_spec(1:kfmax))
! broadcasting the current energy in the forcing range of wavenumbers
count = 1
call MPI_BCAST(energy,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
! now applying the forcing to the RHS for velocities (wrk(:,:,:,1:3))
fac = FAMP / energy
do n = 1,n_forced_nodes
n_shell = k_shell(n)
i = ifn(n)
j = jfn(n)
k = kfn(n)
wrk(i,j,k,1) = wrk(i,j,k,1) + fac * fields(i,j,k,1)
wrk(i,j,k,2) = wrk(i,j,k,2) + fac * fields(i,j,k,2)
wrk(i,j,k,3) = wrk(i,j,k,3) + fac * fields(i,j,k,3)
end do
case default
write(out,*) "WRONG FORCE_TYPE:",force_type
stop
end select
return
end subroutine force_velocity
end module m_force

53
m_io.f90 Normal file
View file

@ -0,0 +1,53 @@
!================================================================================
! M_IO - module that contains the I/O related subroutines.
!
! Time-stamp: <2008-11-04 15:31:31 (chumakov)>
!================================================================================
module m_io
use m_openmpi
implicit none
character*6 :: file_ext
character*80 :: fname
! output file handle
integer :: in=10, out=11
!================================================================================
contains
!================================================================================
!================================================================================
subroutine m_io_init
implicit none
write(fname,"('d',i4.4,'.txt')") myid_world
open(out,file=fname,position="append")
write(out,"('-------------------------------------')")
write(out,"('Process ',i4,' of ',i4,'(',i4.4,') is alive.')") &
myid_world, numprocs_world, numprocs_world-1
write(out,"('My task is ""',a5,'"", my id is',i4)") task, myid
call flush(out)
return
end subroutine m_io_init
!================================================================================
!================================================================================
subroutine m_io_exit
implicit none
write(out,"('Done.')")
close(out)
return
end subroutine m_io_exit
!================================================================================
!================================================================================
!================================================================================
!================================================================================
end module m_io

2396
m_les.f90 Normal file
View file

File diff suppressed because it is too large Load diff

250
m_openmpi.f90 Normal file
View file

@ -0,0 +1,250 @@
!================================================================================
! Module contains interface to OpenMPI
!
! Time-stamp: <2009-08-20 14:22:13 (chumakov)>
!================================================================================
module m_openmpi
!================================================================================
implicit none
include 'mpif.h'
! Uncomment this for the systems that do not have OpenMPI
! In OpenMPI, the parameter MPI_INTEGER_KIND is defined in 'mpif.h'
! With other MPI implementations, this parameter has to be defined manually.
integer MPI_INTEGER_KIND
parameter (MPI_INTEGER_KIND = 4)
! --- MPI variables
logical :: iammaster
integer(kind=MPI_INTEGER_KIND) :: myid_world, numprocs_world
integer(kind=MPI_INTEGER_KIND) :: numprocs_hydro, numprocs_stats, numprocs_parts
integer(kind=MPI_INTEGER_KIND) :: myid, numprocs, master, mpi_err, mpi_info
integer(kind=MPI_INTEGER_KIND) :: id_to, id_from, tag, count
integer(kind=MPI_INTEGER_KIND) :: id_root_hydro, id_root_stats, id_root_parts
! communicator for separate tasks
integer(kind=MPI_INTEGER_KIND) :: MPI_COMM_TASK
! exclusive communicator for root processes of tasks
integer(kind=MPI_INTEGER_KIND) :: MPI_COMM_ROOTS
integer (kind=MPI_INTEGER_KIND) :: sendtag, recvtag
integer (kind=MPI_INTEGER_KIND) :: request, request1, request2, request3, mpi_request
integer (kind=MPI_INTEGER_KIND) :: id_l, id_r
integer (kind=mpi_INTEGER_KIND) :: mpi_status(MPI_STATUS_SIZE)
integer(kind=MPI_INTEGER_KIND) :: color, key
character*5 :: task, split="nevah"
character*10 :: run_name_local
logical :: task_split=.false.
!================================================================================
contains
!================================================================================
subroutine m_openmpi_init
implicit none
integer (kind=mpi_INTEGER_KIND) :: n
integer*4 :: np_local
integer :: i
! first getting the run name form the command line
! (it's local, not global run_name)
! also getting the parameter "split" which governs the process splitting:
! split="split" means that hydro, statistics and particles are assigned three
! separate process groups (they differ by the char*5 parameter "task").
! split="never" (default if the parameter is missing) means that all
! processes do all tasks. (does not work for the particles at this point)
call openmpi_get_command_line
! initializing MPI environment
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)
!--------------------------------------------------------------------------------
! Looking at the command line parameter called "split". If it equals "split"
! then we define task_split=.true. If not, task_split remains .false. (default)
!--------------------------------------------------------------------------------
if (split == "split") task_split = .true.
!--------------------------------------------------------------------------------
! First check if we need to do any task splitting. If we don't (split="never")
! then we define task="hydro" and do a ficticious split with uniform color of
! all processors.
!--------------------------------------------------------------------------------
if (.not. task_split) then
!!$ print *,'not splitting into task groups, all procs are "hydro"'
task = 'hydro'
color = 0
myid = myid_world
goto 1000
end if
!--------------------------------------------------------------------------------
! Definition of processor groups: hydro, stats, parts etc. for task splitting.
!--------------------------------------------------------------------------------
! first finding out if there are any particles involved.
! if there are no particles, then we split the processors in two parts:
! hydro and stats. If there are some particles, we split the processors
! in three parts: "hydro", "stats" and "parts". The variables that
! determines which part the process belongs to is "task".
! first see, how many particles are there
if (myid_world.eq.0) then
! opening the inupt file
open(99,file=run_name_local//'.in')
! skipping the first 35 lines
do i = 1,35
read(99,*)
end do
! reading the number of particles
read(99,*) np_local
close(99)
end if
! broadcasting the number of particles to all processors
count = 1
call MPI_BCAST(np_local,count,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
! now splitting the processors in tasks: hydro, stats and parts
! the curren logic is this:
! - In case if there are no particles, the split between the hydro and
! the stats part is 2/3 and 1/3. This way the total number of processors
! needs to be 3*2^n
! - In case with the particles in the flow, the split is 1/2, 1/4 and 1/4
! first the case when we do not have particles
if (np_local.eq.0) then
! if the numprocs_total is divisible by 3, assign 2/3 of it to hydro
! and the rest to stats
if (int(numprocs_world/3)*3 .eq. numprocs_world) then
numprocs_hydro = numprocs_world * 2/3
numprocs_stats = numprocs_world - numprocs_hydro
numprocs_parts = 0
else if (2**floor(log(real(numprocs_world))/log(2.d0)) .eq. numprocs_world) then
print*, 'numprocs_world is 2^n, allocating half for hydro: ',numprocs_world
numprocs_hydro = numprocs_world / 2
numprocs_stats = numprocs_world - numprocs_hydro
numprocs_parts = 0
else
! if the # of processors N is not 2^n and not divisible by 3, then just take
! the biggest 2^k < N and make these hydro, the rest - stat.
numprocs_hydro = 2**floor(log(real(numprocs_world))/log(2.d0))
numprocs_stats = numprocs_world - numprocs_hydro
numprocs_parts = 0
end if
id_root_hydro = 0
id_root_stats = numprocs_hydro
id_root_parts = 0
else
numprocs_hydro = numprocs_world / 2
numprocs_stats = numprocs_world / 4
numprocs_parts = numprocs_world / 4
id_root_hydro = 0
id_root_stats = numprocs_hydro
id_root_parts = numprocs_hydro + numprocs_stats
end if
! splitting the communicator into several parts
! 1. hydro
! 2. stats
! 3. parts
if (myid_world.lt.numprocs_hydro) then
task = 'hydro'
color = 0
myid = myid_world
elseif (myid_world.ge.numprocs_hydro .and. myid_world .lt. numprocs_hydro+numprocs_stats) then
task = 'stats'
color = 1
myid = myid_world - numprocs_hydro
else
task = 'parts'
color = 2
myid = myid_world - numprocs_hydro - numprocs_stats
end if
!--------------------------------------------------------------------------------
! The actual task splitting happens here
!--------------------------------------------------------------------------------
1000 continue
call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,myid,MPI_COMM_TASK,mpi_err)
call MPI_COMM_SIZE(MPI_COMM_TASK,numprocs,mpi_err)
call MPI_COMM_RANK(MPI_COMM_TASK,myid,mpi_err)
! each task will have its master process
master = 0
iammaster = .false.
if (myid.eq.master) iammaster=.true.
!!$ ! The following is put on hold because it looks like a crazy idea
!!$ ! now creating separate exclusive communicator for the master nodes only
!!$ ! the name of the new communicator is MPI_COMM_ROOTS
!!$ ! if we want quickly broadcast something, then we can use two BCAST calls
!!$ color = 1
!!$ if (iammaster) color = 0
!!$ call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,myid_world,MPI_COMM_ROOTS,mpi_err)
return
end subroutine m_openmpi_init
!================================================================================
subroutine m_openmpi_exit
call MPI_COMM_FREE(MPI_COMM_TASK,mpi_err)
call MPI_FINALIZE(mpi_err)
return
end subroutine m_openmpi_exit
!================================================================================
subroutine openmpi_get_command_line
implicit none
character*80 :: tmp_str
integer :: iargc
! reading the run_name from the command line
if(iargc().eq.0) then
call getarg(0,tmp_str)
print*, 'Format: ',trim(tmp_str),' (run name) ["split"/"never"]'
stop
end if
call getarg(1,run_name_local)
if(len_trim(run_name_local).ne.10) then
print *, 'Run name: "',run_name_local,'"'
print *, ' "1234567890"'
print *, 'Length of run name is less than 10, sorry.'
stop
end if
! getting the split parameter, if it's there
if(iargc().eq.2) call getarg(2,split)
end subroutine openmpi_get_command_line
!================================================================================
end module m_openmpi

515
m_parameters.f90 Normal file
View file

@ -0,0 +1,515 @@
!================================================================================
! M_PARAMETERS - module for all parameters in the calculation:
! such as array dimensions, reynolds numbers, switches/flags etc.
!
! Time-stamp: <2009-08-19 11:50:11 (chumakov)>
! Time-stamp: <2008-11-20 17:27:59 MST (vladimirova)>
!================================================================================
module m_parameters
use m_openmpi
use m_io
implicit none
! --- problem related
character*10 :: run_name
! --- input filep arameters
integer :: nx,ny,nz, nz_all ! Dimensions of the problem
integer :: nxyz, nxyz_all
integer :: n_scalars ! # of scalars
real*8 :: time ! time of simulation
real*8 :: dx, dy, dz
integer :: kmax
integer :: ITIME, ITMIN, ITMAX, IPRINT1, IPRINT2, IWRITE4
real*8 :: TMAX, TRESCALE, TSCALAR, RE, nu, dt
! now many times to rescale teh velocities
integer :: NRESCALE
integer :: flow_type
logical :: variable_dt
integer :: isp_type, ir_exp, force_type
real*8 :: peak_wavenum
real*8 :: famp
integer :: kfmax ! Maximum wavenumber for forcing (integer)
real*8 :: courant
integer :: dealias
integer :: det_rand
real*8 :: RN1, RN2, RN3
! particle-related
! indicator that says which particle tracking scheme to use:
! 0 = trilinear
! 1 = spectral
! 2 = tricubic
! trilinear by default
integer :: particles_tracking_scheme = 0
real*8 :: starttime_particles
! sometimes we want to advect particles by locally averaged field
! the following variables address that concern
real*8 :: particles_filter_size
! number of particles assigned to the processor
! and the total number of particles
integer(kind=MPI_INTEGER_KIND) :: np, np1, nptot
! If using Large Eddy Simulation (LES), the LES model ID is here
integer :: les_model
integer, allocatable :: scalar_type(:)
real*8, allocatable :: pe(:), sc(:), ir_exp_sc(:), peak_wavenum_sc(:), reac_sc(:)
! constants
real*8 :: zip=0.0d0, half=0.5d0
real*8 :: one=1.0d0,two=2.0d0,three=3.d0,four=4.d0,five=5.d0, six=6.d0
integer :: last_dump
! --- supporting stuff
logical :: there
logical :: fos, fov
integer :: ierr
real*8 :: PI, TWO_PI
logical :: int_scalars, int_particles
! --- number of LES variables in the arrays (initialized to zero)
integer :: n_les = 0
! benchmarking tools
logical :: benchmarking=.false.
integer (kind=8) :: i81, i82, bm(12)
!================================================================================
contains
!================================================================================
subroutine m_parameters_init
implicit none
call get_run_name
! constants
PI = four * atan(one)
TWO_PI = two * PI
! switches
int_scalars = .false.
call read_input_file
! maximum resolved wavenumber
if (dealias.eq.0) then
kmax = nx/3
elseif (dealias.eq.1) then
kmax = floor(real(nx,8) / three * sqrt(two))
else
write(out,*) "*** M_PARAMETERS_INIT: wrong dealias flag: ",dealias
call flush(out)
call my_exit(-1)
end if
write(out,*) "kmax = ",kmax
call flush(out)
end subroutine m_parameters_init
!================================================================================
subroutine get_run_name
implicit none
character*80 :: tmp_str
integer :: iargc
! reading the run_name from the command line
if(iargc().eq.0) then
call getarg(0,tmp_str)
write(out,*) 'Format: ',trim(tmp_str),' <run name>'
write(*,*) 'Format: ',trim(tmp_str),' <run name>'
call flush(out)
call MPI_FINALIZE(ierr)
stop
end if
call getarg(1,run_name)
if(len_trim(run_name).ne.10) then
write(out,*) 'Run name: "',run_name,'"'
write(out,*) ' "1234567890"'
write(out,*) 'Length of run name is less than 10, sorry.'
call MPI_FINALIZE(ierr)
stop
end if
write(out,*) 'Run name: "',run_name,'"'
call flush(out)
end subroutine get_run_name
!================================================================================
subroutine read_input_file
implicit none
logical :: there
integer :: n
integer*4 :: passed, passed_all
character*80 :: str_tmp
! making sure the input file is there
inquire(file=run_name//'.in', exist=there)
if(.not.there) then
write(out,*) '*** cannot find the input file'
call flush(out)
call my_exit(-1)
end if
! now the variable "passed" will show if the parameters make sense
passed = 1
! -------------------------------------------------
! reading parameters from the input file
! and checking them for consistency
! -------------------------------------------------
open(in,file=run_name//'.in',form='formatted')
read(in,*)
read(in,*)
read(in,*)
read(in,*,ERR=9000) nx,ny,nz_all
read(in,*)
nz = nz_all/numprocs
if (nz*numprocs.ne.nz_all) then
write(out,*) '*** wrong nz_all:', nz_all, &
'*** should be divisible by numprocs:',numprocs
call flush(out)
passed = 0
end if
write(out,'(70(''=''))')
write(out,"('NX,NY,NZ_ALL', 3i4)") nx,ny,nz_all
write(out,"('NX,NY,NZ ', 3i4)") nx,ny,nz
call flush(out)
dx = 2.0d0 * PI / dble(nx)
dy = 2.0d0 * PI / dble(ny)
dz = 2.0d0 * PI / dble(nz_all)
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) ITMIN
write(out,*) 'ITMIN = ',ITMIN
last_dump = ITMIN
read(in,*,ERR=9000,END=9000) ITMAX
write(out,*) 'ITMAX = ',ITMAX
read(in,*,ERR=9000,END=9000) IPRINT1
write(out,*) 'IPRINT1= ',IPRINT1
read(in,*,ERR=9000,END=9000) IPRINT2
write(out,*) 'IPRINT2= ',IPRINT2
read(in,*,ERR=9000,END=9000) IWRITE4
write(out,*) 'IWRITE4= ',IWRITE4
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) TMAX
write(out,*) 'TMAX =',TMAX
read(in,*,ERR=8000,END=9000) TRESCALE, NRESCALE
100 write(out,*) 'TRESCALE, NRESCALE =',TRESCALE, NRESCALE
read(in,*,ERR=9000,END=9000) TSCALAR
write(out,*) 'TSCALAR =',TSCALAR
read(in,*)
write(out,"(70('-'))")
call flush(out)
! if(TSCALAR.le.TRESCALE) then
! TSCALAR = TRESCALE
! write(out,*) '*** RESET: TSCALAR = ',TSCALAR
! end if
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) flow_type
write(out,*) 'flow_type ', flow_type
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) RE
write(out,*) 'RE = ',RE
nu = 1.0d0/RE
read(in,*,ERR=9000,END=9000) DT
write(out,*) 'DT = ',DT
if (dt.lt.0.0d0) then
variable_dt = .false.
dt = -dt
else
variable_dt = .true.
end if
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) isp_type
write(out,*) 'isp_type= ', isp_type
read(in,*,ERR=9000,END=9000) ir_exp
write(out,*) 'ir_exp = ', ir_exp
read(in,*,ERR=9000,END=9000) peak_wavenum
write(out,*) 'peak_wavenum = ',peak_wavenum
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) force_type
write(out,*) 'force_type', force_type
read(in,*,ERR=9000,END=9000) kfmax
write(out,*) 'kfmax = ',kfmax
read(in,*,ERR=9000,END=9000) FAMP
write(out,*) 'FAMP = ',FAMP
read(in,*)
write(out,"(70('-'))")
call flush(out)
!!$ c------------------------------------------------------------
!!$
!!$ read(in,*,ERR=9000,END=9000) IRESET
!!$ write(out,*) 'IRESET= ',IRESET
!!$
!!$ read(in,*,ERR=9000,END=9000) INEWSC
!!$ write(out,*) 'INEWSC= ',INEWSC
!!$ read(in,*)
!!$
!!$ c------------------------------------------------------------
read(in,*,ERR=9000,END=9000) dealias
write(out,*) 'dealias = ',dealias
read(in,*)
write(out,"(70('-'))")
call flush(out)
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) det_rand
write(out,*) 'det_rand =',det_rand
read(in,*,ERR=9000,END=9000) RN1
write(out,*) 'RN1 =',RN1
read(in,*,ERR=9000,END=9000) RN2
write(out,*) 'RN2 =',RN2
read(in,*,ERR=9000,END=9000) RN3
write(out,*) 'RN3 =',RN3
read(in,*)
write(out,"(70('-'))")
call flush(out)
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) nptot
! DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG
if (.not.task_split .and. nptot > 0) then
write(out,*) "tasks are not split, making nptot=0"
nptot = 0
end if
! DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG
write(out,*) 'nptot =',nptot
read(in,*,ERR=9000,END=9000) particles_tracking_scheme
write(out,*) 'particles_tracking_scheme', particles_tracking_scheme
select case (particles_tracking_scheme)
case (0)
write(out,*) '--- Trilinear tracking'
case (1)
write(out,*) '--- CINT (cubic interpolation on integer nodes)'
case (2)
write(out,*) '--- Spectral tracking (CAUTION: SLOW!)'
case default
write(out,*) 'don''t recognize particle tracking:', &
particles_tracking_scheme
write(out,*) 'reset to zero'
particles_tracking_scheme = 0
end select
call flush(out)
! DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG
if (particles_tracking_scheme .gt. 1) stop 'Cannot do this particle tracking'
! DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG DEBUG
read(in,*,ERR=9000,END=9000) starttime_particles
write(out,*) 'starttime_particles: ',starttime_particles
read(in,*,ERR=9000,END=9000) particles_filter_size
write(out,*) 'particles_filter_size:',particles_filter_size
if (particles_filter_size .gt. zip .and. particles_filter_size .lt. three*dx) then
write(out,*) "particles_filter_size is too small (less than 3*dx)"
write(out,*) particles_filter_size, three*dx
call flush(out)
call my_exit(-1)
end if
read(in,*)
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) les_model
write(out,*) 'les_model =',les_model
read(in,*)
write(out,"(70('-'))")
call flush(out)
! making sure that if the LES mode is on, the dealiasing is 3/2-rule
if (les_model .gt. 0 .and. dealias .ne. 0) then
dealias = 0
write(out,*) "*** LES mode, changing dealias to 0."
call flush(out)
end if
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) n_scalars
write(out,*) '# of scalars:', n_scalars
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
! if there are scalars, then read them one by one
if (n_scalars>0) then
read(in,'(A)',ERR=9000,END=9000) str_tmp
write(out,*) str_tmp
call flush(out)
! reading parameters of each scalar
allocate(scalar_type(n_scalars), pe(n_scalars), sc(n_scalars), &
ir_exp_sc(n_scalars), peak_wavenum_sc(n_scalars), &
reac_sc(n_scalars), stat=ierr)
if (ierr.ne.0) passed = 0
do n = 1,n_scalars
read(in,*,ERR=9000,END=9000) scalar_type(n), sc(n), ir_exp_sc(n), &
peak_wavenum_sc(n), reac_sc(n)
write(out,'(9x,i4,1x,4(f8.3,1x))') scalar_type(n), sc(n), ir_exp_sc(n), &
peak_wavenum_sc(n), reac_sc(n)
call flush(out)
PE(n) = nu/SC(n) ! INVERSE Peclet number
end do
end if
! -------------------------------------------------------------
! closing the input file
close(in)
write(out,'(70(''=''))')
call flush(out)
! defining the rest of the parameters
nxyz = nx * ny * nz
nxyz_all = nx * ny * nz_all
! ------------------------------------------------------------
!--------------------------------------------------------------------------------
! Checking if the task splitting conflicts with particle advection. Currently
! we canot have split=never and have particles. This is to be resolved later,
! now my head is spinning already.
!--------------------------------------------------------------------------------
if (.not.task_split .and. nptot.gt.0) then
write(out,*) "*** READ_INPUT_FILE: Cannot have .not.task_split and nptot > 0. Stopping"
call flush(out)
passed = 0
end if
!--------------------------------------------------------------------------------
count = 1
call MPI_REDUCE(passed,passed_all,count,MPI_INTEGER4,MPI_MIN,0,MPI_COMM_WORLD,mpi_err)
count = 1
call MPI_BCAST(passed_all,count,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
if (passed.lt.one) then
write(out,*) "not passed the check, stopping"
call flush(out)
stop
end if
return
!--------------------------------------------------------------------------------
! ERROR PROCESSING
!--------------------------------------------------------------------------------
8000 continue
NRESCALE = 0
if (TRESCALE.gt.zip) NRESCALE = 1
write(out,*) "*** NRESCALE IS AUTOMATICALLY ASSIGNED to be ONE"
call flush(out)
goto 100
9000 continue
write(out,*)'An error was encountered while reading input file'
call flush(out)
stop
end subroutine read_input_file
!================================================================================
end module m_parameters

1512
m_particles.f90 Normal file
View file

File diff suppressed because it is too large Load diff

155
m_rand_knuth.f90 Normal file
View file

@ -0,0 +1,155 @@
MODULE m_rand_knuth
! Random number generator
! Code converted using TO_F90 by Alan Miller
! Date: 2000-09-10 Time: 16:37:48
! Latest revision - 16 January 2003
! FORTRAN 77 version of "ran_array"
! from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
! including the MODIFICATIONS made in the 9th printing (2002)
! ********* see the book for explanations and caveats! *********
! Author: Steve Kifowit
! http://ourworld.compuserve.com/homepages/steve_kifowit
! with modifications by Alan Miller to rnarry and rnstrt based upon
! Knuth's code.
! For Donald Knuth's Fortran 77 versions, go to:
! http://www-cs-faculty.stanford.edu/~knuth/programs
! Look for frng.f and frngdb.f
IMPLICIT NONE
INTEGER, PARAMETER :: kk=100, ll=37, mm=2**30, tt=70, kkk=kk+kk-1
INTEGER, SAVE :: ranx(kk)
CONTAINS
SUBROUTINE rand_knuth(u, n)
! Generate an array of n real values between 0 and 1.
REAL , INTENT(OUT) :: u(:)
INTEGER, INTENT(IN) :: n
! Local array
INTEGER :: aa(n)
CALL rnarry(aa, n)
u(1:n) = SCALE( REAL(aa), -30)
RETURN
END SUBROUTINE rand_knuth
SUBROUTINE rnarry(aa, n)
! Generate an array of n integers between 0 and 2^30-1.
INTEGER, INTENT(OUT) :: aa(:)
INTEGER, INTENT(IN) :: n
! Local variables
INTEGER :: j
aa(1:kk) = ranx(1:kk)
DO j = kk + 1, n
aa(j) = aa(j-kk) - aa(j-ll)
IF (aa(j) < 0) aa(j) = aa(j) + mm
END DO
DO j=1,ll
ranx(j) = aa(n+j-kk) - aa(n+j-ll)
IF (ranx(j) < 0) ranx(j) = ranx(j) + mm
END DO
DO j=ll+1,kk
ranx(j) = aa(n+j-kk) - ranx(j-ll)
IF (ranx(j) < 0) ranx(j) = ranx(j) + mm
END DO
RETURN
END SUBROUTINE rnarry
SUBROUTINE rand_knuth_start(seed)
! Initialize integer array ranx using the input seed.
INTEGER*8, INTENT(IN) :: seed
! Local variables
INTEGER :: x(kkk), j, ss, sseed, t
IF (seed < 0) THEN
sseed = mm - 1 - MOD(-1-seed, mm)
ELSE
sseed = MOD(seed, mm)
END IF
ss = sseed - MOD(sseed,2) + 2
DO j=1, kk
x(j) = ss
ss = ISHFT(ss, 1)
IF (ss >= mm) ss = ss - mm + 2
END DO
x(kk+1:kkk) = 0
x(2) = x(2)+1
ss = sseed
t = tt - 1
10 DO j=kk, 2, -1
x(j+j-1) = x(j)
END DO
DO j = kkk, kk + 1, -1
x(j-(kk-ll)) = x(j-(kk-ll)) - x(j)
IF (x(j-(kk-ll)) < 0) x(j-(kk-ll)) = x(j-(kk-ll)) + mm
x(j-kk) = x(j-kk) - x(j)
IF (x(j-kk) < 0) x(j-kk) = x(j-kk) + mm
END DO
IF (MOD(ss,2) == 1) THEN
DO j=kk, 1, -1
x(j+1) = x(j)
END DO
x(1) = x(kk+1)
x(ll+1) = x(ll+1) - x(kk+1)
IF (x(ll+1) < 0) x(ll+1) = x(ll+1) + mm
END IF
IF (ss /= 0) THEN
ss = ISHFT(ss, -1)
ELSE
t = t - 1
END IF
IF (t > 0) GO TO 10
DO j=1, ll
ranx(j+kk-ll) = x(j)
END DO
DO j=ll+1,kk
ranx(j-ll) = x(j)
END DO
DO j = 1, 10
CALL rnarry(x,kkk)
END DO
RETURN
END SUBROUTINE rand_knuth_start
! Initialization subroutine
subroutine rand_knuth_init(seed1)
integer*8 :: seed1
call rand_knuth_start(seed1)
return
end subroutine rand_knuth_init
END MODULE m_rand_knuth

445
m_stats.f90 Normal file
View file

@ -0,0 +1,445 @@
module m_stats
use m_parameters
real*8, allocatable :: e_spec(:), e_spec1(:), moments(:,:), sc_diss(:), sc_min(:), sc_max(:)
integer*8, allocatable :: hits(:), hits1(:)
real*8 :: energy, eps_v, eta, etakmax, enstrophy, re_lambda, uvar, x_length
real*8 :: lambda, re_lambda1, tau_e
real*8 :: sctmp
contains
!================================================================================
! This is a small module so the arrays get allocated in all parts of the code
!================================================================================
subroutine m_stats_init
implicit none
allocate(e_spec(kmax), e_spec1(kmax), moments(3+n_scalars,4), hits(kmax), hits1(kmax),&
sc_diss(n_scalars), sc_min(n_scalars), sc_max(n_scalars),stat=ierr)
write(out,*) "Stat allocated.", ierr
e_spec = zip
e_spec1 = zip
hits = 0
hits1 = 0
sc_diss = zip
sc_min = zip
sc_max = zip
return
end subroutine m_stats_init
!================================================================================
!================================================================================
subroutine stat_main
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: n
logical :: there2
real*8 :: fac, fac2
call stat_velocity
if (int_scalars) then
call stat_scalars
! now outputting the scalar statistics
do n = 1, n_scalars
if (myid.eq.0) then
write(fname,"('sc',i2.2,'.gp')") n
inquire(file=fname, exist=there, opened=there2)
if (.not.there) then
open(100+n,file=fname,form='formatted')
write(100+n,'(A)') '# 1.itime 2.time 3.sc.diss 4. mean 5.variance 6.min 7.max'
end if
if(there.and..not.there2) then
open(100+n,file=fname,position='append')
end if
write(100+n,"(i7,10e15.6)") itime, time, sc_diss(n), moments(3+n,1:2), sc_min(n), sc_max(n)
call flush(100+n)
end if
end do
end if
if (task_split) then
write(out,9000) ITIME, TIME
call flush(out)
end if
return
9000 format('ITIME=',i7,' TIME=',e15.7,' Stat files are written.')
end subroutine stat_main
!================================================================================
subroutine stat_velocity
implicit none
logical :: there2
integer :: k, n
! getting the enstrophy
call get_gradient_statistics
! getting the energy spectrum e_spec to the main process
call get_e_spec
! outputting the statistics into files
if (myid.eq.0) then
! getting the total energy
energy = sum(e_spec(1:kmax))
! finding dissipation spectrum and total dissipation
do k = 1,kmax
e_spec1(k) = e_spec(k) * real(k**2,8) * two * nu
end do
eps_v = sum(e_spec1(1:kmax))
! finding Kolmogorov scale
eta = (nu**3/eps_v)**0.25
etakmax = eta * real(kmax,8)
! variance
uvar = two/three*energy
! integral length scale
sctmp = zip
do k = 1, kmax
sctmp = sctmp + e_spec(k) / real(k,8)
end do
x_length = PI / two * sctmp / uvar
! Taylor microscale
lambda = sqrt(15.d0 * uvar * nu / eps_v)
! Taylor-Reynolds number
re_lambda = uvar*sqrt(15.d0/eps_v*RE)
re_lambda1 = sqrt(uvar)*lambda / nu
! Eddy turnover time
tau_e = x_length / sqrt(uvar)
! outputting all this in the stat1 file
inquire(file='stat1.gp', exist=there, opened=there2)
if (.not.there) then
open(69,file='stat1.gp',form='formatted')
write(69,*) '# 1.itime 2.time 3.energy 4.diss 5.eta 6.enstrophy 7.R_lambda'
end if
if(there.and..not.there2) then
open(69,file='stat1.gp',position='append')
end if
write(69,"(i8,20e15.6)") itime, time, energy, eps_v, eta, enstrophy, re_lambda
call flush(69)
! outputting all this in the stat2 file
inquire(file='stat2.gp', exist=there, opened=there2)
if (.not.there) then
open(70,file='stat2.gp',form='formatted')
write(70,'(A)') '# 1.itime 2.time 3.int LS 4. lambda 5.R_lambda1 6.tau_e 7.etakmax'
end if
if(there.and..not.there2) then
open(70,file='stat2.gp',position='append')
end if
write(70,"(i8,20e15.6)") itime, time, x_length, lambda, re_lambda1, tau_e, etakmax
call flush(70)
! outputting the energy spectrum
open(900,file='es.gp',position='append')
write(900,"()")
write(900,"()")
write(900,"('# ITIME=',i7,' TIME=',e17.8)") ITIME, TIME
do k = 1,kmax !min(kmax,nx/3)
write(900,"(i4,2e15.6)") k,e_spec(k), e_spec1(k)
end do
close(900)
end if
return
end subroutine stat_velocity
!================================================================================
!================================================================================
!================================================================================
subroutine get_e_spec
use m_io
use m_fields
use x_fftw
implicit none
real*8 :: sc_rad1, sc_rad2, fac, fac2
integer :: i, j, k, n_shell
real*8 :: energy2
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
e_spec1 = zip
e_spec = zip
! assembling the total energy in each shell and number of hits in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
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) then
fac2 = fac * (fields(i,j,k,1)**2 + fields(i,j,k,2)**2 + fields(i,j,k,3)**2)
if (akx(i).eq.0.d0) fac2 = 0.5d0 * fac2
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
return
end subroutine get_e_spec
!================================================================================-
!================================================================================-
subroutine get_gradient_statistics
use m_fields
use m_work
use x_fftw
implicit none
integer :: i
real*8 :: fac
! normalization factor
fac = one / real(nx*ny*nz_all,8)
do i = 1,3
wrk(:,:,:,i) = fields(:,:,:,i)
end do
! Taking derivatives
call x_derivative(3,'y',6)
call x_derivative(3,'x',5)
call x_derivative(2,'z',4)
call x_derivative(2,'x',3)
call x_derivative(1,'z',2)
call x_derivative(1,'y',1)
!------------------------------------------------------------
! getting vorticity and enstrophy
!------------------------------------------------------------
wrk(:,:,:,3) = wrk(:,:,:,3) - wrk(:,:,:,1) ! omega_3 = v_x - u_y
wrk(:,:,:,2) = wrk(:,:,:,2) - wrk(:,:,:,5) ! omega_2 = u_z - w_x
wrk(:,:,:,1) = wrk(:,:,:,6) - wrk(:,:,:,4) ! omega_1 = w_y - v_z
call xFFT3d(-1,1)
call xFFT3d(-1,2)
call xFFT3d(-1,3)
! getting mean enstrophy
wrk(:,:,:,0) = wrk(:,:,:,1)**2 + wrk(:,:,:,2)**2 + wrk(:,:,:,3)**2
sctmp = sum(wrk(1:nx,:,:,0)) * fac
count = 1
call MPI_REDUCE(sctmp,enstrophy,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
!------------------------------------------------------------
return
end subroutine get_gradient_statistics
!================================================================================
!================================================================================
!================================================================================
subroutine stat_scalars
use m_openmpi
use m_fields
use m_work
use x_fftw
implicit none
integer :: i, j, k, n
real*8 :: q1, q2, fac
if (.not. int_scalars) return
! getting the spectra of the scalar variances
call get_scalar_spectra
! scaling factor
fac = one / real(nx*ny*nz_all)
! --- Calculating moments of scalars
do n = 1, n_scalars
! putting the scalar in wrk0
wrk(:,:,:,0) = fields(:,:,:,3+n)
! taking derivatives
call x_derivative(0,'x',1)
call x_derivative(0,'y',2)
call x_derivative(0,'z',3)
! converting the derivatives to X-space
call xFFT3d(-1,1)
call xFFT3d(-1,2)
call xFFT3d(-1,3)
! getting the dissipation rate of the variance
wrk(:,:,:,4) = wrk(:,:,:,1)**2 + wrk(:,:,:,2)**2 + wrk(:,:,:,3)**2
q1 = two * pe(n) * sum(wrk(1:nx,:,:,4)) * fac
count = 1
call MPI_REDUCE(q1,q2,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
if (myid.eq.0) sc_diss(n) = q2
! converting the scalar itself to X-space
call xFFT3d(-1,0)
! First moment - mean
q1 = sum(wrk(1:nx,:,:,0)) * fac
count = 1
call MPI_REDUCE(q1,q2,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
if (myid.eq.0) moments(3+n,1) = q2
! Second moment - variance
q1 = sum(wrk(1:nx,:,:,0)**2) * fac
count = 1
call MPI_REDUCE(q1,q2,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
if (myid.eq.0) moments(3+n,2) = q2 - moments(3+n,1)**2
! Min and max of the scalar
q1 = minval(wrk(1:nx,:,:,0))
q2 = maxval(wrk(1:nx,:,:,0))
count = 1
call MPI_REDUCE(q1,sc_min(n),count,MPI_REAL8,MPI_MIN,0,MPI_COMM_TASK,mpi_err)
call MPI_REDUCE(q2,sc_max(n),count,MPI_REAL8,MPI_MAX,0,MPI_COMM_TASK,mpi_err)
end do
return
end subroutine stat_scalars
!================================================================================
!================================================================================
subroutine get_scalar_spectra
use m_io
use m_fields
use x_fftw
implicit none
real*8 :: sc_rad1, sc_rad2, fac, fac2
integer :: i, j, k, n, n_shell
real*8 :: energy2
! cycle over the scalars
do n = 1,n_scalars
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
! using the neergy spectra arrays to keep the scalar spectra
e_spec1 = zip
e_spec = zip
hits = 0
hits1 = 0
! assembling the total scalar energy in each shell and number of hits in each shell
do k = 1,nz
do j = 1,ny
do i = 1,nx
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) then
fac2 = fac * fields(i,j,k,3+n)**2
if (akx(i).eq.0.d0) fac2 = 0.5d0 * fac2
e_spec1(n_shell) = e_spec1(n_shell) + fac2
end if
end do
end do
end do
! reducing the energy to two arrays on master node
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! now the master node counts the energy density in each shell
master_node: if (myid.eq.0) then
! multiplying by two because we summed only half of the scalar energy
e_spec = 2.d0 * e_spec
! now the master node puts the scalar energy in the file es_sc##.gp
write(fname,"('es_',i2.2,'.gp')") n
open(900,file=fname,position='append')
write(900,"()")
write(900,"()")
write(900,"('# ITIME=',i7,' TIME=',e17.8)") ITIME, TIME
do k = 1,kmax
write(900,"(i4,4e15.6)") k,e_spec(k)
end do
close(900)
end if master_node
end do
return
end subroutine get_scalar_spectra
!================================================================================
!================================================================================
end module m_stats

75
m_timing.f90 Normal file
View file

@ -0,0 +1,75 @@
module m_timing
integer*8 :: cpu0, cpu1, cpu2, dcpu
integer*4 :: cpu_sec, cpu_min, cpu_hrs, cpu_min_total
integer*4 :: job_runlimit
!================================================================================
contains
!================================================================================
subroutine get_job_runlimit
use m_openmpi
use m_io
implicit none
logical :: there
! 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
write(out,*) 'job_runlimit (default):',job_runlimit
if(myid_world.eq.0) then
inquire(file='job_parameters.txt',exist=there)
if (there) then
open(98,file='job_parameters.txt')
read(98,*,err=5) job_runlimit
5 close(98)
end if
end if
call MPI_BCAST(job_runlimit,1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
write(out,*) 'job_runlimit: ',job_runlimit
call flush(out)
return
end subroutine get_job_runlimit
!================================================================================
!================================================================================
subroutine m_timing_init
call system_clock(cpu0,dcpu)
return
end subroutine m_timing_init
!================================================================================
!================================================================================
subroutine m_timing_check
use m_openmpi
implicit none
call system_clock(cpu1,dcpu)
cpu_sec = (cpu1-cpu0)/dcpu
cpu_min = cpu_sec/60; cpu_min_total = cpu_min
cpu_hrs = cpu_min/60
cpu_min = mod(cpu_min,60)
cpu_sec = mod(cpu_sec,60)
!!$ call MPI_BCAST(cpu_hrs,1,MPI_INTEGER4,master,MPI_COMM_TASK,mpi_err)
!!$ call MPI_BCAST(cpu_min,1,MPI_INTEGER4,master,MPI_COMM_TASK,mpi_err)
return
end subroutine m_timing_check
end module m_timing

149
m_work.f90 Normal file
View file

@ -0,0 +1,149 @@
!================================================================================
! M_WORK - module that contains working arrays wrk1....wrk15.
!
! Time-stamp: <2009-08-18 13:47:59 (chumakov)>
!================================================================================
module m_work
! use m_openmpi
use m_parameters
use m_io
implicit none
! --- work arrays
real*4,allocatable :: tmp4(:,:,:)
real*8, allocatable :: wrk(:,:,:,:)
real*8, allocatable :: wrkp(:,:)
real*8, allocatable :: rhs_old(:,:,:,:)
contains
!================================================================================
!================================================================================
!================================================================================
!================================================================================
subroutine m_work_init
use m_parameters
implicit none
!!$ write(out,*) "Inside m_work_init: ",task
!!$ call flush(out)
! allocating work arrays
if (task.eq.'hydro') then
if (les_model.ge.4) then
! The dynamic structure LES model needs more wrk arrays than usual
call m_work_allocate(max(6,3+n_scalars+n_les+4))
else
call m_work_allocate(max(6,3+n_scalars+n_les+2))
end if
elseif (task.eq.'stats') then
call m_work_allocate(6)
elseif (task.eq.'parts') then
! required recources are different for the "parts" part of the code
! we need several layers of velocities for velocity interpolation
select case (particles_tracking_scheme)
case (0)
allocate(wrk(1:nx+2,1:ny,1:3,0:0),stat=ierr)
write(out,*) "Allocated wrk(1:nx+2,1:ny,1:3,0:0)", ierr
wrk = zip
case (1)
allocate(wrk(1:nx+2,1:ny,1:3,1:3),stat=ierr)
write(out,*) "Allocated wrk(1:nx+2,1:ny,1:3,1:3)", ierr
wrk = zip
case default
stop 'wrong particles_tracking_scheme'
end select
allocate(tmp4(nx,ny,nz), stat=ierr)
write(out,*) "Allocated tmp4."
call flush(out)
else
write(out,*) 'TASK variable is set in such a way that I dont know how to allocate work arrays'
write(out,*) 'task = ',task
call my_exit(-1)
end if
!!$ write(out,*) "Finished m_work_init"
!!$ call flush(out)
return
end subroutine m_work_init
!================================================================================
!================================================================================
! allocating and defining the prescribed number of arrays
!================================================================================
!================================================================================
subroutine m_work_allocate(number)
implicit none
integer :: number
integer :: i, ierr, ierr_total
ierr = 0
!!$ write(out,"('Allocating work: ',i3)") number
!!$ call flush(out)
! array that is needed for output (nx,ny,nz)
allocate(tmp4(nx,ny,nz),stat=i); ierr = ierr + i;
! main working array, needed for FFT etc, so (nx+2,ny,nz)
allocate(wrk(nx+2,ny,nz,0:number),stat=i); ierr = ierr + i
! array for the spare RHS for Adams-Bashforth time-stepping scheme methods
if (task.eq.'hydro') then
allocate(rhs_old(nx+2,ny,nz,3+n_scalars+n_les),stat=i); ierr = ierr + i
end if
if (ierr.ne.0) then
print *,'*** WORK_INIT: error in allocation, stopping. ',ierr
print *,'*** task = ',task
print *,'*** myid = ',myid
call my_exit(-1)
stop
end if
if (allocated(wrk)) write(out,"('allocated wrk(nx+2,ny,nz,0:',i3,')')") number
if (allocated(rhs_old)) write(out,"('allocated rhs_old(nx+2,ny,nz,1:',i3,')')") 3+n_scalars+n_les
call flush(out)
tmp4 = 0.0
wrk = zip
if (allocated(rhs_old)) rhs_old = 0.d0
return
end subroutine m_work_allocate
!================================================================================
!================================================================================
subroutine m_work_exit
implicit none
! write(out,*) 'deallocaing tmp4'; call flush(out)
if (allocated(tmp4)) deallocate(tmp4)
! write(out,*) 'deallocaing wrk'; call flush(out)
if (allocated(wrk)) deallocate(wrk)
! write(out,*) 'deallocaing wrkp'; call flush(out)
if (allocated(wrkp)) deallocate(wrkp)
write(out,*) 'deallocated wrk arrays'; call flush(out)
return
end subroutine m_work_exit
end module m_work

356
main.f90 Normal file
View file

@ -0,0 +1,356 @@
program x_code
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_stats
use m_timing
use m_force
use m_rand_knuth
use m_particles
use m_filter_xfftw
use m_les
implicit none
integer :: n
character :: sym
call m_timing_init ! Setting the time zero
call m_openmpi_init
call m_io_init
call m_parameters_init
call m_les_init
call m_fields_init
call m_work_init
! allocating and initializing FFTW arrays
call x_fftw_allocate(1)
call x_fftw_init
call m_stats_init
call m_force_init
! allocating and initializing particles
if (task.eq.'parts') then
call particles_init
end if
write(out,*) "IN THE PROGRAM."
call flush(out)
! initializing the random number generator
! call rand_knuth_init
! getting the wallclock runlimit for the job
call get_job_runlimit
!-----------------------------------------------------------------------
! Starting from the beginning or from the saved flowfield
!-----------------------------------------------------------------------
if(ITMIN.eq.0) then
call begin_new
else
call begin_restart
endif
! Initializing the LES stuff
if (les) call m_les_begin
! checking divergence
if (task.eq.'hydro') call divergence
! indicators whether to use first-order in time
! for velocities and scalars
fov = .true.
fos = .true.
! need to dealias the fields at the beginning
if (task.eq.'hydro') call dealias_all
!********************************************************************************
! call benchmark
!********************************************************************************
!================================================================================
! MAIN CYCLE
!================================================================================
do 100 ITIME=ITMIN+1,ITMAX
! getting the file extension for current iteration
call get_file_ext
!--------------------------------------------------------------------------------
! now performing the core of the cycle.
! This is done with "if" rather than "select case" because if we're not
! splitting tasks then we want everything to be done consequently by the
! same set of processors.
!
! All the syncronization calls (fields_to_parts, fields_to_stats) will be
! called only if (task_split).
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
! HYDRO PART
! note that even in the case where there is no task splitting,
! 'hydro' part is still there. all processors will have task = 'hydro'
!--------------------------------------------------------------------------------
hydro: if (task.eq.'hydro') then
! ------------------------------------------------------------
! taking care of rescaling when running decaying turbulence
! if the time just was divisible by TRESCALE
! ------------------------------------------------------------
if (flow_type.eq.0 .and. floor((time-dt)/TRESCALE) .lt. floor(time/TRESCALE)) then
! ...and if we haven't rescaled NRESCALE times
if (floor(time/TRESCALE) .le. NRESCALE .and. itime.ne.1) then
write(out,*) "MAIN: Rescaling velocities"
call flush(out)
call velocity_rescale
! after rescaling, the time-sceping needs to be first order
fov = .true.; fos = .true.
if (.not. task_split .and. mod(itime,iprint1).eq.0) call stat_main
end if
end if
! RHS for scalars
call rhs_scalars
! now the velocities in x-space are contained in wrk1...3
! if we are moving particles, then we want to send the velocity field
! to the "parts" part of the code
if (task_split) call fields_to_parts
! advance scalars - either Euler or Adams-Bashforth
if (int_scalars .or. n_les > 0) then
call flush(out)
n = 3 + n_scalars + n_les
if (fos) then
rhs_old(:,:,:,4:n) = wrk(:,:,:,4:n)
fields(:,:,:,4:n) = fields(:,:,:,4:n) + dt * rhs_old(:,:,:,4:n)
fos = .false.
else
fields(:,:,:,4:n) = fields(:,:,:,4:n) + &
dt * ( 1.5d0 * wrk(:,:,:,4:n) - 0.5d0 * rhs_old(:,:,:,4:n) )
rhs_old(:,:,:,4:n) = wrk(:,:,:,4:n)
end if
end if
! RHS for velocities
call rhs_velocity
! adding forcing, if computing a forced flow
if (flow_type.eq.1) call force_velocity
! advance velocity - either Euler or Adams-Bashforth
if (fov) then
rhs_old(:,:,:,1:3) = wrk(:,:,:,1:3)
fields(:,:,:,1:3) = fields(:,:,:,1:3) + dt * rhs_old(:,:,:,1:3)
fov = .false.
else
fields(:,:,:,1:3) = fields(:,:,:,1:3) + &
dt * ( 1.5d0 * wrk(:,:,:,1:3) - 0.5d0 * rhs_old(:,:,:,1:3) )
rhs_old(:,:,:,1:3) = wrk(:,:,:,1:3)
end if
! solve for pressure and update velocities so they are incompressible
call pressure
! advance the time
TIME = TIME + DT
! write the restart file if it's the time
if (mod(itime,IPRINT2).eq.0) call restart_write_parallel
! change the timestep in case we're running with variable timestep
if (variable_dt) call my_dt
! CPU usage statistics
if (mod(itime,iprint1).eq.0) then
call m_timing_check
if (mod(itime,iwrite4).eq.0) then
sym = "*"
else
sym = " "
end if
write(out,9000) itime,time,dt,courant,cpu_hrs,cpu_min,cpu_sec,&
sym,les_model_name
call flush(out)
end if
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
! send the velocities to the "stats" part of the code for statistics
if (task_split) call fields_to_stats
! checking if we need to stop the calculations due to simulation time
if (TIME.gt.TMAX) call my_exit(1)
! checking if we need to start advancing scalars
if (n_scalars.gt.0 .and. .not.int_scalars .and. time.gt.TSCALAR) then
int_scalars = .true.
call init_scalars
write(out,*) "Starting to move the scalars."
call flush(out)
end if
end if
end if hydro
!--------------------------------------------------------------------------------
! STATISTICS PART
!--------------------------------------------------------------------------------
stats: if (task.eq.'stats' .or. .not.task_split) then
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
! if this is a separate set of processors, then...
stats_task_split: if (task_split) then
! checking if we need to stop the calculations due to simulation time
if (TIME.gt.TMAX) call my_exit(1)
end if stats_task_split
! these are executed regardless of the processor configuration
if (task_split) call fields_to_stats
if (mod(itime,iprint1).eq.0) call stat_main
if (mod(itime,iwrite4).eq.0) call io_write_4
end if
end if stats
!--------------------------------------------------------------------------------
! PARTICLE PARTS
! NOTE: This is not enabled to work when not task_split.
! Need to return to it later.
! Currently the particles can be calculated only if we split the tasks due to
! requirements on the wrk array sizes in the particle interpolation routines.
!--------------------------------------------------------------------------------
particles: if (task.eq.'parts') then
call fields_to_parts
if (int_particles) then
call particles_move
if (mod(itime,iwrite4).eq.0) call particles_restart_write_binary
end if
if (mod(itime,iprint1).eq.0 .or. mod(itime,iwrite4).eq.0) then
if (TIME.gt.TMAX) call my_exit(1)
end if
end if particles
!!$!--------------------------------------------------------------------------------
!!$! OTHER PARTS
!!$!--------------------------------------------------------------------------------
!!$
!!$
!!$ write(out,*) "skipping the time step",ITIME
!!$ call flush(out)
!!$
!!$
!--------------------------------------------------------------------------------
! COMMON PARTS
!--------------------------------------------------------------------------------
! every 10 iterations checking
! 1) for the run time: are we getting close to the job_runlimit?
! 2) for the user termination: is there a file "stop" in directory?
if (mod(ITIME,10).eq.0) then
! synchronize all processors, hard
!!$ call MPI_BARRIER(MPI_COMM_WORLD,mpi_err)
if (myid_world.eq.0) call m_timing_check
count = 1
call MPI_BCAST(cpu_min_total,count,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
! allowing 5 extra minutes for writing restart file
! note that for large-scale calculations (e.g. 1024^3)
! the restart writing time can be long (up to 20 minutes or so).
! this should be taken care of in the job submission script
! via file job_parameters.txt
if (cpu_min_total+5 .gt. job_runlimit) call my_exit(2)
! user termination. If the file "stop" is in the directory, stop
inquire(file='stop',exist=there)
if (there) call my_exit(3)
end if
100 continue
!================================================================================
!--------------------------------------------------------------------------------
! In a case when we've gone to ITMAX, write the restart file
!--------------------------------------------------------------------------------
ITIME = ITIME-1
if (task.eq.'hydro') call restart_write_parallel
call my_exit(0)
call m_openmpi_exit
stop
9000 format('ITIME=',i6,3x,'TIME=',f8.4,4x,'DT=',f8.5,3x,'Courn= ',f6.4, &
2x,'CPU:(',i4.4,':',i2.2,':',i2.2,')',x,a1,x,a3)
end program x_code
!=============================================================================
subroutine benchmark
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_stats
use m_timing
use m_force
use m_rand_knuth
use m_particles
use m_filter_xfftw
use m_les
implicit none
integer :: n, nmax
call m_timing_init
benchmarking = .true.
bm = 0
wrk(:,:,:,1) = fields(:,:,:,1)
do n = 1, nmax
call xfft3d(1,1)
call xfft3d(-1,1)
end do
if (myid.eq.0) then
write(out,*) "BENF: statistics on forward transform"
write(out,*) "BENF: R2C: ", bm(1)/nmax
write(out,*) "BENF: T13: ", bm(2)/nmax
write(out,*) "BENF: C2C: ", bm(3)/nmax
write(out,*) "BENF: ====="
write(out,*) "BENF: TOT: ", bm(11)/nmax
write(out,*) "BENB: statistics on backward transform"
write(out,*) "BENB: C2C: ", bm(4)/nmax
write(out,*) "BENB: T13: ", bm(5)/nmax
write(out,*) "BENB: C2R: ", bm(6)/nmax
write(out,*) "BENB: ====="
write(out,*) "BENB: TOT: ", bm(12)/nmax
end if
close(out)
stop
end subroutine benchmark

34
my_dt.f90 Normal file
View file

@ -0,0 +1,34 @@
subroutine my_dt
use m_openmpi
use m_parameters
implicit none
integer :: n
real*8 :: courant_max=0.1
if (flow_type.le.1) then
! Courant number is calculated every iteration in rhsv.f
if (courant.lt.courant_max) then
DT = DT * 1.01d0
DT = min(DT,0.01)
else
DT = DT/1.05d0
end if
!----------------------------------------------------------------------
! make sure DT is appropriate for scalars
! DT < 0.09 * dx^2*Pe
!----------------------------------------------------------------------
if (int_scalars) then
do n=1,n_scalars
!!$ DT = min(DT,0.09*dx**2/PE(n))
DT = min(DT,DT*courant_max/courant/sc(n))
end do
end if
end if
return
end subroutine my_dt

123
my_exit.f90 Normal file
View file

@ -0,0 +1,123 @@
!================================================================================
! MY_EXT - de-initializes everything, called at the end the main program
! or whenever we need to stop the program anywhere
!
! reason - the passed variable that describes the particular reason why
! we want to stop the program.
!================================================================================
subroutine my_exit(reason)
use m_openmpi
use m_io
use m_fields
use m_work
! use x_fftw
use m_particles
implicit none
integer :: reason
write(out,"('my_exit with reason ',i4)") reason
select case (reason)
case(0)
write(out,*) '---------------------------------------------'
write(out,*) ' NORMAL TERMIATION'
write(out,*) '---------------------------------------------'
case(1)
write(out,*) '---------------------------------------------'
write(out,*) ' TIME TERMINATION'
write(out,*) '---------------------------------------------'
call flush(out)
case(2)
write(out,*) '---------------------------------------------'
write(out,*) ' RUN-TIME TERMINATION'
write(out,*) '---------------------------------------------'
case(3)
write(out,*) '---------------------------------------------'
write(out,*) ' USER TERMINATION'
write(out,*) '---------------------------------------------'
case default
write(out,*) '---------------------------------------------'
write(out,*) ' TERMINATION FOR NO APPARENT REASON'
write(out,*) '---------------------------------------------'
end select
if (reason.ge.0) then
if (task.eq.'hydro') call restart_write_parallel
if (task.eq.'parts') call particles_restart_write_binary
end if
write(out,*) "Done."
call flush(out)
close(out)
stop
! call m_fields_exit
! call m_work_exit
! call x_fftw_allocate(-1)
! call m_io_exit
! call m_openmpi_exit
return
end subroutine my_exit
!!$!================================================================================
!!$! MY_INIT - initializes everything, called at the beginning of the main program
!!$!================================================================================
!!$subroutine my_init
!!$
!!$! --- modules used
!!$ use m_openmpi
!!$ use m_io
!!$ use m_parameters
!!$ use m_fields
!!$ use m_work
!!$ use x_fftw
!!$ use m_filter_xfftw
!!$ implicit none
!!$ integer :: iargc
!!$
!!$! --- initializing
!!$ call openmpi_init
!!$ call io_init
!!$ call parameters_init
!!$ call fields_init
!!$ call work_init(15)
!!$ call x_fftw_allocate(1)
!!$ call x_fftw_init
!!$
!!$! --- getting the filter size
!!$ if(iargc().eq.0) then
!!$ call getarg(0,tmp_str)
!!$ write(out,*) 'Format: ',trim(tmp_str),' <filter_size>'
!!$ write(*,*) 'Format: ',trim(tmp_str),' <filter_size>'
!!$ call my_exit(-1)
!!$ end if
!!$
!!$ ftype = 2
!!$ call getarg(1,txt7)
!!$ read(txt7,*) filter_size
!!$ write(out,*) 'Filter size = ',filter_size
!!$ call flush(out)
!!$
!!$ call filter_xfftw_init
!!$
!!$return
!!$end subroutine my_init

96
pressure.f90 Normal file
View file

@ -0,0 +1,96 @@
subroutine pressure
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: i, j, k, n
real*8 :: div1, div2, lapl1, lapl2, p1, p2
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! getting the divergence (i*k*\hat(u))
! remember that in the Fourier space indicies go as (ix,iz,iy)
! getting divergence
div2 = akx(i )*fields(i ,j,k,1) + aky(k)*fields(i ,j,k,2) + akz(j)*fields(i ,j,k,3)
div1 = - ( akx(i+1)*fields(i+1,j,k,1) + aky(k)*fields(i+1,j,k,2) + akz(j)*fields(i+1,j,k,3) )
! inverce laplace operator
lapl1 = akx(i )**2 + aky(k)**2 + akz(j)**2
lapl2 = akx(i+1)**2 + aky(k)**2 + akz(j)**2
if (lapl1.eq.0.d0) lapl1 = 9e20
if (lapl2.eq.0.d0) lapl2 = 9e20
! calculating pressure
p1 = - div1 / lapl1
p2 = - div2 / lapl2
! Taking derivatives of the pressure and subtracting from the corresponding velocities
fields(i ,j,k,1) = fields(i ,j,k,1) + p2 * akx(i+1)
fields(i+1,j,k,1) = fields(i+1,j,k,1) - p1 * akx(i )
fields(i ,j,k,2) = fields(i ,j,k,2) + p2 * aky(k)
fields(i+1,j,k,2) = fields(i+1,j,k,2) - p1 * aky(k)
fields(i ,j,k,3) = fields(i ,j,k,3) + p2 * akz(j)
fields(i+1,j,k,3) = fields(i+1,j,k,3) - p1 * akz(j)
end do
end do
end do
return
end subroutine pressure
subroutine divergence
use m_openmpi
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: i,j,k
real*8 :: dmin, dmax, d1
wrk(:,:,:,1:3) = fields(:,:,:,1:3)
call x_derivative(1,'x',4)
call x_derivative(2,'y',5)
call x_derivative(3,'z',6)
call xFFT3d(-1,4)
call xFFT3d(-1,5)
call xFFT3d(-1,6)
wrk(:,:,:,0) = wrk(:,:,:,4) + wrk(:,:,:,5) + wrk(:,:,:,6)
d1 = minval(wrk(1:nx,:,:,0))
call MPI_REDUCE(d1,dmin,1,MPI_REAL8,MPI_MIN,0,MPI_COMM_TASK,mpi_err)
d1 = maxval(wrk(1:nx,:,:,0))
call MPI_REDUCE(d1,dmax,1,MPI_REAL8,MPI_MAX,0,MPI_COMM_TASK,mpi_err)
if (myid.eq.0) then
write(out,*) 'divergence:',dmin,dmax
!!$ print *, 'divergence:',dmin,dmax
call flush(out)
end if
return
end subroutine divergence

506
restart_io.f90 Normal file
View file

@ -0,0 +1,506 @@
!!$!================================================================================
!!$!================================================================================
!!$!================================================================================
!!$subroutine restart_read
!!$
!!$ use m_openmpi
!!$ use m_parameters
!!$ use m_io
!!$ use m_fields
!!$ use m_work
!!$ use x_fftw
!!$ use m_particles
!!$ implicit none
!!$
!!$ integer*4 :: nx1,ny1,nz1, nums1, MST1, nums_read
!!$ integer :: i, j, k, n
!!$
!!$ real*8 :: ST
!!$
!!$ fname = run_name//'.64.'//file_ext
!!$ inquire(file=fname,exist=there)
!!$ if(.not.there) then
!!$ write(out,*) '*** error: Cannot find file : '//trim(fname)
!!$ stop
!!$ end if
!!$
!!$ write(out,*) 'Reading from the file (seq): ',trim(fname)
!!$ call flush(out)
!!$
!!$ ! root process reads parameters from the file
!!$ if (myid_world.eq.0) then
!!$ open(91,file=fname,form='unformatted',access='stream')
!!$ read(91) nx1, ny1, nz1, nums1, MST1, TIME, DT
!!$ end if
!!$
!!$ call MPI_BCAST(nx1, 1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
!!$ call MPI_BCAST(ny1, 1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
!!$ call MPI_BCAST(nz1, 1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
!!$ call MPI_BCAST(nums1,1,MPI_INTEGER4,0,MPI_COMM_WORLD,mpi_err)
!!$
!!$ call MPI_BCAST(TIME,1,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
!!$ call MPI_BCAST( DT,1,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
!!$
!!$ ! everyone checks of the parameters coinside with what's in the .in file
!!$ if (nx.ne.nx1 .or. ny.ne.ny1 .or. nz_all.ne.nz1) then
!!$ write(out,*) '*** error: Dimensions are different'
!!$ write(out,*) '*** .in file: ',nx,ny,nz_all
!!$ write(out,*) '*** restart file: ',nx1,ny1,nz1
!!$ call flush(out)
!!$ stop
!!$ end if
!!$
!!$!-----------------------------------------------------------------------
!!$! dealing with scalars.
!!$!
!!$! The number of scalars can be varied throughout the simulation.
!!$! If the restart file has fewer scalars than the
!!$! .in file, the scalars are added and initialized according to
!!$! their description in the .in-file.
!!$! If the restart file has more scalars than the .in file, the
!!$! extra scalars are dropped.
!!$!
!!$! in short, whatever is specfied in the .in file, prevails.
!!$!-----------------------------------------------------------------------
!!$
!!$ if (n_scalars.lt.nums1) then
!!$
!!$ write(out,*) ' WARNING: nums in restart file:',nums1
!!$ write(out,*) ' nums in .in file :',n_scalars
!!$ write(out,'(''Losing '',i3,'' scalars.'')') nums1-n_scalars
!!$ call flush(out)
!!$ nums_read = n_scalars
!!$
!!$ else if (n_scalars.gt.nums1) then
!!$
!!$ write(out,*) ' WARNING: nums in restart file:',nums1
!!$ write(out,*) ' nums in .in file :',n_scalars
!!$ write(out,'(''Adding '',i3,'' scalars.'')') n_scalars-nums1
!!$ call flush(out)
!!$ nums_read = nums1
!!$
!!$
!!$ else
!!$
!!$ nums_read = n_scalars
!!$
!!$ end if
!!$
!!$!----------------------------------------------------------------------
!!$! ------------ reading the stuff from the restart file ---------------
!!$!----------------------------------------------------------------------
!!$
!!$ ! only the hydro part of processors is involved in this
!!$ hydro_only: if (task.eq.'hydro') then
!!$ count = (nx+2) * ny * nz
!!$ ! the root reads everything and sends to the slaves
!!$ if (myid.eq.0) then
!!$ do n = 1, 3 + nums_read + n_les
!!$ ! first chunk belongs to the root
!!$ read(91) (((fields(i,j,k,n),i=1,nx),j=1,ny),k=1,nz)
!!$ ! the rest gets read and sent to the appropriate porcess
!!$ do id_to = 1,numprocs-1
!!$ read(91) (((wrk(i,j,k,1),i=1,nx),j=1,ny),k=1,nz)
!!$ tag = (3+nums_read) * id_to + n-1
!!$ call MPI_SEND(wrk(1,1,1,1),count,MPI_REAL8,id_to,tag,MPI_COMM_TASK,mpi_err)
!!$
!!$
!!$ end do
!!$ end do
!!$ ! then the root closes the restart file
!!$ close(91)
!!$ else
!!$ ! the slaves receive and put it into ss array
!!$ do n = 1, 3 + nums_read + n_les
!!$ tag = (3+nums_read) * myid + n-1
!!$ call MPI_RECV(fields(1,1,1,n),count,MPI_REAL8,0,tag,MPI_COMM_TASK,mpi_status,mpi_err)
!!$ end do
!!$ end if
!!$
!!$ end if hydro_only
!!$
!!$ return
!!$end subroutine restart_read
!!$
!!$
!!$
!!$
!!$!================================================================================
!!$!================================================================================
!!$!================================================================================
!!$subroutine restart_write
!!$
!!$! This routine is called to GENERATE RESTART FILES
!!$! The file is written without MPI-2 tricks
!!$! the stuff gets sent to the root process adn written out in
!!$! some orderly fashion
!!$
!!$ use m_openmpi
!!$ use m_parameters
!!$ use m_io
!!$ use m_fields
!!$ use m_work
!!$ use x_fftw
!!$ implicit none
!!$
!!$ integer :: n, nums_out, i, j, k
!!$
!!$ real*8 :: ST
!!$ integer :: MST
!!$
!!$ if (itime.eq.last_dump) return
!!$
!!$
!!$!---------------------------------------------------------------------
!!$! dumping the restart file with particles
!!$!---------------------------------------------------------------------
!!$! if (int_particles) call particles_restart_write
!!$! if (int_particles) call particles_restart_write_binary
!!$
!!$ ! how many scalars to write
!!$ nums_out = 0
!!$ if (int_scalars) nums_out = n_scalars
!!$
!!$ ! first FFT everything to real space
!!$ wrk(:,:,:,1:3+nums_out+n_les) = fields(:,:,:,1:3+nums_out+n_les)
!!$
!!$ ! --------------- writing process ------------------
!!$
!!$ fname = run_name//'.64.'//file_ext
!!$
!!$ ! if not the root, just send the stuff to the root
!!$ if (myid.ne.0) then
!!$
!!$ do n = 1 , 3 + nums_out + n_les
!!$
!!$ wrk(:,:,:,0) = wrk(:,:,:,n)
!!$ tag = (3+nums_out+n_les) * myid + n-1
!!$ count = (nx+2) * ny * nz
!!$ call MPI_ISEND(wrk(1,1,1,0),count,MPI_REAL8,0,tag,MPI_COMM_TASK,request,mpi_err)
!!$
!!$! write(out,'(''sending var '',i3,'' to '',i4,'': '',i3)') n,0,mpi_err
!!$! call flush(out)
!!$
!!$ call MPI_WAIT(request,mpi_status,mpi_err)
!!$
!!$! write(out,'(''sent var '',i3,'' to '',i4,'': '',i3)') n,0,mpi_err
!!$! call flush(out)
!!$
!!$ end do
!!$
!!$ else
!!$ ! if it's the root, then write the restart file
!!$!! open(91,file=fname,form='binary')
!!$ open(91,file=fname,form='unformatted', access='stream')
!!$
!!$ ! first write the parameters
!!$ ST = zip
!!$ write(91) int(nx,4),int(ny,4),int(nz*numprocs,4),int(nums_out,4),int(MST,4),TIME,DT
!!$
!!$ ! then write the variables, one by one
!!$ do n = 1 , 3 + nums_out + n_les
!!$
!!$ ! first write the chunk from the root process
!!$ wrk(:,:,:,0) = wrk(:,:,:,n)
!!$ write(91) (((wrk(i,j,k,0),i=1,nx),j=1,ny),k=1,nz)
!!$ ! then receive chinks of the same variable from each process and write it out
!!$ do id_from=1,numprocs-1
!!$ tag = (3+nums_out+n_les) * id_from + n-1
!!$ count = (nx+2) * ny * nz
!!$ call MPI_RECV(wrk(1,1,1,0),count,MPI_REAL8,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err)
!!$
!!$! write(out,'(''received var '',i3,'' from '',i4,'': '',i3)') n,id_from,mpi_err
!!$! call flush(out)
!!$
!!$ write(91) (((wrk(i,j,k,0),i=1,nx),j=1,ny),k=1,nz)
!!$ end do
!!$ end do
!!$ close(91)
!!$
!!$ end if
!!$
!!$ write(out,*) '------------------------------------------------'
!!$ write(out,*) 'Restart file written (seq): '//trim(fname)
!!$ write(out,"(' Velocities and ',i3,' scalars (incl. LES)')") nums_out+n_les
!!$ write(out,"(' Restart file time = ',f15.10,i7)") time,itime
!!$ write(out,*) '------------------------------------------------'
!!$ call flush(out)
!!$
!!$
!!$ ! setting the variable last_dump to current timestep number
!!$ last_dump = ITIME
!!$
!!$ return
!!$end subroutine restart_write
!!$
!!$
!================================================================================
!================================================================================
!================================================================================
subroutine restart_write_parallel
! This routine is called to GENERATE RESTART FILES
! The file is written using the collective write (MPI-2 standard)
use m_openmpi
use m_parameters
use m_io
use m_fields
use m_work
use x_fftw
implicit none
integer :: n, nums_out, i, j, k
integer :: MST
integer(kind=MPI_INTEGER_KIND) :: fh
integer(kind=MPI_OFFSET_KIND) :: offset
real*8, allocatable :: sctmp8(:,:,:)
integer*4 :: nx1, ny1, nz1, nums1, nles1
real*8 :: ST
if (itime.eq.last_dump) return
! how many scalars to write
nums_out = 0
if (int_scalars) nums_out = n_scalars
! using wrk array
wrk(:,:,:,1:3+nums_out) = fields(:,:,:,1:3+nums_out)
if (n_les>0) wrk(:,:,:,3+nums_out+1:3+nums_out+n_les) = fields(:,:,:,3+n_scalars+1:3+n_scalars+n_les)
! --------------- writing process ------------------
fname = run_name//'.64.'//file_ext
! allocating the temporary array sctmp8
allocate(sctmp8(nx,ny,nz),stat=ierr)
if (ierr.ne.0) stop '*** RESTART_READ_PARALLEL: cannot allocate sctmp8'
sctmp8 = zip
! 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
nx1 = nx; ny1 = ny; nz1 = nz_all; nles1 = n_les; nums1 = nums_out;
count = 1
ST = zip
call MPI_FILE_WRITE(fh, nx1, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, ny1, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, nz1, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, nums1, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, nles1, count, MPI_INTEGER4, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, TIME, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_WRITE(fh, DT, count, MPI_REAL8, mpi_status, mpi_err)
end if
! all nodes write their stuff into the file
! first writing the velocities and passive scalars
writing_fields: do n = 1, 3 + nums_out + n_les
offset = 36 + (n-1)*nx*ny*nz_all*8 + myid*nx*ny*nz*8
count = nx * ny * nz
! note that we want the data from the restart to have dimensions (nx,ny,nz),
! while the fields array has fimensions (nx+2,ny,nz).
! this is an artefact of the times when the code used to write the variables in real space.
! that is why we need to duplicate each field in the sctmp array first, and then
! write sctmp8 into the file with appropriate offset
sctmp8(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,n)
call MPI_FILE_WRITE_AT(fh, offset, sctmp8, count, MPI_REAL8, mpi_status, mpi_err)
end do writing_fields
call MPI_FILE_CLOSE(fh, mpi_err)
call MPI_INFO_FREE(mpi_info, mpi_err)
deallocate(sctmp8)
write(out,*) '------------------------------------------------'
write(out,*) 'Restart file written (par): '//trim(fname)
write(out,"(' Velocities and ',i3,' passive scalars')") nums_out
if (n_les>0) write(out,"(' Also wrote',i3,' LES scalars)')") n_les
write(out,"(' Restart file time = ',f15.10,i7)") time, itime
write(out,*) '------------------------------------------------'
call flush(out)
last_dump = ITIME
return
end subroutine restart_write_parallel
!================================================================================
!================================================================================
!================================================================================
!================================================================================
subroutine restart_read_parallel
use m_openmpi
use m_parameters
use m_io
use m_fields
use m_work
use x_fftw
use m_particles
implicit none
integer*4 :: nx1,ny1,nz1, nums1, nles1, nums_read
integer :: i, j, k, n, n_skip
integer(kind=MPI_INTEGER_KIND) :: fh
integer(kind=MPI_OFFSET_KIND) :: offset
real*8, allocatable :: sctmp8(:,:,:)
! checking if the restart file exists
fname = run_name//'.64.'//file_ext
inquire(file=fname,exist=there)
if(.not.there) then
write(out,*) '*** error: Cannot find file : '//trim(fname)
stop
end if
write(out,*) 'Reading from the file (par): ',trim(fname)
call flush(out)
! ----------------------------------------------------------------------
! first reading the parameters from the restart file.
! the root process opens it and reads the parameters, then broadcasts
! the parameters. After that it's decided if the parameters make sense,
! how many scalars to read etc.
! ----------------------------------------------------------------------
if (myid.eq.0) then
open(91,file=fname,form='unformatted',access='stream')
read(91) nx1, ny1, nz1, nums1, nles1, TIME, DT
close(91)
end if
call MPI_BCAST(nx1, 1,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST(ny1, 1,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST(nz1, 1,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST(nums1,1,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST(nles1,1,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST(TIME,1,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
call MPI_BCAST( DT,1,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
! checking if the array sizes are the same in .in file and restart file
if (nx.ne.nx1 .or. ny.ne.ny1 .or. nz_all.ne.nz1) then
write(out,*) '*** error: Dimensions are different'
write(out,*) '*** .in file: ',nx,ny,nz_all
write(out,*) '*** restart file: ',nx1,ny1,nz1
call flush(out)
call my_exit(-1)
end if
! checking if the number of LES quantities is the same in .in and restart files
if (n_les .ne. nles1) then
write(out,*) '*** WARNING : Different values of n_les:'
write(out,*) '*** .in file: ',n_les
write(out,*) '*** restart file: ',nles1
write(out,*) '*** Make sure you are running the same simulation.'
call flush(out)
end if
!-----------------------------------------------------------------------
! dealing with scalars.
!
! The number of scalars can be varied throughout the simulation.
! If the restart file has fewer scalars than the
! .in file, the scalars are added and initialized according to
! their description in the .in-file.
! If the restart file has more scalars than the .in file, the
! extra scalars are dropped.
!
! in short, whatever is specfied in the .in file, prevails.
!-----------------------------------------------------------------------
if (n_scalars.lt.nums1) then
write(out,*) ' WARNING: nums in restart file:',nums1
write(out,*) ' nums in .in file :',n_scalars
write(out,'(''Losing '',i3,'' scalars.'')') nums1-n_scalars
call flush(out)
nums_read = n_scalars
else if (n_scalars.gt.nums1) then
write(out,*) ' WARNING: nums in restart file:',nums1
write(out,*) ' nums in .in file :',n_scalars
write(out,'(''Adding '',i3,'' scalars.'')') n_scalars-nums1
call flush(out)
nums_read = nums1
! initializing the added scalars
do n=nums1+1,n_scalars
call init_scalar(n)
end do
else
nums_read = n_scalars
end if
!----------------------------------------------------------------------
! ------------ reading the stuff from the restart file ---------------
!----------------------------------------------------------------------
! allocating the temporary array sctmp8
allocate(sctmp8(nx,ny,nz),stat=ierr)
if (ierr.ne.0) stop '*** RESTART_READ_PARALLEL: cannot allocate sctmp8'
sctmp8 = zip
! 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)
! note that the data from the restart file has dimensions (nx,ny,nz),
! while the fields array has fimensions (nx+2,ny,nz).
! that is why we need to read each field in the sctmp array first, and then
! rearrange it and put into the fields array.
reading_fields: do n = 1, 3 + nums_read
write(out,"('Reading variable # ',i3)") n
call flush(out)
offset = 36 + (n-1)*nx*ny*nz_all*8 + myid*nx*ny*nz*8
count = nx * ny * nz
! call MPI_FILE_READ_AT(fh, offset, sctmp8, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_READ_AT_ALL(fh, offset, sctmp8, count, MPI_REAL8, mpi_status, mpi_err)
fields(1:nx,1:ny,1:nz,n) = sctmp8(1:nx,1:ny,1:nz)
end do reading_fields
! now reading the LES variables. They are stored after the fields
! u,v,w,sc(1...nums1), so we are applying the offset based on the
! number of scalars in the restart file.
reading_les_fields: do n = 1, n_les
write(out,"('Reading LES variable # ',i3)") n
call flush(out)
n_skip = 3 + nums1
offset = 36 + (n_skip+n-1)*nx*ny*nz_all*8 + myid*nx*ny*nz*8
count = nx * ny * nz
! call MPI_FILE_READ_AT(fh, offset, sctmp8, count, MPI_REAL8, mpi_status, mpi_err)
call MPI_FILE_READ_AT_ALL(fh, offset, sctmp8, count, MPI_REAL8, mpi_status, mpi_err)
fields(1:nx,1:ny,1:nz,3+n_scalars+n) = sctmp8(1:nx,1:ny,1:nz)
end do reading_les_fields
call MPI_FILE_CLOSE(fh, mpi_err)
call MPI_INFO_FREE(mpi_info, mpi_err)
deallocate(sctmp8)
write(out,*) "Restart file successfully read."
call flush(out)
return
end subroutine restart_read_parallel

594
rhs_scalars.f90 Normal file
View file

@ -0,0 +1,594 @@
subroutine rhs_scalars
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_timing
use m_les
implicit none
integer :: i, j, k, n, n1, n2, nv, ns_lo, ns_hi
real*8 :: rtmp1, rtmp2, wnum2, r11, r12, r21, r22, r31, r32
! calculate turbulent viscosity, if any
if (les) call les_get_turb_visc
! If we're not advancing scalars, put the real-space velocities
! in wrk1...3 and return
! same is done if dealias=0, that is, we use 2/3 rule.
! if dealias=1 (phase shifts) then this is done later in the subroutine
! also if doing LES and n_les (the # of les-related auxilary scalars) is
! greater than 0, then we need to transport these LES-related scalars,
! even if n_scalars=0.
! thus we do the obligatory part (tranfer velocities to X-space) and quit
! only when we are not transporting any scalars and if there are no
! LES-related scalars.
if (.not.int_scalars .and. n_les .eq. 0) then
! converting velocities to the real space and returning
wrk(:,:,:,1:3) = fields(:,:,:,1:3)
do n = 1,3
call xFFT3d(-1,n)
end do
return
end if
! making the RHS for all scalars zero
wrk(:,:,:,4:3+n_scalars+n_les) = zip
!--------------------------------------------------------------------------------
! If dealias=0, performing the 2/3 rule dealiasing on scalars
!--------------------------------------------------------------------------------
if (dealias.eq.0) then
! converting velocities to the real space
wrk(:,:,:,1:3) = fields(:,:,:,1:3)
do n = 1,3
call xFFT3d(-1,n)
end do
! Do each scalar one at a time. Keep the velocities in wrk1:3 intact
! because they are needed later.
! first we need to know which scalars do we want to transport.
! There are three cases:
! (0) No LES extra scalars, and passive scalars have not been initialized yet
! Thus this subroutine calculates the IFFT of velocities and exits.
! This is taken care of earlier.
! (1) Both passive scalars and LES extra scalars are transported
! This is possible when n_les > 0 and int_scalars=.true.
! (2) Only LES extra scalars are transported
! This is possible if and only if (.not.int_scalars .and. n_les>0)
!
! The last two cases are taken care of by prescribing ns_lo and ns_hi, the
! smallest and largest number of the scalar that needs to be transported.
ns_lo = 1;
ns_hi = n_scalars + n_les
if (.not.int_scalars) ns_lo = n_scalars + 1
do n = ns_lo, ns_hi
wrk(:,:,:,0) = fields(:,:,:,3+n)
call xFFT3d(-1,0)
! Products of the scalar and velocities
do i = 1,3
wrk(:,:,:,n+2+i) = wrk(:,:,:,0) * wrk(:,:,:,i)
call xFFT3d(1,n+2+i)
end do
! Assembling the RHS in wrk(:,:,:,3+n)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! If the dealiasing option is 2/3-rule (dealias=0) then we retain the modes
! inside the cube described by $| k_i | \leq k_{max}$, $i=1,2,3$.
! The rest of the modes is purged
if (ialias(i,j,k) .gt. 0) then
! all the wavenumbers that are greater than kmax get zeroed out
wrk(i ,j,k,3+n) = zip
wrk(i+1,j,k,3+n) = zip
else
! taking the convective term, multiply it by "i"
! (see how it's done in x_fftw.f90)
! and adding the diffusion term
! also using the fact that the waveunmbers for (i,j,k) are the same
! as wavenumbers for (i+1,j,k)
! i * (a + ib) + d = -b + ia + d
rtmp1 = akx(i+1)*wrk(i+1,j,k,3+n) + aky(k)*wrk(i+1,j,k,4+n) + akz(j)*wrk(i+1,j,k,5+n)
rtmp2 = akx(i )*wrk(i ,j,k,3+n) + aky(k)*wrk(i ,j,k,4+n) + akz(j)*wrk(i ,j,k,5+n)
wnum2 = akx(i)**2 + aky(k)**2 + akz(j)**2
wrk(i ,j,k,3+n) = rtmp1 - pe(n) * wnum2*fields(i ,j,k,3+n)
wrk(i+1,j,k,3+n) = - rtmp2 - pe(n) * wnum2*fields(i+1,j,k,3+n)
end if
end do
end do
end do
! Now adding the reaction part (for scalars only, not for LES-related quantities
! that are formally scalars with indicies n_scalars+1...n_scalars+n_les )
if (n .le. n_scalars) then
if (scalar_type(n).ge.100) then
call add_reaction(n)
call dealias_rhs(3+n)
end if
end if
end do
end if
!--------------------------------------------------------------------------------
! If dealias=1, performing the phase shift and truncation on scalars
! The main ideology is as follows. First we evaluate the phase-shifted
! quantities, then not phase shifted. This is done in order to have more
! working space: at the beginning we have all work arrays available, but at the
! end we must put sin/cos factors in wrk0 and u.v.w in real space in wrk1...3.
! They are going to be used again in rhs_velocity later.
!--------------------------------------------------------------------------------
phase_shifting_dealiasing: if (dealias.eq.1) then
! define the sin/cos factors that are used in phase shifting.
! computing sines and cosines for the phase shift of dx/2,dy/2,dz/2
! and putting them into wrk0
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
wrk(i ,j,k,0) = cos(half*(akx(i )+aky(k)+akz(j))*dx)
wrk(i+1,j,k,0) = sin(half*(akx(i+1)+aky(k)+akz(j))*dx)
end do
end do
end do
! Doing all scalars at the same time. We can do this because the
! number of work arrays that are available to us is 3+n+2. First
! three later will be taken by velocities, and the last two are
! primary work arrays here.
! First, get phase shifted velocities and scalars
do n = 1, 3 + n_scalars + n_les
! phase-shifting the quantity
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
wrk(i ,j,k,n) = fields(i ,j,k,n) * wrk(i,j,k,0) - fields(i+1,j,k,n) * wrk(i+1,j,k,0)
wrk(i+1,j,k,n) = fields(i+1,j,k,n) * wrk(i,j,k,0) + fields(i ,j,k,n) * wrk(i+1,j,k,0)
end do
end do
end do
! transforming it to real space
call xFFT3d(-1,n)
end do
! now we have two vacant arrays: n_scalars+n_les+4 and n_scalars+n_les+5. Work in them
n1 = n_scalars + n_les + 4
n2 = n_scalars + n_les + 5
! do one scalar at a time
phase_shifted_rhs: do n = 4, 3 + n_scalars + n_les
! getting all three products of phase-shifted scalar and phase-shifted velocities
! using three work arrays: n, n1 and n2
wrk(:,:,:,n1) = wrk(:,:,:,n) * wrk(:,:,:,1)
wrk(:,:,:,n2) = wrk(:,:,:,n) * wrk(:,:,:,2)
wrk(:,:,:,n ) = wrk(:,:,:,n) * wrk(:,:,:,3)
! transforming them to Fourier space
call xFFT3d(1,n1)
call xFFT3d(1,n2)
call xFFT3d(1,n )
! phase shifting them back using wrk0 and adding -0.5*(ik) to the RHS for the scalar
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
if (ialias(i,j,k) .gt. 1) then
wrk(i:i+1,j,k,n) = zip
else
! (u+)*(phi+) phase shifted back
r11 = wrk(i ,j,k,n1) * wrk(i,j,k,0) + wrk(i+1,j,k,n1) * wrk(i+1,j,k,0)
r12 = wrk(i+1,j,k,n1) * wrk(i,j,k,0) - wrk(i ,j,k,n1) * wrk(i+1,j,k,0)
! (v+)*(phi+) phase shifted back
r21 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
r22 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
! (w+)*(phi+) phase shifted back
r31 = wrk(i ,j,k,n ) * wrk(i,j,k,0) + wrk(i+1,j,k,n ) * wrk(i+1,j,k,0)
r32 = wrk(i+1,j,k,n ) * wrk(i,j,k,0) - wrk(i ,j,k,n ) * wrk(i+1,j,k,0)
! adding -0.5*(ik)*(the result) to the RHSs for the scalar
wrk(i ,j,k,n) = + 0.5d0 * ( akx(i+1)*r12 + aky(k)*r22 + akz(j)*r32 )
wrk(i+1,j,k,n) = - 0.5d0 * ( akx(i )*r11 + aky(k)*r21 + akz(j)*r31 )
end if
end do
end do
end do
end do phase_shifted_rhs
! at this moment wrk4...3+n_scalars+n_les contain the half of the convective term, which was
! obtained from the phase shifted quantities. Now we need to add the other half of the
! convective term (the one that is obtained by multiplication of no-phase-shifted stuff)
! first get the velocities into the real space and put them in wrk1...3
! this should remain in there untouched, to be used in rhs_velocity later
do n = 1,3
wrk(:,:,:,n) = fields(:,:,:,n)
call xFFT3d(-1,n)
end do
! now do scalars one at a time since we don't have enough storage to do them
! all at once
not_phase_shifted_rhs: do n = 4, 3 + n_scalars + n_les
! get the scalar into the real space and put it in the wrk(0)
wrk(:,:,:,0) = fields(:,:,:,n)
call xFFT3d(-1,0)
! calculate the product of the scalar with velocitiy components
wrk(:,:,:,n1) = wrk(:,:,:,0) * wrk(:,:,:,1)
wrk(:,:,:,n2) = wrk(:,:,:,0) * wrk(:,:,:,2)
wrk(:,:,:, 0) = wrk(:,:,:,0) * wrk(:,:,:,3)
! transform it to the Fourier space
call xFFT3d(1,n1)
call xFFT3d(1,n2)
call xFFT3d(1, 0)
! add the -0.5*(ik)*results to the RHS, along with the diffusion term
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
if (ialias(i,j,k) .lt. 2) then
rtmp1 = akx(i+1)*wrk(i+1,j,k,n1) + aky(k)*wrk(i+1,j,k,n2) + akz(j)*wrk(i+1,j,k,0)
rtmp2 = akx(i )*wrk(i ,j,k,n1) + aky(k)*wrk(i ,j,k,n2) + akz(j)*wrk(i ,j,k,0)
wnum2 = akx(i)**2 + aky(k)**2 + akz(j)**2
wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * rtmp1 - pe(n-3) * wnum2*fields(i ,j,k,n)
wrk(i+1,j,k,n) = wrk(i+1,j,k,n) - 0.5d0 * rtmp2 - pe(n-3) * wnum2*fields(i+1,j,k,n)
end if
end do
end do
end do
end do not_phase_shifted_rhs
! ------------------------------------------------------------
! Add the reaction rates to the RHS
!
! The LES-related scalars are not affected because the
! upper bound for n is 3+n_scalars, not 3+n_scalars+n_les
! -----------------------------------------------------------
reaction_rates: do n = 4, 3+n_scalars
if (scalar_type(n-3) .gt. 300) then
!!$ ! putting phase shifted scalar in wrk(n1)
!!$ do k = 1,nz
!!$ do j = 1,ny
!!$ do i = 1,nx+1,2
!!$ wrk(i ,j,k,n1) = fields(i ,j,k,n) * wrk(i,j,k,0) - fields(i+1,j,k,n) * wrk(i+1,j,k,0)
!!$ wrk(i+1,j,k,n1) = fields(i+1,j,k,n) * wrk(i,j,k,0) + fields(i ,j,k,n) * wrk(i+1,j,k,0)
!!$ end do
!!$ end do
!!$ end do
!!$ ! transforming it to real space
!!$ call xFFT3d(-1,n1)
!!$ ! putting the phase-shifted reaction rate in wrk(n2)
!!$ call scalar_reaction_rate(n1,n2)
!!$ ! transforming to Fourier space
!!$ call xFFT3d(1,n2)
!!$ ! phase shifting back and adding a half of it to the RHS
!!$ do k = 1,nz
!!$ do j = 1,ny
!!$ do i = 1,nx+1,2
!!$ if (ialias(i,j,k) .le. 1) then
!!$ ! phase shifting back
!!$ r11 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
!!$ r12 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
!!$ ! adding 0.5*(the result) to the RHSs for the scalar
!!$ wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * r11
!!$ wrk(i+1,j,k,n) = wrk(i+1,j,k,n) + 0.5d0 * r12
!!$ end if
!!$ end do
!!$ end do
!!$ end do
!!$
!!$ ! scond part: doing the same thing with not-phase-shifted scalar
!!$ ! putting it in wrk(n1)
!!$ wrk(:,:,:,n1) = fields(:,:,:,n)
!!$ ! transforming it to real space
!!$ call xFFT3d(-1,n1)
!!$ ! putting the phase-shifted reaction rate in wrk(n2)
!!$ call scalar_reaction_rate(n1,n2)
!!$ ! transforming to Fourier space
!!$ call xFFT3d(1,n2)
!!$ ! phase shifting back and adding a half of it to the RHS
!!$ do k = 1,nz
!!$ do j = 1,ny
!!$ do i = 1,nx+1,2
!!$ if (ialias(i,j,k) .le. 1) then
!!$ ! phase shifting back
!!$ r11 = wrk(i ,j,k,n2) * wrk(i,j,k,0) + wrk(i+1,j,k,n2) * wrk(i+1,j,k,0)
!!$ r12 = wrk(i+1,j,k,n2) * wrk(i,j,k,0) - wrk(i ,j,k,n2) * wrk(i+1,j,k,0)
!!$ ! adding 0.5*(the result) to the RHSs for the scalar
!!$ wrk(i ,j,k,n) = wrk(i ,j,k,n) + 0.5d0 * r11
!!$ wrk(i+1,j,k,n) = wrk(i+1,j,k,n) + 0.5d0 * r12
!!$ end if
!!$ end do
!!$ end do
!!$ end do
if (scalar_type(n-3) .eq. 311) then
! since the reaction rate is cubic, we need to apply some severe truncation.
! putting the scalar in wrk(n1)
wrk(:,:,:,n1) = fields(:,:,:,n)
! truncating it so only the modes with |k_i| < nx/4 remain
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
if (abs(akx(i)).gt.nx/4 .or. abs(aky(k)).gt.nx/4 .or. abs(akz(j)).gt.nx/4) then
wrk(i ,j,k,n1) = zip
wrk(i+1,j,k,n1) = zip
end if
end do
end do
end do
! transforming it to real space
call xFFT3d(-1,n1)
! self-adjusting bistable reaction needs the mean value of the scalar
rtmp1 = fields(1,1,1,n)/nxyz_all
call MPI_BCAST(rtmp1, 1, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
! getting the reaction rate
wrk(:,:,:,n2) = reac_sc(n-3) * (one - wrk(:,:,:,n1)**2) * (wrk(:,:,:,n1) - rtmp1)
! transforming it to Fourier space
call xFFT3d(1,n2)
! adding to the RHS
wrk(:,:,:,n) = wrk(:,:,:,n) + wrk(:,:,:,n2)
else
write(out,*) "The scalar type has a reaction rate that is not supported yet:", scalar_type(n-3)
call flush(out)
call my_exit(-1)
end if
end if
end do reaction_rates
end if phase_shifting_dealiasing
! special case - passive scalar with the uniform gradient as a source
! adding the source term - the first component of velocity, because we assume
! that the uniform gradient has slope 1 and direction in the x-direction
gradient_source: do n = 1, n_scalars
if (scalar_type(n) .eq. 0) wrk(:,:,:,n+3) = wrk(:,:,:,n+3) - fields(:,:,:,1)
end do gradient_source
!--------------------------------------------------------------------------------
! Add LES to the RHS of all the scalars
!--------------------------------------------------------------------------------
les_active: if (les) then
call les_rhs_scalars
end if les_active
return
end subroutine rhs_scalars
!================================================================================
!================================================================================
subroutine add_reaction(n)
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: n, rtype
real*8 :: scmean, rrate
! reaction type
rtype = scalar_type(n)/100
! raction rate
rrate = reac_sc(n)
select case (rtype)
case (1)
! KPP reaction rate
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2)
case (2)
! symmetric bistable
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2) * wrk(:,:,:,0)
case (3)
! self-adjusting bistable
scmean = fields(1,1,1,3+n)/nxyz_all
call MPI_BCAST(scmean, 1, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
wrk(:,:,:,0) = rrate * (1.d0 - wrk(:,:,:,0)**2) * &
(wrk(:,:,:,0) - scmean)
case default
write(out,*) "Unknown reaction rate"
call flush(out)
stop
end select
! FFT the reaction into the Fourier space
call xFFT3d(1,0)
! Adding reaction to the RHS in wrk(:,:,:,3+n)
wrk(:,:,:,3+n) = wrk(:,:,:,3+n) + wrk(:,:,:,0)
end subroutine add_reaction
!================================================================================
subroutine dealias_rhs(n)
use m_io
use m_parameters
use m_work
use x_fftw
implicit none
integer :: i, j, k, n
real*8 :: wnum2, akmax
akmax = real(kmax,8)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
if ( abs(akx(i)).gt.akmax .or. &
abs(aky(k)).gt.akmax .or. &
abs(akz(j)).gt.akmax ) then
wrk(i ,j,k,n) = zip
wrk(i+1,j,k,n) = zip
end if
end do
end do
end do
return
end subroutine dealias_rhs
!================================================================================
!================================================================================
!================================================================================
!================================================================================
subroutine test_rhs_scalars
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: i,j,k, n
real*8 :: a,b,c, x,y,z
if (task.eq.'hydro') then
a = 1.d0
b = 5.d0
c = 17.d0
do k = 1,nz
do j = 1,ny
do i = 1,nx
x = dx*real(i-1)
y = dx*real(j-1)
z = dx*real(myid*nz + k-1)
wrk(i,j,k,1) = sin(a * x)
wrk(i,j,k,2) = sin(b * y)
wrk(i,j,k,3) = sin(c * z)
wrk(i,j,k,4) = cos(a * x)
end do
end do
end do
do n = 1,4
call xFFT3d(1,n)
fields(:,:,:,n) = wrk(:,:,:,n)
end do
nu = .5d0
call rhs_scalars
print *,'got rhs'
call xFFT3d(-1,4)
do k = 1,nz
do j = 1,ny
do i = 1,nx
x = dx*real(i-1)
y = dx*real(j-1)
z = dx*real(myid*nz + k-1)
! checking
wrk(i,j,k,0) = -a*cos(2.*a*x) - cos(a*x)*(b*cos(b*y) + c*cos(c*z) + nu*a**2)
end do
end do
end do
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,4)
!!$ fname = 'r1.arr'
!!$ call write_tmp4
!!$
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,0)
!!$ fname = 'r0.arr'
!!$ call write_tmp4
wrk(:,:,:,0) = abs(wrk(1:nx,:,:,0) - wrk(1:nx,:,:,4))
print *,'Maximum error is ',maxval(wrk(1:nx,:,:,0))
!!$ tmp4(:,:,:) = wrk(1:nx,:,:,3) - wrk(1:nx,:,:,6)
!!$ fname = 'e3.arr'
!!$ call write_tmp4
!!$
end if
return
end subroutine test_rhs_scalars

553
rhs_velocity.f90 Normal file
View file

@ -0,0 +1,553 @@
subroutine rhs_velocity
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
use m_les
implicit none
integer :: i, j, k, n, nx3
real*8 :: t1(0:6), rtmp, wnum2, rnx3
! The IFFT of velocities has been done earlier in rhs_scalars
! the velocities were kept in wrk1...wrk3, intact.
!!$ ! putting the velocity field in the wrk array
!!$ wrk(:,:,:,1:3) = fields(:,:,:,1:3)
!!$ ! performing IFFT to convert them to the X-space
!!$ call xFFT3d(-1,1)
!!$ call xFFT3d(-1,2)
!!$ call xFFT3d(-1,3)
!-------------------------------------------------------------------------
! getting the Courant number (on the master process only)
wrk(:,:,:,4) = abs(wrk(:,:,:,1)) + abs(wrk(:,:,:,2)) + abs(wrk(:,:,:,3))
rtmp = maxval(wrk(1:nx,:,:,4))
call MPI_REDUCE(rtmp,courant,1,MPI_REAL8,MPI_MAX,0,MPI_COMM_TASK,mpi_err)
if (variable_dt) then
count = 1
call MPI_BCAST(courant,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
end if
courant = courant * dt / dx
!-------------------------------------------------------------------------
!--------------------------------------------------------------------------------
! Calculating the right-hand side for the velocities
!
! There are two options available: the standard 2/3 rule (dealias=0) and
! combination of phase shift and truncation (dealias=1). The latter retains
! more modes but requires more calculations thus slowing down the simulation.
! These are treated separately in two different "if" blocks. This is done in
! order not to complicate the logic. Also this way both blocks can be
! optimized separately.
!--------------------------------------------------------------------------------
two_thirds_rule: if (dealias.eq.0) then
! getting all 6 products of velocities
do k = 1,nz
do j = 1,ny
do i = 1,nx
t1(1) = wrk(i,j,k,1) * wrk(i,j,k,1)
t1(2) = wrk(i,j,k,1) * wrk(i,j,k,2)
t1(3) = wrk(i,j,k,1) * wrk(i,j,k,3)
t1(4) = wrk(i,j,k,2) * wrk(i,j,k,2)
t1(5) = wrk(i,j,k,2) * wrk(i,j,k,3)
t1(6) = wrk(i,j,k,3) * wrk(i,j,k,3)
do n = 1,6
wrk(i,j,k,n) = t1(n)
end do
end do
end do
end do
! converting the products to the Fourier space
do n = 1,6
call xFFT3d(1,n)
end do
! Building the RHS.
! First, put into wrk arrays the convectove terms (that will be multiplyed by "i"
! later) and the factor that corresponds to the diffusion
! Do not forget that in Fourier space the indicies are (ix, iz, iy)
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
t1(1) = - ( akx(i) * wrk(i,j,k,1) + aky(k) * wrk(i,j,k,2) + akz(j) * wrk(i,j,k,3) )
t1(2) = - ( akx(i) * wrk(i,j,k,2) + aky(k) * wrk(i,j,k,4) + akz(j) * wrk(i,j,k,5) )
t1(3) = - ( akx(i) * wrk(i,j,k,3) + aky(k) * wrk(i,j,k,5) + akz(j) * wrk(i,j,k,6) )
t1(4) = - nu * ( akx(i)**2 + aky(k)**2 + akz(j)**2 )
do n = 1,4
wrk(i,j,k,n) = t1(n)
end do
end do
end do
end do
! now take the actual fields from fields(:,:,:,:) and calculate the RHSs
! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
! which are not yet multiplied by "i"
! wrk(:,:,:,4) contains the Laplace operator in Fourier space. To get the diffusion term
! we need to take wrk(:,:,:,4) and multiply it by the velocity
t1(6) = real(kmax,8)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! If the dealiasing option is 2/3-rule (dealias=0) then we retain the modes
! inside the cube described by $| k_i | \leq k_{max}$, $i=1,2,3$.
! The rest of the modes is purged
if (ialias(i,j,k) .gt. 0) then
! setting the Fourier components to zero
wrk(i ,j,k,1:3) = zip
wrk(i+1,j,k,1:3) = zip
else
! RHS for u, v and w
do n = 1,3
! taking the convective term, multiply it by "i"
! (see how it's done in x_fftw.f90)
! and adding the diffusion term
rtmp = - wrk(i+1,j,k,n) + wrk(i ,j,k,4) * fields(i ,j,k,n)
wrk(i+1,j,k,n) = wrk(i ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n)
wrk(i ,j,k,n) = rtmp
end do
end if
end do
end do
end do
end if two_thirds_rule
!--------------------------------------------------------------------------------
! The second option (dealias=1). All pairwise products of velocities are
! dealiased using one phase shift of (dx/2,dy/2,dz/2).
!--------------------------------------------------------------------------------
phase_shifting: if (dealias.eq.1) then
! work parameters
wrk(:,:,:,0) = zip
! getting all 6 products of velocities
do k = 1,nz
do j = 1,ny
do i = 1,nx
t1(1) = wrk(i,j,k,1) * wrk(i,j,k,1)
t1(2) = wrk(i,j,k,1) * wrk(i,j,k,2)
t1(3) = wrk(i,j,k,1) * wrk(i,j,k,3)
t1(4) = wrk(i,j,k,2) * wrk(i,j,k,2)
t1(5) = wrk(i,j,k,2) * wrk(i,j,k,3)
t1(6) = wrk(i,j,k,3) * wrk(i,j,k,3)
do n = 1,6
wrk(i,j,k,n) = t1(n)
end do
end do
end do
end do
! converting the products to the Fourier space
do n = 1,6
call xFFT3d(1,n)
end do
! Building the RHS.
! First, put into wrk arrays the convectove terms (that will be multiplyed by "i"
! later) and the factor that corresponds to the diffusion
! Do not forget that in Fourier space the indicies are (ix, iz, iy)
do k = 1,nz
do j = 1,ny
do i = 1,nx+2
t1(1) = - ( akx(i) * wrk(i,j,k,1) + aky(k) * wrk(i,j,k,2) + akz(j) * wrk(i,j,k,3) )
t1(2) = - ( akx(i) * wrk(i,j,k,2) + aky(k) * wrk(i,j,k,4) + akz(j) * wrk(i,j,k,5) )
t1(3) = - ( akx(i) * wrk(i,j,k,3) + aky(k) * wrk(i,j,k,5) + akz(j) * wrk(i,j,k,6) )
! putting a factor from the diffusion term into t1(4) (and later in wrk4)
t1(4) = - nu * ( akx(i)**2 + aky(k)**2 + akz(j)**2 )
do n = 1,4
wrk(i,j,k,n) = t1(n)
end do
end do
end do
end do
! now use the actual fields from fields(:,:,:,:) to calculate the RHSs
! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
! which are not yet multiplied by "i"
! wrk(:,:,:,4) contains the Laplace operator in Fourier space. To get the diffusion term
! we need to take wrk(:,:,:,4) and multiply it by the velocity
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! If the dealiasing option is (dealias=1) then we retain the modes
! for which no more than one component of the k-vector is larger than nx/3.
! The rest of the modes is purged.
if (ialias(i,j,k) .gt. 1) then
! setting the Fourier components to zero
wrk(i ,j,k,1:3) = zip
wrk(i+1,j,k,1:3) = zip
else
! RHS for u, v and w
do n = 1,3
! taking the HALF of the convective term, multiply it by "i"
! and adding the diffusion term
rtmp = - 0.5d0 * wrk(i+1,j,k,n) + wrk(i ,j,k,4) * fields(i ,j,k,n)
wrk(i+1,j,k,n) = 0.5d0 * wrk(i ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n)
wrk(i ,j,k,n) = rtmp
end do
end if
end do
end do
end do
!--------------------------------------------------------------------------------
! Second part of the phase shifting technique
!--------------------------------------------------------------------------------
! since wrk1...3 are taken by parts of RHS constructed earlier, we can use
! only wrk0 and wrk4...6.
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! computing sines and cosines for the phase shift of dx/2,dy/2,dz/2
! and putting them into wrk0
wrk(i ,j,k,0) = cos(half*(akx(i )+aky(k)+akz(j))*dx)
wrk(i+1,j,k,0) = sin(half*(akx(i+1)+aky(k)+akz(j))*dx)
! wrk4 will have phase-shifted u
wrk(i ,j,k,4) = fields(i ,j,k,1) * wrk(i,j,k,0) - fields(i+1,j,k,1) * wrk(i+1,j,k,0)
wrk(i+1,j,k,4) = fields(i+1,j,k,1) * wrk(i,j,k,0) + fields(i ,j,k,1) * wrk(i+1,j,k,0)
! wrk5 will have phase-shifted v
wrk(i ,j,k,5) = fields(i ,j,k,2) * wrk(i,j,k,0) - fields(i+1,j,k,2) * wrk(i+1,j,k,0)
wrk(i+1,j,k,5) = fields(i+1,j,k,2) * wrk(i,j,k,0) + fields(i ,j,k,2) * wrk(i+1,j,k,0)
end do
end do
end do
! transforming u+ and v+ into X-space
call xFFT3d(-1,4)
call xFFT3d(-1,5)
! now wrk4 and wrk5 contain u+ and v+
! getting (u+)*(u+) in real space, converting it to Fourier space,
! phase shifting back and adding -0.5*(the results) to the RHS for u
wrk(:,:,:,6) = wrk(:,:,:,4)**2
call xFFT3d(1,6)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
wrk(i ,j,k,6) = rtmp
end do
end do
end do
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * akx(i+1) * wrk(i+1,j,k,6)
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * akx(i ) * wrk(i ,j,k,6)
end do
end do
end do
! getting (u+)*(v+) in real space, converting it to Fourier space,
! phase shifting back and adding -0.5*(the results) to the RHSs for u and v
wrk(:,:,:,6) = wrk(:,:,:,4)*wrk(:,:,:,5)
call xFFT3d(1,6)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
wrk(i ,j,k,6) = rtmp
end do
end do
end do
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * aky(k) * wrk(i+1,j,k,6)
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * aky(k) * wrk(i ,j,k,6)
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * akx(i+1) * wrk(i+1,j,k,6)
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * akx(i ) * wrk(i ,j,k,6)
end do
end do
end do
! getting (v+)*(v+) in real space, converting it to Fourier space,
! phase shifting back and adding -0.5*(the results) to the RHS for v
wrk(:,:,:,6) = wrk(:,:,:,5)**2
call xFFT3d(1,6)
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
wrk(i ,j,k,6) = rtmp
end do
end do
end do
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * aky(k) * wrk(i+1,j,k,6)
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * aky(k) * wrk(i ,j,k,6)
end do
end do
end do
! now get the (w+) in wrk6
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! wrk6 will have phase-shifted w
wrk(i ,j,k,6) = fields(i ,j,k,3) * wrk(i,j,k,0) - fields(i+1,j,k,3) * wrk(i+1,j,k,0)
wrk(i+1,j,k,6) = fields(i+1,j,k,3) * wrk(i,j,k,0) + fields(i ,j,k,3) * wrk(i+1,j,k,0)
end do
end do
end do
! transforming w+ into X-space
call xFFT3d(-1,6)
! at this point wrk4..6 contain (u+), (v+) and (w+) in real space.
! the combinations that we have not dealt with are: uw, vw and ww.
! we'll deal with all three of them at once.
! first get all three of these in wrk4...6 and
wrk(:,:,:,4) = wrk(:,:,:,4) * wrk(:,:,:,6)
wrk(:,:,:,5) = wrk(:,:,:,5) * wrk(:,:,:,6)
wrk(:,:,:,6) = wrk(:,:,:,6)**2
! transform them into Fourier space
call xFFT3d(1,4)
call xFFT3d(1,5)
call xFFT3d(1,6)
! phase shift back to origianl grid and add to corresponding RHSs
do n = 4,6
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
rtmp = wrk(i ,j,k,n) * wrk(i,j,k,0) + wrk(i+1,j,k,n) * wrk(i+1,j,k,0)
wrk(i+1,j,k,n) = wrk(i+1,j,k,n) * wrk(i,j,k,0) - wrk(i ,j,k,n) * wrk(i+1,j,k,0)
wrk(i ,j,k,n) = rtmp
end do
end do
end do
end do
! adding to corresponding RHSs
do k = 1,nz
do j = 1,ny
do i = 1,nx+1,2
! If the dealiasing option is (dealias=1) then we retain the modes
! for which no more than one component of the k-vector is larger than nx/3.
! The rest of the modes is purged.
if (ialias(i,j,k) .lt. 2) then
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * akz(j) * wrk(i+1,j,k,4)
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * akz(j) * wrk(i ,j,k,4)
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * akz(j) * wrk(i+1,j,k,5)
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * akz(j) * wrk(i ,j,k,5)
wrk(i ,j,k,3) = wrk(i ,j,k,3) + 0.5d0 * &
(akx(i+1)*wrk(i+1,j,k,4) + aky(k)*wrk(i+1,j,k,5) + akz(j)*wrk(i+1,j,k,6))
wrk(i+1,j,k,3) = wrk(i+1,j,k,3) - 0.5d0 * &
(akx(i )*wrk(i ,j,k,4) + aky(k)*wrk(i ,j,k,5) + akz(j)*wrk(i ,j,k,6))
else
wrk(i:i+1,j,k,1) = zip
wrk(i:i+1,j,k,2) = zip
wrk(i:i+1,j,k,3) = zip
end if
end do
end do
end do
end if phase_shifting
! if performing large eddy simulations, call LES subroutine to augment
! the right hand side for velocioties
les_active: if (les) then
call les_rhs_velocity
end if les_active
return
end subroutine rhs_velocity
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
!================================================================================
subroutine test_rhs_velocity
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
integer :: i,j,k, n
real*8 :: a,b,c, x,y,z
! defining very particular velocities so the RHS can be computed analytically
if (task.eq.'hydro') then
write(out,*) 'inside.'
call flush(out)
a = 1.d0
b = 1.d0
c = 1.d0
do k = 1,nz
do j = 1,ny
do i = 1,nx
x = dx*real(i-1)
y = dx*real(j-1)
z = dx*real(myid*nz + k-1)
wrk(i,j,k,1) = sin(a * x)
wrk(i,j,k,2) = sin(b * y)
wrk(i,j,k,3) = sin(c * z)
end do
end do
end do
write(out,*) 'did work'
call flush(out)
do n = 1,3
call xFFT3d(1,n)
fields(:,:,:,n) = wrk(:,:,:,n)
end do
write(out,*) 'did FFTs'
call flush(out)
nu = 0.d0
call rhs_velocity
write(out,*) 'got rhs'
call flush(out)
do n = 1,3
call xFFT3d(-1,n)
end do
write(out,*) 'did FFTs'
call flush(out)
do k = 1,nz
do j = 1,ny
do i = 1,nx
x = dx*real(i-1)
y = dx*real(j-1)
z = dx*real(myid*nz + k-1)
! checking u
wrk(i,j,k,4) = -sin(a*x) * ( two*a*cos(a*x) + b*cos(b*y) + c*cos(c*z) + nu*a**2)
! checking v
wrk(i,j,k,5) = -sin(b*y) * ( two*b*cos(b*y) + a*cos(a*x) + c*cos(c*z) + nu*b**2)
! checking w
wrk(i,j,k,6) = -sin(c*z) * ( two*c*cos(c*z) + b*cos(b*y) + a*cos(a*x) + nu*c**2)
end do
end do
end do
!!$ do k = 1,nz
!!$ write(out,"(3e15.6)") wrk(1,1,k,3),wrk(1,1,k,5),wrk(1,1,k,4)
!!$ end do
wrk(:,:,:,0) = &
abs(wrk(:,:,:,1) - wrk(:,:,:,4)) + &
abs(wrk(:,:,:,2) - wrk(:,:,:,5)) + &
abs(wrk(:,:,:,3) - wrk(:,:,:,6))
print *,'Maximum error is ',maxval(wrk(1:nx,:,:,0))
tmp4(:,:,:) = wrk(1:nx,:,:,1) - wrk(1:nx,:,:,4)
fname = 'e1.arr'
call write_tmp4
tmp4(:,:,:) = wrk(1:nx,:,:,2) - wrk(1:nx,:,:,5)
fname = 'e2.arr'
call write_tmp4
tmp4(:,:,:) = wrk(1:nx,:,:,3) - wrk(1:nx,:,:,6)
fname = 'e3.arr'
call write_tmp4
end if
return
end subroutine test_rhs_velocity

159
velocity_rescale.f90 Normal file
View file

@ -0,0 +1,159 @@
!================================================================================
! Subroutine that rescales the current velocities. The rescales velocity
! will have the spectrum that is defined in the input file via parameters
! isp_type (spectrum type) and peak_wavenum (peak wavenumber).
!
! The program is copies from a part of init_velocity.f90. Some time in the
! future we should make it one routine that is called in init_velocity.
!
! Time-stamp: <2010-01-25 17:01:55 (chumakov)>
!================================================================================
subroutine velocity_rescale
use m_openmpi
use m_parameters
use m_io
use m_fields
use x_fftw
implicit none
integer :: i, j, k, n
real, allocatable :: rr(:)
real*8, allocatable :: e_spec(:), e_spec1(:)
integer *8, allocatable :: hits(:), hits1(:)
integer :: n_shell
real*8 :: sc_rad1, sc_rad2
real*8 :: wmag, wmag2, ratio, fac
! if Taylor-Green, return
if (isp_type.eq.-1) return
!================================================================================
allocate( e_spec(kmax), e_spec1(kmax), rr(nx+2), hits(kmax), hits1(kmax), stat=ierr)
if (ierr.ne.0) stop "cannot allocate the init_velocity arrays"
write(out,*) 'Rescaling the velocities'
call flush(out)
!-------------------------------------------------------------------------------
! Making the spectrum to be what is prescribed in the input file <...>.in
!-------------------------------------------------------------------------------
! --- first get the energy spectrum (copied from m_stat.f90)
! need this normalization factor because the FFT is unnormalized
fac = one / real(nx*ny*nz_all)**2
! zeroing out the arrays
e_spec1 = zip
e_spec = zip
hits = 0
hits1 = 0
! finding the total energy in each shell and number of hits in each shell
! storing them in arrays hits1 and e_spec1
do k = 1,nz
do j = 1,ny
do i = 1,nx
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) then
hits1(n_shell) = hits1(n_shell) + 1
e_spec1(n_shell) = e_spec1(n_shell) + &
fac * (fields(i,j,k,1)**2 + fields(i,j,k,2)**2 + fields(i,j,k,3)**2)
end if
end do
end do
end do
! reducing the number of hits and energy to two arrays on master node
count = kmax
call MPI_REDUCE(hits1,hits,count,MPI_INTEGER8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
count = kmax
call MPI_REDUCE(e_spec1,e_spec,count,MPI_REAL8,MPI_SUM,0,MPI_COMM_TASK,mpi_err)
! now the master node counts the energy density in each shell
if (myid.eq.0) then
fac = four/three * PI / two
do k = 1,kmax
sc_rad1 = real(k,8) + half
sc_rad2 = real(k,8) - half
if (k.eq.1) sc_rad2 = 0.d0
if (hits(k).gt.0) then
e_spec(k) = e_spec(k) / hits(k) * fac * (sc_rad1**3 - sc_rad2**3)
else
e_spec(k) = zip
end if
end do
end if
! broadcasting the spectrum
count = kmax
call MPI_BCAST(e_spec,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
!-------------------------------------------------------------------------------
! Now make the spectrum to be as desired
!-------------------------------------------------------------------------------
! first, define the desired spectrum
do k = 1,kmax
wmag = real(k, 8)
ratio = wmag / peak_wavenum
if (isp_type.eq.0) then
! Plain Kolmogorov spectrum
e_spec1(k) = wmag**(-5.d0/3.d0)
else if (isp_type.eq.1) then
! Exponential spectrum
e_spec1(k) = ratio**3 / peak_wavenum * exp(-3.0D0*ratio)
else if (isp_type.eq.3) then
! Von Karman spectrum
fac = two * PI * ratio
e_spec1(k) = fac**4 / (one + fac**2)**3
else
write(out,*) "ERROR: WRONG INITIAL SPECTRUM TYPE: ",isp_type
call flush(out)
stop
end if
end do
! normalize it so it has the unit total energy
e_spec1 = e_spec1 / sum(e_spec1(1:kmax))
! now go over all Fourier shells and multiply the velocities in a shell by
! the sqrt of ratio of the resired to the current spectrum
do k = 1,nz
do j = 1,ny
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
fields(i,j,k,1:3) = fields(i,j,k,1:3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
else
fields(i,j,k,1:3) = zip
end if
end do
end do
end do
write(out,*) "Rescaled the velocities."
call flush(out)
! deallocate work arrays
deallocate(e_spec, e_spec1, rr, hits, hits1, stat=ierr)
return
end subroutine velocity_rescale

73
write_tmp4.f90 Normal file
View file

@ -0,0 +1,73 @@
subroutine write_tmp4
use m_parameters
use m_io
use m_work
implicit none
integer :: my_out=13, i,j,k
integer*4 :: sizes(3)
!======================================================================
! --- defining the size of whole array
sizes(1)=nx; sizes(2)=ny; sizes(3)=nz*numprocs;
count = nx*ny*nz
if (myid.ne.master) then
id_to = master
tag = myid
call MPI_SEND(tmp4,count,MPI_REAL4,master,tag,MPI_COMM_TASK,mpi_err)
else
!! open(my_out,file=fname,form='binary')
open(my_out,file=fname,form='unformatted', access='stream')
write(my_out) sizes(1:3)
write(my_out) (((tmp4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
do id_from=1,numprocs-1
tag = id_from
call MPI_RECV(tmp4,count,MPI_REAL4,id_from,tag,MPI_COMM_TASK,mpi_status,mpi_err)
write(my_out) (((tmp4(i,j,k),i=1,nx),j=1,ny),k=1,nz)
end do
close(my_out)
end if
!======================================================================
return
end subroutine write_tmp4
!======================================================================
subroutine write_tmp4_all
use m_parameters
use m_io
use m_work
implicit none
integer :: my_out=13, i,j,k
integer(kind=MPI_INTEGER_KIND) :: sizes(3), fh
integer(kind=MPI_OFFSET_KIND) :: offset
!======================================================================
! --- defining the size of whole array
sizes(1)=nx; sizes(2)=ny; sizes(3)=nz*numprocs;
! --- writing into the file with appropriate offset
call MPI_INFO_CREATE(mpi_info,mpi_err)
if(ierr.ne.0) stop '*** WRITE_TMP4_ALL: cannot create mpi_info'
call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_WRONLY+MPI_MODE_CREATE,mpi_info,fh,mpi_err)
if (myid.eq.0) call MPI_FILE_WRITE_AT(fh,0,sizes,3,MPI_INTEGER4,mpi_status,mpi_err)
offset = 12 + myid*nx*ny*nz * 4
count = nx * ny * nz
call MPI_FILE_WRITE_AT_ALL(fh,offset,tmp4,count,MPI_REAL4,mpi_status,mpi_err)
call MPI_FILE_CLOSE(fh,mpi_err)
call MPI_INFO_FREE(mpi_info,mpi_err)
!======================================================================
return
end subroutine write_tmp4_all

1286
x_fftw.f90 Normal file
View file

File diff suppressed because it is too large Load diff