! $Id$
!
!***********************************************************************
!
! The Pencil Code is a high-order finite-difference code for compressible
! hydrodynamic flows with magnetic fields and particles. It is highly
! modular and can easily be adapted to different types of problems.
!
! MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
! MMMMMMM7MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
! MMMMMMMIMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
! MMMMMMMIIMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
! MMMMMMM7IMMMMMMMMMMMMMMMMMMM7IMMMMMMMMMMMMMMMMMMMMDMMMMMMM
! MMMMMMZIIMMMMMMMMMMMMMMMMMMMIIMMMMMMMMMMMMMMMMMMMMIMMMMMMM
! MMMMMMIIIZMMMMMMMMMMMMMMMMMMIIMMMMMMMMMMMMMMMMMMMMIMMMMMMM
! MMMMMMIIIIMMMMMMMMMMMMMMMMMNII$MMMMMMMMMMMMMMMMMM$IMMMMMMM
! MMMMM8IIIIMMMMMMMMMMMMMMMMM$IIIMMMMMMMMMMMMMMMMMMII7MMMMMM
! MMMMMD7II7MMMMMMMMMMMMMMMMMIIIIMMMMMMMMMMMMMMMMMMIIIMMMMMM
! MMMMN..:~=ZMMMMMMMMMMMMMMMMIIIIDMMMMMMMMMMMMMMMMDIIIMMMMMM
! MMMM8.,:~=?MMMMMMMMMMMMMMMMOII7NMMMMMMMMMMMMMMMMZIIIDMMMMM
! MMMM. ,::=+MMMMMMMMMMMMMMM..,~=?MMMMMMMMMMMMMMMMIIII$MMMMM
! MMMM..,:~=+MMMMMMMMMMMMMM8 .,~=+DMMMMMMMMMMMMMMM8II78MMMMM
! MMMM .,:~==?MMMMMMMMMMMMMN .,~=+NMMMMMMMMMMMMMM8 ,~~+MMMMM
! MMM7 ,:~+=?MMMMMMMMMMMMM .,~==?MMMMMMMMMMMMMM..,~~+DMMMM
! MMM. ,:~==?MMMMMMMMMMMMM .,~~=?MMMMMMMMMMMMMM. ,~~+?MMMM
! MMM. ,:~~=??MMMMMMMMMMMM ,,~~=?MMMMMMMMMMMMM8 .,~~=?NMMM
! MMM. ,:~~=+?MMMMMMMMMMMI ,,~~=+?MMMMMMMMMMMM. .,~~=?MMMM
! MM~. .,:~:==?MMMMMMMMMMM ,,~~==?MMMMMMMMMMMM .,~~=+?MMM
! MMN8 .D,D+=M=8MMMMMMMMMM ,,~~==?MMMMMMMMMMMM. .,~~==?MMM
! MM==8. .8===MMMMMMMMMM .,,~~==?$MMMMMMMMMM= .,~~==?MMM
! MM==D. +===MMMMMMMMMO$ .I?7$=7=NMMMMMMMMMM. ,,~~==?$MM
! MM==D. +===MMMMMMMMM==O?.. ====MMMMMMMMMM. ,,~~==??MM
! MM==D. +===MMMMMMMMM===. .===MMMMMMMMMM+$.?=,7==?7MM
! MM==D. +===MMMMMMMMM===. .===MMMMMMMMMZ==8 .8==MM
! MM==D. +===MMMMMMMMM===. .===MMMMMMMMMZ==I .Z==MM
! MM==D. . +===MMMMMMMMM===.... .===MMMMMMMMMZ==I .Z==MM
! MM==D. +===MMMMMMMMM===. .===MMMMMMMMMZ==I .Z==MM
! MM==D. +===MMMMMMMMM===. .===MMMMMMMMMZ==I .Z==MM
! MM==D. +===MMMMMMMMM===.. .===MMMMMMMMMZ==I .Z==MM
! MM==D. . +===MMMMMMMMM===.... .===MMMMMMMMMZ==I .Z==MM
!
! More information can be found in the Pencil Code manual and at the
! website http://www.nordita.org/software/pencil-code/.
!
!***********************************************************************
program start
!
use Boundcond, only: update_ghosts, initialize_boundcond
use Cdata
use Chemistry, only: init_chemistry
use Chiral, only: init_chiral
use Cosmicrayflux, only: init_fcr
use Cosmicray, only: init_ecr
use Density, only: init_lnrho, write_z_stratification
use Diagnostics
use Dustdensity, only: init_nd
use Dustvelocity, only: init_uud
use Energy, only: init_energy
use EquationOfState, only: init_eos
use FArrayManager, only: farray_clean_up
use Filter
use General
use Gravity, only: init_gg
use Grid
use HDF5_IO, only: init_hdf5, initialize_hdf5
use Hydro, only: init_uu
use Hyperresi_strict, only: hyperresistivity_strict
use Hypervisc_strict, only: hyperviscosity_strict
use ImplicitPhysics, only: init_implicit_physics
use Initcond
use InitialCondition, only: initial_condition_all, initial_condition_clean_up
use Interstellar, only: init_interstellar
use IO, only: wgrid, directory_names, wproc_bounds, output_globals
use HDF5_IO, only: wdim
use Lorenz_gauge, only: init_lorenz_gauge
use Magnetic, only: init_aa
use Messages
use Mpicomm
use NeutralDensity, only: init_lnrhon
use NeutralVelocity, only: init_uun
use OMP_lib
use Param_IO
use Particles_main
use Polymer, only: init_poly
use PointMasses
use Pscalar, only: init_lncc
use Radiation, only: init_rad, radtransfer
use Register
use Selfgravity, only: calc_selfpotential
use SharedVariables, only: sharedvars_clean_up
use Slices, only: setup_slices, wvid_prepare, wvid
use Snapshot
use Solid_Cells, only: init_solid_cells
use Special, only: init_special, initialize_mult_special
use Sub
use Syscalls, only: memusage
use Ascalar, only: init_acc
use Testfield, only: init_aatest
use Testflow, only: init_uutest
!
implicit none
!
real, allocatable, dimension (:,:,:,:) :: f, df
integer :: i, ifilter, stat
logical :: lnoerase=.false.
real :: dang
integer :: memory, memuse, memcpu
!
lstart = .true.
!
! Check processor layout, get processor numbers and define whether we are root.
!
call mpicomm_init
!
! Initialize OpenMP use
!
!$ include 'omp_init.h'
!
! Identify version.
!
if (lroot) call svn_id( &
'$Id$')
!
! Initialize the message subsystem, eg. color setting etc.
!
call initialize_messages
!
! Initialize use of multiple special modules if relevant.
!
call initialize_mult_special
!
! Allocate large arrays. We need to make them allocatable in order to
! avoid segfaults at 128^3 (7 variables) with Intel compiler on 32-bit
! Linux boxes. Not clear why they crashed (we _did_ increase stacksize
! limits), but they did and now don't. Also, the present approach runs
! up to nx=ny=nz=135, but not for even slightly larger grids.
!
if (.not.lnowrite) then
allocate(f(mx,my,mz,mfarray),STAT=stat)
if (stat>0) call fatal_error('start','Could not allocate f')
endif
!
! Pre-initialize f and df to absurd value (to crash the code should we
! later use uninitialized slots of those fields).
!
f=huge(1.0)
if (lmodify .or. nfilter/=0) then
allocate(df(mx,my,mz,mvar),STAT=stat)
if (stat>0) call fatal_error('start','Could not allocate df')
df=huge(1.0)
endif
!
! Read initialization parameters from "start.in".
!
call read_all_init_pars
!
call set_coorsys_dimmask
!
! Set up directory names and check whether the directories exist.
!
call directory_names
!
! Initialise MPI communication.
!
call initialize_mpicomm
!
! Initialise HDF5 communication.
!
call init_hdf5
call initialize_hdf5
!
! Register variables in the f array.
!
call register_modules
!allocate( f(mx,my,mz,nvar+naux+nscratch+nglobal),STAT=stat)
!allocate(df(mx,my,mz,nvar) ,STAT=stat)
if (lparticles) call particles_register_modules
!
! The logical headtt is sometimes referred to in start.x, even though it is
! not yet defined. So we set it simply to lroot here.
!
headtt=lroot
!
! Initialize start time, unless lread_oldsnap=T
!
if (.not.lread_oldsnap) t=tstart
!
! Will we write all slots of f?
!
if (lwrite_aux) then
mvar_io=mvar+maux
else
mvar_io=mvar
endif
!
! Print resolution.
!
if (lroot) print*, 'nxgrid, nygrid, nzgrid=', nxgrid, nygrid, nzgrid
!
! Unfortunately the following test for existence of directory fails under
! OSF1:
! inquire(FILE=trim(directory_snap), EXIST=exist)
! if (.not. exist) &
! call stop_it('Need directory <' // trim(directory_snap) // '>')
!
!
! Set box dimensions, make sure Lxyz and xyz1 are in sync.
! Box defaults to [-pi,pi] for all directions if none of xyz1 or Lxyz are set.
! If luniform_z_mesh_aspect_ratio=T, the default Lz scales with nzgrid/nxgrid.
! If wav1 is set, we initialize Lxyz for a cubic domain of size 2*pi/wav1.
! If, in addition to wav1, also wav1z is set, we initialize Lz separately.
!
if (wav1/=impossible) then
Lxyz=twopi/wav1
xyz0=-Lxyz/2.
if (wav1z/=impossible) then
Lxyz(3)=twopi/wav1z
xyz0(3)=-Lxyz(3)/2.
endif
endif
!
do i=1,3
if (Lxyz(i) == impossible) then
if (xyz1(i) == impossible) then
if (i==3.and.luniform_z_mesh_aspect_ratio) then
Lxyz(i)=2*pi*real(nzgrid)/real(nxgrid)
xyz0(i)=-pi*real(nzgrid)/real(nxgrid)
!
! FG: force theta coordinate to span 0:pi for periodic across pole
!
elseif (lpole(i)) then
if (lperi(i)) call fatal_error('start', &
'lperi and lpole cannot be used together in same component')
if (.not.lspherical_coords) call fatal_error('start', &
'lpole only implemented for spherical coordinates')
if (i==2) then
xyz0(i) = 0.
Lxyz(i) = pi
xyz0(3) = 0.
Lxyz(3) = 2.*pi
else
call fatal_error('start','origin/pole not included for components or coordinates')
endif
else
Lxyz(i)=2.*pi ! default value
endif
else
if (lpole(i)) then ! overwrite xyz0 and xyz1 to 0:pi
if (.not. lspherical_coords) call fatal_error('start', &
'lpole only implemented for spherical coordinates')
if (lperi(i)) call fatal_error('start', &
'lperi and lpole cannot be used together in same component')
if (i==2) then
xyz0(i) = 0.
Lxyz(i) = pi
xyz0(3) = 0.
xyz1(3) = 2.*pi
else
call fatal_error('start','origin/pole not included for components or coordinates')
endif
else
Lxyz(i)=xyz1(i)-xyz0(i)
endif
endif
else ! Lxyz was set
if (xyz1(i)/=impossible) & ! both Lxyz and xyz1 are set
call fatal_error('start','Cannot set Lxyz and xyz1 at the same time')
endif
!
! Possibility that mesh size was given in units of pi
!
if (xyz_units(i)=='pi') then
xyz0(i)=xyz0(i)*pi
Lxyz(i)=Lxyz(i)*pi
endif
enddo
!
if (lyinyang) then
if (lroot) print*, 'Setting latitude and longitude intervals for Yin-Yang grid, ignoring input'
!
! Min(dy,dz) put between Yin and Yang grid at closest distance to minimize overlap.
!
dang=rel_dang*min(1./max(1,(nygrid-1)),3./max(1,(nzgrid-1)))*0.5*pi
!dang=.99*min(1./max(1,(nygrid-1)),3./max(1,(nzgrid-1)))*0.5*pi ! only valid for equidistant grid!!
!dang=-1.8*min(1./max(1,(nygrid-1)),3./max(1,(nzgrid-1)))*0.5*pi ! only valid for equidistant grid!!
!dang=-2.8*min(1./max(1,(nygrid-1)),3./max(1,(nzgrid-1)))*0.5*pi ! only valid for equidistant grid!!
xyz0(2:3) = (/ 1./4., 1./4. /)*pi+0.5*dang
Lxyz(2:3) = (/ 1./2., 3./2. /)*pi-dang
endif
xyz1=xyz0+Lxyz
!
! Abbreviations
!
x0 = xyz0(1)
y0 = xyz0(2)
z0 = xyz0(3)
Lx = Lxyz(1)
Ly = Lxyz(2)
Lz = Lxyz(3)
!
! Position of equator (if any).
!
if (lequatory) yequator=xyz0(2)+0.5*Lxyz(2)
if (lequatorz) zequator=xyz0(3)+0.5*Lxyz(3)
!
! Check consistency.
!
if (.not.lperi(1).and.nxgrid<2) call fatal_error('start','for lperi(1)=F: must have nxgrid>1')
if (.not.lperi(2).and.nygrid<2) call fatal_error('start','for lperi(2)=F: must have nygrid>1')
if (.not.lperi(3).and.nzgrid<2) call fatal_error('start','for lperi(3)=F: must have nzgrid>1')
!
! Initialise random number generator in processor-independent fashion
!
call get_nseed(nseed) ! get state length of random number generator
call random_seed_wrapper(GET=seed,CHANNEL=1)
if (ichannel2>1) call random_seed_wrapper(GET=seed2,CHANNEL=2)
!
! Generate grid and initialize specific grid variables.
!
call construct_grid(x,y,z,dx,dy,dz)
!
! Write grid.dat file.
!
call wgrid('grid.dat',lwrite=.true.)
call wproc_bounds(trim(datadir) // "/proc_bounds.dat")
!if (lparticles.or.lpointmasses.or.lshear) call wproc_bounds(trim(directory)//'/proc_bounds.dat')
!
! Call rprint_list to initialize diagnostics and write indices to file.
!
call rprint_list(.false.)
if (lparticles) call particles_rprint_list(.false.)
!
if (lreport_undefined_diagnostics) call report_undefined_diagnostics
!
! Update the list of neighboring processes.
!
call update_neighbors
!
! Write .general file for data explorer.
!
if (lroot) call write_dx_general(trim(datadir)//'/var.general', &
x0-nghost*dx,y0-nghost*dy,z0-nghost*dz)
!
! Parameter dependent initialization of module variables and final
! pre-timestepping setup (must be done before need_XXXX can be used, for
! example).
!
call initialize_boundcond
call initialize_modules(f)
if (lparticles) call particles_initialize_modules(f)
!
! Initial conditions: by default, we put f=0 (ss=lnrho=uu=0, etc).
! alternatively: read existing snapshot and overwrite only some fields
! by the various procedures below.
!
! If lmodify is true, read var.dat into auxiliary array df
! and read modify_filename into f
!
if (lmodify) then
call rsnap('var.dat',df,mvar,lread_nogrid)
call rsnap(modify_filename,f,mvar,lread_nogrid)
else
if (lread_oldsnap) then
call rsnap('var.dat',f,mvar,lread_nogrid)
else
f(:,:,:,1:mvar)=0.0
endif
endif
!
if (lyinyang) call warning('start', &
'Most initial conditions depending on y or z will be set correctly only on Yin grid')
!
! Initialise random number generator in processor-dependent fashion for
! random initial data.
! Slightly tricky, since setting seed=(/iproc,0,0,0,0,0,0,0,.../)
! would produce very short-period random numbers with the Intel compiler;
! so we need to retain most of the initial entropy of the generator.
! Different initial seed (seed0) and random numbers on different CPUs
! The default is seed0=1812 for some obscure Napoleonic reason
! Fred: NB, when using persistent variables take care whether seeds need
! to be synchronous or not for init_*
! Maybe we want to allow the 2 channels to be initialized separately.
!
seed(1)=-((seed0-1812+1)*10+iproc_world)
if (ichannel1<2) call random_seed_wrapper(PUT=seed,CHANNEL=1)
if (ichannel2>1) call random_seed_wrapper(PUT=seed2,CHANNEL=2)
!
! Set random seed independent of processor prior to initial conditions.
! Do this only if seed0 is modified from its original value.
! NB Default is proc dependent seed during init_* then reverts proc independ
! seed0\=1812 proc independent seed throughout
!
if (seed0/=1812) then
seed(1)=seed0
if (ichannel1<2) call random_seed_wrapper(PUT=seed,CHANNEL=1)
if (ichannel2>1) call random_seed_wrapper(PUT=seed2,CHANNEL=2)
endif
!
! calculate scale factor of the universe
!
if (lread_scl_factor_file_new) call calc_scl_factor
!
! The following init routines only need to add to f.
! wd: also in the case where we have read in an existing snapshot??
!
do i=1,init_loops
if (lroot .and. init_loops/=1) &
print '(A33,i3,A25)', 'start: -- performing loop number', i, ' of initial conditions --'
call init_gg(f)
!
! This is a temporary solution for calculating the correct initial
! condition for anelastic case where we need to set div(rho u)=0
! The dust-vortex auto-test requires that uu is initialised before
! lnrho
! DM: it seems to me that the sequence of calls below is important.
! so this is just a caution : Please do not modify the sequence of
! calls to 'init' routines below.
!
if (lanelastic) then
call init_lnrho(f)
call init_uu(f)
else
call init_uu(f)
call init_lnrho(f)
endif
call init_energy(f)
call init_implicit_physics(f)
call init_aa(f)
call init_lorenz_gauge(f)
call init_poly(f)
call init_aatest(f)
call init_uutest(f)
call init_rad(f)
call init_lncc(f)
call init_acc(f)
call init_chiral(f)
call init_chemistry(f)
call init_uud(f)
call init_nd(f)
call init_uun(f)
call init_lnrhon(f)
call init_ecr(f)
call init_fcr(f)
call init_interstellar(f)
call init_solid_cells(f)
call init_pointmasses(f)
call init_special(f)
!
! If desired, the f array can be initialized in one call.
!
if (linitial_condition) call initial_condition_all(f)
call init_eos(f) ! has to come last
enddo
!
! Initialize particles.
!
if (lparticles) call particles_init(f)
!
if (lradiation_ray) then
call update_ghosts(f)
call radtransfer(f)
endif
!
! Filter initial velocity.
! NOTE: this procedure is currently not very efficient,
! because for all variables, boundary conditions and communication
! are done again and again for each variable.
!
if (nfilter/=0) then
do ifilter=1,nfilter
call rmwig(f,df,iux,iuz,awig)
enddo
if (lroot) print*,'DONE: filter initial velocity, nfilter=',nfilter
endif
!
! Calculate the potential of the self-gravity (mostly for debugging).
!
call calc_selfpotential(f)
!
! For sixth order momentum-conserving, symmetric hyperviscosity with positive
! definite heating rate we need to precalculate the viscosity term. The
! restivitity term for sixth order hyperresistivity with positive definite
! heating rate must also be precalculated.
!
if (lhyperviscosity_strict) call hyperviscosity_strict(f)
if (lhyperresistivity_strict) call hyperresistivity_strict(f)
!
! Set random seed independent of processor after initial conditions.
! Do this only if seed0 is not already changed from its original value.
! This seems to be a bug. If we want to set seed0 globally, why should it
! then always be 1812? If added "lseed_global.and." so we can turn it off.
!
if (lseed_global.and.seed0==1812) then
seed(1)=seed0
if (ichannel1<2) call random_seed_wrapper(PUT=seed,CHANNEL=1)
if (ichannel2>1) call random_seed_wrapper(PUT=seed2,CHANNEL=2)
elseif (lseed_procdependent) then
seed(1)=-((seed0-1812+1)*10+iproc_world)
if (ichannel1<2) call random_seed_wrapper(PUT=seed,CHANNEL=1)
if (ichannel2>1) call random_seed_wrapper(PUT=seed2,CHANNEL=2)
endif
!
! Final update of ghost cells, after that 'f' is only altered by smoothing if requested.
!
call update_ghosts(f)
if (lsmooth_farray .and. dimensionality>0 .and. (iux/=0.or.iax/=0)) then
!Why not all variables?
call smooth(f,merge(iux,iax,iux>0),merge(iaz,iuz,iaz>0),.true.,farray_smooth_width)
call update_ghosts(f)
endif
!
! Write initial condition to disk.
!
if (lwrite_ic) then
if (lparticles) call write_snapshot_particles(f,ENUM=.false.,snapnum=0)
!
call wsnap('VAR0',f,mvar_io,ENUM=.false.,FLIST='varN.list',noghost=.true.)
call pointmasses_write_snapshot('QVAR0',ENUM=.false.,FLIST='qvarN.list')
endif
!
! The option lnowrite writes everything except the actual var.dat file.
! This can be useful if auxiliary files are outdated, and don't want
! to overwrite an existing var.dat.
!
lnoerase = control_file_exists("NOERASE")
if ( (.not.lnowrite .and. .not.lnoerase) .or. lwrite_var_anyway) then
if (ip<12) print*,'START: writing to '//trim(directory_snap)//'/var.dat'
if (lparticles) call write_snapshot_particles(f,ENUM=.false.)
call pointmasses_write_snapshot('qvar.dat',ENUM=.false.)
call wsnap('var.dat',f,mvar_io,ENUM=.false.,noghost=.true.)
elseif (lmodify) then
call wsnap(modify_filename,f,mvar_io,ENUM=.false.,noghost=.true.)
endif
call wdim('dim.dat')
if (lroot) then
if (lparticles) call write_dim_particles(trim(datadir))
call pointmasses_write_qdim(trim(datadir)//'/qdim.dat')
endif
!
! Added here possibility two write slice file from start.
! This is possible for the radiation module, for example.
!
if (dvid/=0.) then
call setup_slices
call wvid_prepare
call wvid(f)
endif
!
! Write global variables.
!
if (mglobal/=0) call output_globals('global.dat',f(:,:,:,mvar+maux+1:mvar+maux+mglobal),mglobal)
!
! If requested, write original (z-dependent) stratification to file.
!
call write_z_stratification(f)
!
! Write input parameters to a parameter file (for run.x and IDL).
! Do this late enough, so initialization routines can adjust them.
!
call write_all_init_pars('IDL')
!
! Write information about pencils to disc.
!
call write_pencil_info
memuse=memusage()
call mpireduce_max_int(memuse,memcpu)
call mpireduce_sum_int(memuse,memory)
if (lroot) then
print'(1x,a,f9.3)', 'Maximum used memory per cpu [MBytes] = ', memcpu/1024.
if (memory>1e6) then
print'(1x,a,f12.3)', 'Maximum used memory [GBytes] = ', memory/1024.**2
else
print'(1x,a,f12.3)', 'Maximum used memory [MBytes] = ', memory/1024.
endif
print*
endif
!
! Gvie all modules the possibility to exit properly.
!
call finalize_modules(f)
!
! Stop MPI.
!
call mpifinalize
!
! Free any allocated memory.
! MR: do we need all these cleanups? The process ends immediately afterwards and all memory is released anyway.
!
call farray_clean_up
call sharedvars_clean_up
call initial_condition_clean_up
if (lparticles) call particles_cleanup
!
! Deallocate the arrays allocated for 1-D and 2-D diagnostics.
!
call diagnostics_clean_up
!
! announce completion.
!
if (lroot) print*
if (lroot) print*,'start.x has completed successfully'
if (lroot) print*
!
endprogram