! $Id: mhs_equilibrium.f90 10874 2009-05-17 16:34:17Z wdobler $
!
! Initial condition (density, magnetic field, velocity)
! for magnetohydrostatical equilibrium in a global accretion
! disk with an imposed (cylindrically symmetric) sound speed
! profile in spherical coordinates.
!
! 07-may-09/wlad: adapted from noinitial_condition.f90
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: linitial_condition = .true.
!
!***************************************************************
module InitialCondition
!
use Cdata
use Cparam
use Messages
use Sub, only: keep_compiler_quiet
!
implicit none
!
include '../initial_condition.h'
!
real :: density_power_law,temperature_power_law,plasma_beta
logical :: lnumerical_mhsequilibrium=.true.
!
namelist /initial_condition_pars/ &
density_power_law,temperature_power_law,plasma_beta,&
lnumerical_mhsequilibrium
!
contains
!***********************************************************************
subroutine register_initial_condition()
!
! Configure pre-initialised (i.e. before parameter read) variables
! which should be know to be able to evaluate
!
! 07-oct-09/wlad: coded
!
! Identify CVS/SVN version information.
!
if (lroot) call svn_id( &
"$Id: mhs_equilibrium.f90 10874 2009-05-17 16:34:17Z wdobler $")
!
endsubroutine register_initial_condition
!***********************************************************************
subroutine initial_condition_uu(f)
!
! Initialize the velocity field.
!
! 07-may-09/wlad: coded
!
use Gravity, only: acceleration
use Sub, only: get_radial_distance
use FArrayManager, only: farray_use_global
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
real, dimension (mx) :: rr_sph,rr_cyl,g_r
real, dimension (mx) :: OOK2,OO2,tmp1,H2
real :: p,q,tmp2,ksi
integer, pointer :: iglobal_cs2
integer :: ics2
!
! Set the sound speed
!
call set_sound_speed(f)
!
! Get the sound speed
!
nullify(iglobal_cs2)
call farray_use_global('cs2',iglobal_cs2)
ics2=iglobal_cs2
!
! Analytical expression that leads to an equilibrium configuration.
! Commented out because the numerical equilibrium enforced below
! is a lot better when it comes to reproduce what the code actually
! solves.
!
if (.not.lnumerical_mhsequilibrium) then
!
if (lmagnetic) then
ksi=(1.+plasma_beta)/plasma_beta
else
ksi=1.
endif
!
p=-density_power_law
q=-temperature_power_law
!
do m=1,my;do n=1,mz
call acceleration(g_r)
call get_radial_distance(rr_sph,rr_cyl)
OOK2=max(-g_r/(rr_sph*sinth(m)**3),0.)
!
H2=f(:,m,n,ics2)/OOK2
!
tmp1=H2/rr_cyl**2*(ksi*(p+q-2.) + 2.)
tmp2=ksi*(q+1.)*(1.-sinth(m))
!
OO2=OOK2*(tmp1+tmp2+sinth(m))
!
f(:,m,n,iuz) = f(:,m,n,iuz) + rr_cyl*sqrt(OO2)
enddo;enddo
!
endif
!
endsubroutine initial_condition_uu
!***********************************************************************
subroutine initial_condition_lnrho(f)
!
! Initialize logarithmic density. init_lnrho
! will take care of converting it to linear
! density if you use ldensity_nolog
!
! 07-may-09/wlad: coded
!
use EquationOfState, only: rho0
use FArrayManager, only: farray_use_global
use Gravity, only: potential
use Sub, only: get_radial_distance
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
real, dimension (mx) :: rr_sph,rr_cyl
real, dimension (mx) :: lnrhomid,strat
real, dimension (mx) :: cs2,tmp1,tmp2
real :: p,q
integer :: ics2
integer, pointer :: iglobal_cs2,iglobal_glnTT
!
p=-density_power_law
q=-temperature_power_law
!
if (lroot) print*,&
'initial_condition_lnrho: locally isothermal approximation'
if (lroot) print*,'Radial density stratification with power law=',p
if (lroot) print*,'Radial temperature stratification with power law=',q
!
! Get the sound speed globals
!
nullify(iglobal_cs2)
call farray_use_global('cs2',iglobal_cs2)
ics2=iglobal_cs2
!
! Pencilize the density allocation.
!
do n=1,mz
do m=1,my
!
! Midplane density
!
call get_radial_distance(rr_sph,rr_cyl)
lnrhomid=log(rho0)+p*log(rr_cyl)
!
! Vertical stratification
!
cs2=f(:,m,n,ics2)
call potential(POT=tmp1,RMN=rr_sph)
call potential(POT=tmp2,RMN=rr_cyl)
strat=-(tmp1-tmp2)/cs2
f(:,m,n,ilnrho) = f(:,m,n,ilnrho)+lnrhomid+strat
!
enddo
enddo
!
! If the run is non-magnetic, this is the last routine called.
! Enforce numerical equilibrium then
!
if ((.not.lmagnetic).and.lnumerical_mhsequilibrium) &
call enforce_numerical_equilibrium(f,lhd=.true.)
!
endsubroutine initial_condition_lnrho
!***********************************************************************
subroutine initial_condition_aa(f)
!
! Initialize the magnetic vector potential. Constant plasma
! beta magnetic field.
!
! 07-may-09/wlad: coded
!
use FArrayManager, only: farray_use_global
use Sub, only: gij,curl_mn,get_radial_distance
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
real, dimension(mx) :: pressure,Bphi,Atheta
real, dimension(mx) :: rr_sph,rr_cyl,tmp,tmp2
real :: integral,dr
integer, pointer :: iglobal_cs2
integer :: ics2,irho,i
!
! Get the sound speed globals
!
nullify(iglobal_cs2)
call farray_use_global('cs2',iglobal_cs2)
ics2=iglobal_cs2
!
! Density is already in linear after init_lnrho, so
!
irho=ilnrho
!
do n=1,mz
do m=1,my
!
pressure=f(:,m,n,irho)*f(:,m,n,ics2)
!
! The following line assumes mu0=1
!
Bphi = sqrt(2*pressure/plasma_beta)
!
! Bphi = 1/r*d/dr(r*Atheta), so integrate: Atheta=1/r*Int(B*r)dr
!
call get_radial_distance(rr_sph,rr_cyl)
dr=rr_sph(2)-rr_sph(1)
tmp=Bphi*rr_sph
!
tmp2(1)=0.
do i=2,mx
tmp2(i)=tmp2(i-1) + tmp(i)*dr
enddo
Atheta=tmp2/rr_sph
!
f(:,m,n,iay)=f(:,m,n,iay)+Atheta
!
enddo
enddo
!
! All quantities are set. Enforce numerical equilibrium.
!
if (lnumerical_mhsequilibrium) &
call enforce_numerical_equilibrium(f,lhd=.false.)
!
endsubroutine initial_condition_aa
!***********************************************************************
subroutine set_sound_speed(f)
!
! Set the thermo-related quantities. Illustrates that
! the user can define as many internal routines as wanted.
!
! 10-may-09/wlad : moved from initial_condition_lnrho
!
use FArrayManager, only: farray_use_global
use EquationOfState, only: cs20
use Sub, only: get_radial_distance
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx) :: rr_sph,rr_cyl,cs2
integer, pointer :: iglobal_cs2,iglobal_glnTT
real :: q
integer :: ics2,iglnTT
!
q=-temperature_power_law
!
! Get the globals needed to store sound speed and temperature gradient
!
nullify(iglobal_cs2)
call farray_use_global('cs2',iglobal_cs2);ics2=iglobal_cs2
nullify(iglobal_glnTT)
call farray_use_global('glnTT',iglobal_glnTT);iglnTT=iglobal_glnTT
!
! Set the sound speed - a power law in cylindrical radius.
!
do m=1,my
do n=1,mz
call get_radial_distance(rr_sph,rr_cyl)
cs2=cs20*rr_cyl**q
!
! Store cs2 in one of the free slots of the f-array
!
f(:,m,n,ics2)=cs2
!
! Same for the temperature gradient
!
f(:,m,n,iglnTT )=q/rr_sph
f(:,m,n,iglnTT+1)=q/rr_sph*cotth(m)
f(:,m,n,iglnTT+2)=0.
!
enddo
enddo
!
endsubroutine set_sound_speed
!***********************************************************************
subroutine enforce_numerical_equilibrium(f,lhd)
!
! This subroutine does exactly what's done in runtime to the
! velocity field, in order to ensure precise numerical
! magnetohydrostatical equilibrium.
!
use Sub, only: grad,get_radial_distance,&
gij,curl_mn,gij_etc,cross_mn,multsv_mn
use FArrayManager, only: farray_use_global
! use BoundCond, only: update_ghosts
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,3,3) :: aij,bij
real, dimension (nx,3) :: glnTT,grho,glnrho
real, dimension (nx,3) :: aa,bb,jj,jxb,jxbr
real, dimension (nx) :: cs2,rr_sph,rr_cyl,rho1
real, dimension (nx) :: fpres_thermal,fpres_magnetic
integer, pointer :: iglobal_cs2, iglobal_glnTT
integer :: j,irho,ics2,iglnTT,i
logical :: lhd
!
! Get the temperature globals
!
nullify(iglobal_cs2)
call farray_use_global('cs2',iglobal_cs2);ics2=iglobal_cs2
nullify(iglobal_glnTT)
call farray_use_global('glnTT',iglobal_glnTT);iglnTT=iglobal_glnTT
!
! Use ilnrho as rho - works only for ldensity_nolog (which isn't passed
! yet, but, hey, it's MY own custom initial condition file. I know what
! it does. :-)
!
if (lhd) then
irho=iglnTT+2 !take an empty slot of f to put irho
f(l1:l2,m1:m2,n1:n2,irho)=exp(f(l1:l2,m1:m2,n1:n2,ilnrho))
else
irho=ilnrho
endif
!
!call update_ghosts(f)
!
! Azimuthal speed that perfectly balances the pressure gradient.
!
do m=m1,m2
do n=n1,n2
!
call grad(f,irho,grho)
do j=1,3
glnrho(:,j)=grho(:,j)/f(l1:l2,m,n,irho)
enddo
cs2=f(l1:l2,m,n,ics2)
glnTT(:,1:2)=f(l1:l2,m,n,iglnTT:iglnTT+1)
!
fpres_thermal=cs2*(glnrho(:,2)+glnTT(:,2))
!
if (lmagnetic) then
aa=f(l1:l2,m,n,iax:iaz) !aa
call gij(f,iaa,aij,1) !aij
call curl_mn(aij,bb,aa) !bb
call gij_etc(f,iaa,aa,aij,bij) !bij
call curl_mn(bij,jj,bb) !jj
call cross_mn(jj,bb,jxb) !jxb
rho1=1./f(l1:l2,m,n,irho) !1/rho
call multsv_mn(rho1,jxb,jxbr) !jxb/rho
fpres_magnetic=-jxbr(:,2)
else
fpres_magnetic=0.
endif
!
call get_radial_distance(rr_sph,rr_cyl)
f(l1:l2,m,n,iuz)=f(l1:l2,m,n,iuz)+&
sqrt(rr_sph*(fpres_thermal+fpres_magnetic)/cotth(m))
!
enddo
enddo
!
! Revert the free slot used for irho to its original value.
!
if (lhd) f(:,:,:,iglnTT+2)=0.
!
endsubroutine enforce_numerical_equilibrium
!***********************************************************************
subroutine read_initial_condition_pars(unit,iostat)
!
! 07-may-09/wlad: coded
!
integer, intent(in) :: unit
integer, intent(inout), optional :: iostat
!
if (present(iostat)) then
read(unit,NML=initial_condition_pars,ERR=99, IOSTAT=iostat)
else
read(unit,NML=initial_condition_pars,ERR=99)
endif
!
99 return
!
endsubroutine read_initial_condition_pars
!***********************************************************************
subroutine write_initial_condition_pars(unit)
!
integer, intent(in) :: unit
!
write(unit,NML=initial_condition_pars)
!
endsubroutine write_initial_condition_pars
!********************************************************************
!************ DO NOT DELETE THE FOLLOWING **************
!********************************************************************
!** This is an automatically generated include file that creates **
!** copies dummy routines from noinitial_condition.f90 for any **
!** InitialCondition routines not implemented in this file **
!** **
include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition