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

648 lines
18 KiB
Fortran
Raw Permalink Blame History

!================================================================================
! 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 :: lx,ly,lz,spx,spy,spz !mx
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, fdm_sw
real*8 :: TMAX, TRESCALE, TSCALAR, RE, nu, dt
! now many times to rescale teh velocities
integer :: NRESCALE
integer :: flow_type
real*8 :: t_length,delta_k
real*8 :: t_u_prime,c_f
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)
!! For model spectrum (4)-------------------------------------------201403
real*8, allocatable :: model_e_spec(:)
real*8 :: model_A,k0
!!----------------------------------------------------------------------
! FDM variables
integer :: swadtv,restartnum
integer :: inx1,inx2,yrsw
real*8 :: ac,bc,c_cut,c_ref,min_wr,prof_wr,hx,hy,diff,pre
real*8 :: dt_fullsave,dt_fdmsave,t_fullsave,t_fdmsave
real*8 :: fdmcfl,startad,dummyu_,schmidt
real*8 :: bgn_save
! files for statistics
integer, parameter :: tp_stat = 690
integer, parameter :: tp_field = 101
integer, parameter :: tp_plane = 506
integer, parameter :: model_spectrum = 404
integer, parameter :: tp_spec = 501
integer, parameter :: tp_tke_spec = 502
integer, parameter :: tp_eps_spec = 503
integer, parameter :: tp_max_min_uvw = 504
integer, parameter :: St_data = 505
integer, parameter :: restart_pair = 506
!================================================================================
contains
!================================================================================
subroutine m_parameters_init
implicit none
integer :: mx
call get_run_name
! constants
PI = four * atan(one)
TWO_PI = two * PI
! switches
int_scalars = .false.
call read_input_file
! maximum resolved wavenumber
mx = 2 * nx / lx ! Nx/2PI
if (dealias.eq.0) then
kmax = mx/3
elseif (dealias.eq.1) then
kmax = floor(real(mx,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,'(a30,3x,i5)')"SUB_m_parameters_init 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
real*8 :: t_tke
! 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
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,'(a30,3x,3i5)') 'NX,NY,NZ_ALL',nx,ny,nz_all
write(out,'(a30,3x,3i5)') 'NX,NY,NZ',nx,ny,nz
call flush(out)
read(in,*,ERR=9000) lx,ly,lz
write(out,'(a30,3x,3i5)') 'LX,LY,LZ[PI]',lx,ly,lz
dx = PI * dble(lx) / dble(nx)
dy = PI * dble(ly) / dble(ny)
dz = PI * dble(lz) / dble(nz_all)
call flush(out)
! -------------------------------------------------------------
read(in,*)
read(in,*,ERR=9000,END=9000) ITMIN
write(out,'(a30,3x,i5)') 'ITMIN =',ITMIN
last_dump = ITMIN
read(in,*,ERR=9000,END=9000) IPRINT1
write(out,'(a30,3x,i5)') 'IPRINT1 =',IPRINT1
read(in,*,ERR=9000,END=9000) IPRINT2
write(out,'(a30,3x,i10)') 'IPRINT2 =',IPRINT2
read(in,*,ERR=9000,END=9000) IWRITE4
write(out,'(a30,3x,i5)') 'IWRITE4 =',IWRITE4
read(in,*,ERR=9000,END=9000) spx,spy,spz
write(out,'(a30,3x,3i5)') 'spx,spy,spz',spx,spy,spz
read(in,*,ERR=9000,END=9000) fdm_sw
write(out,'(a30,3x,i5)') 'fdm_sw= ',fdm_sw
read(in,*)
write(out,"(70('-'))")
call flush(out)
read(in,*,ERR=9000,END=9000) TMAX
write(out,'(a30,3x,e15.5)') 'TMAX =',TMAX
read(in,*,ERR=8000,END=9000) TRESCALE, NRESCALE
100 write(out,'(a30,3x,f10.4,i5)') 'TRESCALE, NRESCALE',TRESCALE,NRESCALE
read(in,*,ERR=9000,END=9000) TSCALAR
write(out,'(a30,3x,e15.5)') 'TSCALAR =',TSCALAR
read(in,*)
write(out,"(70('-'))")
call flush(out)
read(in,*,ERR=9000,END=9000) flow_type
write(out,'(a30,3x,i5)') 'flow_type =',flow_type
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) nu
write(out,'(a30,3x,f10.5)') 'nu(kinematic vis) =',nu
RE = 1.0d0/nu
read(in,*,ERR=9000,END=9000) DT
write(out,'(a30,3x,f10.5)') 'DT =',DT
if (dt.lt.0.0d0) then
variable_dt = .false.
dt = -dt
else
variable_dt = .true.
end if
read(in,*)
write(out,"(70('-'))")
ITMAX=int(TMAX/dt)
write(out,'(a30,3x,i8)') 'ITMAX =',ITMAX
call flush(out)
!ksj*************************************************************
read(in,*,ERR=9000,END=9000) t_length
write(out,'(a30,3x,f10.5)') 'target integral length =',t_length
read(in,*,ERR=9000,END=9000) t_u_prime
write(out,'(a30,3x,f10.5)') 'target u` =',t_u_prime
t_tke = three/two*(t_u_prime**2.0)
write(out,'(a30,3x,f10.5)') 'target tke =',t_tke
read(in,*,ERR=9000,END=9000) delta_k
write(out,'(a30,3x,f10.5)') 'del_k =',delta_k
read(in,*)
write(out,"(70('-'))")
call flush(out)
!****************************************************************
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) isp_type
write(out,'(a30,3x,i5)') 'isp_type =',isp_type
read(in,*,ERR=9000,END=9000) ir_exp
write(out,'(a30,3x,i5)') 'ir_exp =',ir_exp
read(in,*,ERR=9000,END=9000) peak_wavenum
write(out,'(a30,3x,f10.5)') 'peak_wavenum =',peak_wavenum
read(in,*)
write(out,"(70('-'))")
call flush(out)
! ------------------------------------------------------------
read(in,*,ERR=9000,END=9000) force_type
write(out,'(a30,3x,i5)') 'force_type =',force_type
read(in,*,ERR=9000,END=9000) kfmax
write(out,'(a30,3x,i5)') 'kfmax =',kfmax
read(in,*,ERR=9000,END=9000) c_f
write(out,'(a30,3x,f10.5)') 'c_f =',c_f
famp=1.d0/c_f
read(in,*)
write(out,"(70('-'))")
call flush(out)
read(in,*,ERR=9000,END=9000) dealias
write(out,'(a30,3x,i5)') 'dealias =',dealias
read(in,*)
write(out,"(70('-'))")
call flush(out)
! -------------------------------------------------------------
read(in,*,ERR=9000,END=9000) det_rand
write(out,'(a30,3x,i5)') 'det_rand =',det_rand
read(in,*,ERR=9000,END=9000) RN1
write(out,'(a30,3x,e20.10)') 'RN1 =',RN1
read(in,*,ERR=9000,END=9000) RN2
write(out,'(a30,3x,e20.10)') 'RN2 =',RN2
read(in,*,ERR=9000,END=9000) RN3
write(out,'(a30,3x,e20.10)') '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,'(a30,3x,i5)') 'nptot =',nptot
read(in,*,ERR=9000,END=9000) particles_tracking_scheme
write(out,'(a30,3x,i5)') '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,'(a30,3x,e20.10)') 'starttime_particles: ',starttime_particles
read(in,*,ERR=9000,END=9000) particles_filter_size
write(out,'(a30,3x,e20.10)') '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,'(a30,3x,i5)') '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,'(a30,3x,i5)') '# 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
! -------------------------------------------------------------
! FDM input
read(in,*)
read(in,*)
read(in,*)
read(in,*)
read(in,*)
read(in,*,ERR=9000,END=9000) bgn_save
write(out,'(a30,3x,f8.2)') 'bgn_save =',bgn_save
read(in,*,ERR=9000,END=9000) swadtv
write(out,'(a30,3x,i5)') 'swadtv =',swadtv
read(in,*,ERR=9000,END=9000) restartnum
write(out,'(a30,3x,i5)') 'restartnum =',restartnum
read(in,*)
write(out,"(70('-'))")
call flush(out)
read(in,*,ERR=9000,END=9000) dt_fdmsave
write(out,'(a30,3x,e20.10)') 'dt_fdmsave =',dt_fdmsave
read(in,*,ERR=9000,END=9000) dt_fullsave
write(out,'(a30,3x,e20.10)') 'dt_fullsave =',dt_fullsave
read(in,*,ERR=9000,END=9000) schmidt
write(out,'(a30,3x,e20.10)') 'Sc# =',schmidt
read(in,*,ERR=9000,END=9000) pre
write(out,'(a30,3x,e20.10)') 'pre =',pre
read(in,*,ERR=9000,END=9000) ac
write(out,'(a30,3x,e20.10)') 'ac =',ac
read(in,*,ERR=9000,END=9000) bc
write(out,'(a30,3x,e20.10)') 'bc =',bc
read(in,*,ERR=9000,END=9000) c_cut
write(out,'(a30,3x,e20.10)') 'c_cut =',c_cut
read(in,*,ERR=9000,END=9000) c_ref
write(out,'(a30,3x,e20.10)') 'c_ref =',c_ref
read(in,*,ERR=9000,END=9000) min_wr
write(out,'(a30,3x,e20.10)') 'min_wr =',min_wr
read(in,*,ERR=9000,END=9000) prof_wr
write(out,'(a30,3x,e20.10)') 'prof_wr =',prof_wr
read(in,*,ERR=9000,END=9000) dummyu_
write(out,'(a30,3x,e20.10)') 'dummyu_ =',dummyu_
read(in,*,ERR=9000,END=9000) startad
write(out,'(a30,3x,e20.10)') 'startad =',startad
read(in,*,ERR=9000,END=9000) inx1
write(out,'(a30,3x,i5)') 'inx1 =',inx1
read(in,*,ERR=9000,END=9000) inx2
write(out,'(a30,3x,i5)') 'inx2 =',inx2
read(in,*,ERR=9000,END=9000) yrsw
if (inx1.gt.0) then
if (yrsw.eq.0) then
yrsw=1
elseif (yrsw.eq.1) then
yrsw=inx1
endif
endif
write(out,'(a30,3x,i5)') 'yrsw =',yrsw
read(in,*,ERR=9000,END=9000) fdmcfl
write(out,'(a30,3x,e20.10)') 'FDM_CFL =',fdmcfl
! initialize ----------------------
t_fdmsave=dt_fdmsave
t_fullsave=dt_fullsave
diff=nu/schmidt
hx=dx
hy=hx
! 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
Reference in a new issue View git blame Copy permalink