! $Id: equ.f90 13579 2010-03-31 11:58:14Z AxelBrandenburg $
!
! This module adds evolution terms to the dynamical equations, calling
! subroutines in the chosen set of physics modules.
!
module Equ
!
use Cdata
use Messages
!
implicit none
!
private
!
public :: pde, debug_imn_arrays, initialize_pencils
!
contains
!***********************************************************************
include 'pencil_init.inc' ! defines subroutine initialize_pencils()
!***********************************************************************
subroutine pde(f,df,p)
!
! Call the different evolution equations.
!
! 10-sep-01/axel: coded
!
use Boundcond
use Density
use Diagnostics
use Entropy
use EquationOfState
use Gravity
use Grid, only: calc_pencils_grid
use Hydro
use Magnetic
use Mpicomm
use Shear
use Shock, only: calc_pencils_shock, calc_shock_profile, &
calc_shock_profile_simple
use Sub
use Viscosity, only: calc_viscosity, calc_pencils_viscosity
!
logical :: early_finalize
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
real, dimension (nx) :: maxadvec,advec2,maxdiffus,maxdiffus3
intent(inout) :: f ! inout due to lshift_datacube_x,
! density floor, or velocity ceiling
intent(out) :: df, p
!
! Print statements when they are first executed.
!
headtt = headt .and. lfirst .and. lroot
!
if (headtt.or.ldebug) print*,'pde: ENTER'
if (headtt) call svn_id( &
"$Id: equ.f90 13579 2010-03-31 11:58:14Z AxelBrandenburg $")
!
! Initialize counter for calculating and communicating print results.
! Do diagnostics only in the first of the 3 (=itorder) substeps.
!
ldiagnos =lfirst.and.lout
l1davgfirst=lfirst.and.l1davg
l2davgfirst=lfirst.and.l2davg
!
! Record times for diagnostic and 2d average output.
!
if (ldiagnos) tdiagnos=t ! (diagnostics are for THIS time)
if (l1davgfirst) t1ddiagnos=t ! (1-D averages are for THIS time)
if (l2davgfirst) t2davgfirst=t ! (2-D averages are for THIS time)
!
! Grid spacing.
!
dxyz_2 = dx_1(l1:l2)**2+dy_1(m1)**2+dz_1(n1)**2
dxyz_4 = dx_1(l1:l2)**4+dy_1(m1)**4+dz_1(n1)**4
dxyz_6 = dx_1(l1:l2)**6+dy_1(m1)**6+dz_1(n1)**6
!
! Write crash snapshots to the hard disc if the time-step is very low.
! The user must have set crash_file_dtmin_factor>0.0 in &run_pars for
! this to be done.
!
if (crash_file_dtmin_factor > 0.0) call output_crash_files(f)
!
! Call "before_boundary" hooks (for f array precalculation)
!
if (ldensity) call density_before_boundary(f)
if (lshear) call shear_before_boundary(f)
!
! Initiate shock profile calculation and use asynchronous to handle
! communication along processor/periodic boundaries.
!
if (lshock) call calc_shock_profile(f)
!
! Prepare x-ghost zones; required before f-array communication
! AND shock calculation
!
call boundconds_x(f)
!
! Initiate (non-blocking) communication and do boundary conditions.
! Required order:
! 1. x-boundaries (x-ghost zones will be communicated) - done above
! 2. communication
! 3. y- and z-boundaries
!
if (ldebug) print*,'pde: bef. initiate_isendrcv_bdry'
call initiate_isendrcv_bdry(f)
if (early_finalize) then
call finalize_isendrcv_bdry(f)
call boundconds_y(f)
call boundconds_z(f)
endif
!
! Set inverse timestep to zero before entering loop over m and n.
!
if (lfirst.and.ldt) then
if (dtmax/=0.0) then
dt1_max=1/dtmax
else
dt1_max=0.0
endif
endif
!
! Calculate shock profile (simple).
!
if (lshock) call calc_shock_profile_simple(f)
!
! Calculate averages, currently only required for certain settings
! in hydro of the testfield procedure (only when lsoca=.false.)
!
call timing('pde','before mn-loop')
!
!------------------------------------------------------------------------------
! Do loop over m and n.
!
mn_loop: do imn=1,ny*nz
n=nn(imn)
m=mm(imn)
lfirstpoint=(imn==1) ! true for very first m-n loop
llastpoint=(imn==(ny*nz)) ! true for very last m-n loop
!
! Store the velocity part of df array in a temporary array
! while solving the anelastic case.
!
call timing('pde','before ldensity_anelastic',mnloop=.true.)
!
! Make sure all ghost points are set.
!
if (.not.early_finalize.and.necessary(imn)) then
call finalize_isendrcv_bdry(f)
call boundconds_y(f)
call boundconds_z(f)
endif
call timing('pde','finished boundconds_z',mnloop=.true.)
!
! For each pencil, accumulate through the different modules
! advec_XX and diffus_XX, which are essentially the inverse
! advective and diffusive timestep for that module.
! (note: advec_cs2 and advec_va2 are inverse _squared_ timesteps)
!
if (lfirst.and.ldt.and.(.not.ldt_paronly)) then
if (lhydro) then
advec_uu=0.0
endif
if (ldensity) then
diffus_diffrho=0.0; advec_lnrho=0.0
if (leos) advec_cs2=0.0
endif
if (lentropy) then
diffus_chi=0.0
endif
if (lmagnetic) then
advec_va2=0.0
diffus_eta=0.0
endif
if (lviscosity) then
diffus_nu=0.0
endif
if (lshear) then
advec_shear=0.0
endif
endif
!
! Grid spacing.
!
dxyz_2 = dx_1(l1:l2)**2+dy_1(m1)**2+dz_1(n1)**2
dxyz_4 = dx_1(l1:l2)**4+dy_1(m1)**4+dz_1(n1)**4
dxyz_6 = dx_1(l1:l2)**6+dy_1(m1)**6+dz_1(n1)**6
!
! Calculate grid/geometry related pencils.
!
call calc_pencils_grid(f,p)
!
! Calculate pencils for the pencil_case.
! Note: some no-modules (e.g. nohydro) also calculate some pencils,
! so it would be wrong to check for lhydro etc in such cases.
!
call calc_pencils_hydro(f,p)
call calc_pencils_density(f,p)
call calc_pencils_eos(f,p)
if (lshock) call calc_pencils_shock(f,p)
if (lviscosity) call calc_pencils_viscosity(f,p)
if (lmagnetic) call calc_pencils_magnetic(f,p)
if (lgrav) call calc_pencils_gravity(f,p)
if (lshear) call calc_pencils_shear(f,p)
!
! --------------------------------------------------------
! NO CALLS MODIFYING PENCIL_CASE PENCILS BEYOND THIS POINT
! --------------------------------------------------------
!
! hydro, density, entropy, and magnetic evolution
! Note that pressure gradient is added in dss_dt to momentum,
! even if lentropy=.false.
!
call duu_dt(df,p)
call dlnrho_dt(df,p)
call dss_dt(df,p)
call daa_dt(df,p)
!
! Add gravity, if present
!
if (lgrav) then
if (lhydro) call duu_dt_grav(f,df,p)
endif
!
! Add shear if present
!
if (lshear) call shearing(f,df,p)
!
! In max_mn maximum values of u^2 (etc) are determined sucessively
! va2 is set in magnetic (or nomagnetic)
! In rms_mn sum of all u^2 (etc) is accumulated
! Calculate maximum advection speed for timestep; needs to be done at
! the first substep of each time step
! Note that we are (currently) accumulating the maximum value,
! not the maximum squared!
!
! The dimension of the run ndim (=0, 1, 2, or 3) enters the viscous time step.
! This has to do with the term on the diagonal, cdtv depends on order of scheme
!
if (lfirst.and.ldt.and.(.not.ldt_paronly)) then
!
! sum or maximum of the advection terms?
! (lmaxadvec_sum=.false. by default)
!
maxadvec=0.0
maxdiffus=0.0
if (lhydro) maxadvec=maxadvec+advec_uu
if (lshear) maxadvec=maxadvec+advec_shear
if (ldensity.or.lmagnetic) then
advec2=0.0
if (ldensity) advec2=advec2+advec_cs2+advec_lnrho**2
if (lmagnetic) advec2=advec2+advec_va2
maxadvec=maxadvec+sqrt(advec2)
endif
!
! Time step constraints from each module.
! (At the moment, magnetic and testfield use the same variable.)
!
if (lviscosity) maxdiffus=max(diffus_nu,maxdiffus)
if (ldensity) maxdiffus=max(diffus_diffrho,maxdiffus)
if (lentropy) maxdiffus=max(diffus_chi,maxdiffus)
if (lmagnetic) maxdiffus=max(diffus_eta,maxdiffus)
!
! cdt, cdtv, and cdtc are empirical non-dimensional coefficients
!
dt1_advec = maxadvec/cdt
dt1_diffus = maxdiffus/cdtv
dt1_max = max(dt1_max,sqrt(dt1_advec**2+dt1_diffus**2))
!
! Diagnostics showing how close to advective and diffusive time steps we are
!
if (ldiagnos.and.idiag_dtv/=0) then
call max_mn_name(maxadvec/cdt,idiag_dtv,l_dt=.true.)
endif
if (ldiagnos.and.idiag_dtdiffus/=0) then
call max_mn_name(maxdiffus/cdtv,idiag_dtdiffus,l_dt=.true.)
endif
!
! Regular and hyperdiffusive mesh Reynolds numbers
!
if (ldiagnos) then
if (idiag_Rmesh/=0) &
call max_mn_name(pi_1*maxadvec/(maxdiffus+tini),idiag_Rmesh)
if (idiag_Rmesh3/=0) &
call max_mn_name(pi5_1*maxadvec/(maxdiffus3+tini),idiag_Rmesh3)
if (idiag_maxadvec/=0) &
call max_mn_name(maxadvec,idiag_maxadvec)
endif
endif
!
call timing('pde','end of mn loop',mnloop=.true.)
!
! End of loops over m and n.
!
headtt=.false.
enddo mn_loop
call timing('pde','at the end of the mn_loop')
!
! ----------------------------------------
! NO CALLS MODIFYING DF BEYOND THIS POINT
! ----------------------------------------
!
! Check for NaNs in the advection time-step.
!
if (notanumber(dt1_advec)) then
print*, 'pde: dt1_advec contains a NaN at iproc=', iproc
if (lhydro) print*, 'advec_uu =',advec_uu
if (lshear) print*, 'advec_shear=',advec_shear
if (lentropy) print*, 'advec_cs2 =',advec_cs2
if (lmagnetic) print*, 'advec_va2 =',advec_va2
call fatal_error_local('pde','')
endif
!
! Diagnostics.
!
if (ldiagnos) call diagnostic
!
! 1-D diagnostics.
!
if (l1davgfirst) then
call xyaverages_z
call xzaverages_y
call yzaverages_x
endif
!
! 2-D averages.
!
if (l2davgfirst) then
if (lwrite_yaverages) call yaverages_xz
if (lwrite_zaverages) call zaverages_xy
endif
!
! Note: zaverages_xy are also needed if bmx and bmy are to be calculated
! (of course, yaverages_xz does not need to be calculated for that).
!
if (.not.l2davgfirst.and.ldiagnos.and.ldiagnos_need_zaverages) then
if (lwrite_zaverages) call zaverages_xy
endif
!
! Reset lwrite_prof.
!
lwrite_prof=.false.
!
endsubroutine pde
!***********************************************************************
subroutine debug_imn_arrays
!
! For debug purposes: writes out the mm, nn, and necessary arrays.
!
! 23-nov-02/axel: coded
!
open(1,file=trim(directory)//'/imn_arrays.dat')
do imn=1,ny*nz
if (necessary(imn)) write(1,'(a)') '----necessary=.true.----'
write(1,'(4i6)') imn,mm(imn),nn(imn)
enddo
close(1)
!
endsubroutine debug_imn_arrays
!***********************************************************************
subroutine output_crash_files(f)
!
! Write crash snapshots when time-step is low.
!
! 15-aug-2007/anders: coded
!
use Snapshot
!
real, dimension(mx,my,mz,mfarray) :: f
!
integer, save :: icrash=0
character (len=10) :: filename
character (len=1) :: icrash_string
!
if ( (it>1) .and. (itsub==1) .and. (dt<=crash_file_dtmin_factor*dtmin) ) then
write(icrash_string, fmt='(i1)') icrash
filename='crash'//icrash_string//'.dat'
call wsnap(trim(directory_snap)//'/'//filename,f,mvar_io,.false.)
if (lroot) then
print*, 'Time-step is very low - writing '//trim(filename)
print*, '(it, itsub=', it, itsub, ')'
print*, '(t, dt=', t, dt, ')'
endif
!
! Next crash index, cycling from 0-9 to avoid excessive writing of
! snapshots to the hard disc.
!
icrash=icrash+1
icrash=mod(icrash,10)
endif
!
endsubroutine output_crash_files
!***********************************************************************
endmodule Equ