! $Id$
!
! This modules solves two reactive scalar advection equations
! This is used for modeling the spatial evolution of left and
! right handed aminoacids.
!
!** 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 :: lchiral = .true.
!
! MVAR CONTRIBUTION 2
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Chiral
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
!
implicit none
!
include 'chiral.h'
!
integer :: iXX_chiral=0, iYY_chiral=0, iZZ_chiral=0
integer :: iXX2_chiral=0, iYY2_chiral=0
character (len=labellen) :: initXX_chiral='zero', initYY_chiral='zero', initZZ_chiral='zero'
logical :: llorentzforceEP=.false., lZZ_chiral=.false.
logical :: linitialize_aa_from_EP=.false.
logical :: linitialize_aa_from_EP_alpgrad=.false.
logical :: linitialize_aa_from_EP_betgrad=.false.
real :: tinitialize_aa_from_EP=0.
real :: amplXX_chiral=.1, widthXX_chiral=.5
real :: amplYY_chiral=.1, widthYY_chiral=.5
real :: amplZZ_chiral=.1, widthZZ_chiral=.5
real :: kx_XX_chiral=1.,ky_XX_chiral=1.,kz_XX_chiral=1.,radiusXX_chiral=0.
real :: kx_YY_chiral=1.,ky_YY_chiral=1.,kz_YY_chiral=1.,radiusYY_chiral=0.
real :: xposXX_chiral=0.,yposXX_chiral=0.,zposXX_chiral=0.
real :: xposYY_chiral=0.,yposYY_chiral=0.,zposYY_chiral=0.
real :: initpowerYY_chiral=2., cutoffYY_chiral=0.
!
namelist /chiral_init_pars/ &
initXX_chiral, amplXX_chiral, kx_XX_chiral, ky_XX_chiral, kz_XX_chiral, &
initYY_chiral, amplYY_chiral, kx_YY_chiral, ky_YY_chiral, kz_YY_chiral, &
initZZ_chiral, amplZZ_chiral, lZZ_chiral, &
radiusXX_chiral, widthXX_chiral, &
radiusYY_chiral, widthYY_chiral, &
xposXX_chiral, yposXX_chiral, zposXX_chiral, &
xposYY_chiral, yposYY_chiral, zposYY_chiral, &
initpowerYY_chiral, cutoffYY_chiral
!
real :: chiral_diffXX=impossible, chiral_diff=0., chiral_crossinhibition=1.,chiral_fidelity=1.
real :: chiral_diffZZ=impossible
real :: chiral_fishernu=0., chiral_fisherK=1. !fishers equation growth rate, carrying capac.
real :: chiral_fishermu=0., chiral_fisherH=1. !fishers equation growth rate, carrying capac.
real :: chiral_fisherR=0., chiral_fisherR2=0. !(reinfections, second model corresponds to SIRS)
real :: chiral_fisherR2_tend=0., chiral_fisherR2_tstart=0. !(second model, tend and start time)
real, dimension(3) :: gradX0=(/0.0,0.0,0.0/), gradY0=(/0.0,0.0,0.0/)
character (len=labellen) :: chiral_reaction='BAHN_model'
logical :: limposed_gradient=.false.
logical :: lupw_chiral=.false.
!
namelist /chiral_run_pars/ &
chiral_diffXX, chiral_diff, chiral_diffZZ, chiral_crossinhibition, chiral_fidelity, &
chiral_reaction, limposed_gradient, gradX0, gradY0, &
chiral_fishernu, chiral_fisherK, chiral_fisherR, chiral_fisherR2, &
chiral_fishermu, chiral_fisherH, &
lupw_chiral, llorentzforceEP, &
linitialize_aa_from_EP,tinitialize_aa_from_EP, &
linitialize_aa_from_EP_alpgrad,linitialize_aa_from_EP_betgrad, &
chiral_fisherR2_tend, chiral_fisherR2_tstart
!
integer :: idiag_XX_chiralmax=0, idiag_XX_chiralm=0
integer :: idiag_YY_chiralmax=0, idiag_YY_chiralm=0
integer :: idiag_ZZ_chiralmax=0, idiag_ZZ_chiralm=0
integer :: idiag_QQm_chiral=0, idiag_QQ21m_chiral=0, idiag_QQ21QQm_chiral=0
integer :: idiag_brmsEP=0, idiag_bmaxEP=0
integer :: idiag_jrmsEP=0, idiag_jmaxEP=0
integer :: idiag_jbmEP=0
integer :: idiag_R2tdep=0
!
! Auxiliary variables.
!
real, dimension (nx,3) :: gXX_chiral, gYY_chiral, gZZ_chiral
real, dimension (nx) :: del2XX_chiral, del2YY_chiral,del2ZZ_chiral
real :: chiral_fisherR2_tdep
!
contains
!***********************************************************************
subroutine register_chiral()
!
! Initialise variables which should know that we solve for passive
! scalar: iXX_chiral and iYY_chiral; increase nvar accordingly
!
! 28-may-04/axel: adapted from pscalar
!
use FArrayManager
!
call farray_register_pde('XX_chiral',iXX_chiral)
call farray_register_pde('YY_chiral',iYY_chiral)
if (lZZ_chiral) call farray_register_pde('ZZ_chiral',iZZ_chiral)
!
! extra fields for Higgs
!
if (initYY_chiral=='power_randomphase_higgs' .or. &
initYY_chiral=='higgs_mono') then
call farray_register_pde('XX2_chiral',iXX2_chiral)
call farray_register_pde('YY2_chiral',iYY2_chiral)
endif
!
! Identify version number.
!
if (lroot) call svn_id( &
"$Id$")
!
endsubroutine register_chiral
!***********************************************************************
subroutine initialize_chiral(f)
!
! Perform any necessary post-parameter read initialization
! Dummy routine
!
! 28-may-04/axel: adapted from pscalar
!
real, dimension (mx,my,mz,mfarray) :: f
!
! default: chiral_diffXX=chiral_diff
!
if (chiral_diffXX==impossible) chiral_diffXX=chiral_diff
if (chiral_diffZZ==impossible) chiral_diffZZ=chiral_diff
if (headt) print*,'chiral_diffXX, chiral_diff, chiral_diffZZ =', &
chiral_diffXX, chiral_diff, chiral_diffZZ
!
! set f to zero and then call the same initial condition
! that was used in start.csh
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_chiral
!***********************************************************************
subroutine init_chiral(f)
!
! initialise passive scalar field; called from start.f90
!
! 28-may-04/axel: adapted from pscalar
!
use Sub
use General
use Initcond
use InitialCondition, only: initial_condition_chiral
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz) :: r, pom
real, dimension (mx,my,mz) :: norm1
!
! check first for initXX_chiral
!
select case (initXX_chiral)
case ('zero'); f(:,:,:,iXX_chiral)=0.
case ('const'); f(:,:,:,iXX_chiral)=amplXX_chiral
case ('random_minus1_to1'); call random_number_wrapper(r); f(:,:,:,iXX_chiral)=2.*r-1.
case ('blob'); call blob(amplXX_chiral,f,iXX_chiral,radiusXX_chiral,xposXX_chiral,yposXX_chiral,zposXX_chiral)
case ('hat-x'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,kx=kx_XX_chiral)
case ('hat-y'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,ky=ky_XX_chiral)
case ('hat-z'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,kz=kz_XX_chiral)
case ('gaussian-x'); call gaussian(amplXX_chiral,f,iXX_chiral,kx=kx_XX_chiral)
case ('gaussian-y'); call gaussian(amplXX_chiral,f,iXX_chiral,ky=ky_XX_chiral)
case ('gaussian-z'); call gaussian(amplXX_chiral,f,iXX_chiral,kz=kz_XX_chiral)
case ('gaussian-noise'); call gaunoise(amplXX_chiral,f,iXX_chiral)
case ('positive-noise'); call posnoise(amplXX_chiral,f,iXX_chiral)
case ('wave-x'); call wave(amplXX_chiral,f,iXX_chiral,kx=kx_XX_chiral)
case ('wave-y'); call wave(amplXX_chiral,f,iXX_chiral,ky=ky_XX_chiral)
case ('wave-z'); call wave(amplXX_chiral,f,iXX_chiral,kz=kz_XX_chiral)
case ('cos(x-cosz)'); call cosxz_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,kz_XX_chiral)
case ('cosy_sinz'); call cosy_sinz(amplXX_chiral,f,iXX_chiral,ky_XX_chiral,kz_XX_chiral)
case ('cosx_cosz'); call cosx_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,kz_XX_chiral)
case ('cosx_cosy_cosz'); call cosx_cosy_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,ky_XX_chiral,kz_XX_chiral)
case ('cosx_siny_cosz'); call cosx_siny_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,ky_XX_chiral,kz_XX_chiral)
case default; call fatal_error('init_chiral','no such initXX_chiral: '//trim(initXX_chiral))
endselect
!
! check next for initYY_chiral
!
select case (initYY_chiral)
case ('zero'); f(:,:,:,iYY_chiral)=0.
case ('const'); f(:,:,:,iYY_chiral)=amplYY_chiral
case ('random_0to_pi'); call random_number_wrapper(r); f(:,:,:,iYY_chiral)=2.*pi*r
case ('blob'); call blob(amplYY_chiral,f,iYY_chiral,radiusYY_chiral,xposYY_chiral,yposYY_chiral,zposYY_chiral)
case ('hat-x'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,kx=kx_YY_chiral)
case ('hat-y'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,ky=ky_YY_chiral)
case ('hat-z'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,kz=kz_YY_chiral)
case ('gaussian-x'); call gaussian(amplYY_chiral,f,iYY_chiral,kx=kx_YY_chiral)
case ('gaussian-y'); call gaussian(amplYY_chiral,f,iYY_chiral,ky=ky_YY_chiral)
case ('gaussian-z'); call gaussian(amplYY_chiral,f,iYY_chiral,kz=kz_YY_chiral)
case ('gaussian-noise'); call gaunoise(amplYY_chiral,f,iYY_chiral)
case ('positive-noise'); call posnoise(amplYY_chiral,f,iYY_chiral)
case ('wave-x'); call wave(amplYY_chiral,f,iYY_chiral,kx=kx_YY_chiral)
case ('wave-y'); call wave(amplYY_chiral,f,iYY_chiral,ky=ky_YY_chiral)
case ('wave-z'); call wave(amplYY_chiral,f,iYY_chiral,kz=kz_YY_chiral)
case ('cos(y-sinz)'); call cosyz_sinz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
case ('cosy_sinz'); call cosy_sinz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
case ('cosy_cosz'); call cosy_cosz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
case ('cosx_cosy_cosz'); call cosx_cosy_cosz(amplYY_chiral,f,iYY_chiral,kx_YY_chiral,ky_YY_chiral,kz_YY_chiral)
case ('cosx_siny_cosz'); call cosx_siny_cosz(amplYY_chiral,f,iYY_chiral,kx_YY_chiral,ky_YY_chiral,kz_YY_chiral)
case ('chiral_list'); call chiral_list(amplYY_chiral,f,iYY_chiral,'chiral_list')
case ('power_randomphase')
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iXX_chiral,iXX_chiral,lscale_tobox=.false.)
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iYY_chiral,iYY_chiral,lscale_tobox=.false.)
case ('power_randomphase_higgs')
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iXX_chiral,iXX_chiral,lscale_tobox=.false.)
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iYY_chiral,iYY_chiral,lscale_tobox=.false.)
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iXX2_chiral,iXX2_chiral,lscale_tobox=.false.)
call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,0.,cutoffYY_chiral,&
f,iYY2_chiral,iYY2_chiral,lscale_tobox=.false.)
!
! normalize
!
r=1./sqrt(f(:,:,:,iXX_chiral)**2 &
+f(:,:,:,iYY_chiral)**2 &
+f(:,:,:,iXX2_chiral)**2 &
+f(:,:,:,iYY2_chiral)**2)
f(:,:,:,iXX_chiral) =r*f(:,:,:,iXX_chiral)
f(:,:,:,iYY_chiral) =r*f(:,:,:,iYY_chiral)
f(:,:,:,iXX2_chiral)=r*f(:,:,:,iXX2_chiral)
f(:,:,:,iYY2_chiral)=r*f(:,:,:,iYY2_chiral)
!
! Higgs monopole
!
case ('higgs_mono')
pom=sqrt(spread(spread(x**2,2,my),3,mz) &
+spread(spread(y**2,1,mx),3,mz))
r=sqrt(spread(spread(x**2,2,my),3,mz) &
+spread(spread(y**2,1,mx),3,mz) &
+spread(spread(z**2,1,mx),2,my))
f(:,:,:,iXX_chiral) =sqrt(.5*(1.+spread(spread(z,1,mx),2,my)/r))
f(:,:,:,iYY_chiral) =0.
f(:,:,:,iXX2_chiral)=sqrt(.5*(1.-(spread(spread(z,1,mx),2,my)/r))) &
!-- *cos(atan2(spread(spread(y,1,mx),3,mz),spread(spread(x,2,my),3,mz)))
*spread(spread(x,2,my),3,mz)/pom
f(:,:,:,iYY2_chiral)=sqrt(.5*(1.-(spread(spread(z,1,mx),2,my)/r))) &
!-- *sin(atan2(spread(spread(y,1,mx),3,mz),spread(spread(x,2,my),3,mz)))
*spread(spread(y,1,mx),3,mz)/pom
case default; call fatal_error('init_chiral','no such initYY_chiral: '//trim(initYY_chiral))
endselect
!
! check next for initZZ_chiral
!
if (lZZ_chiral) then
select case (initZZ_chiral)
case ('zero'); f(:,:,:,iZZ_chiral)=0.
case ('const'); f(:,:,:,iZZ_chiral)=amplZZ_chiral
case ('chiral_list'); call chiral_list(amplZZ_chiral,f,iZZ_chiral,'chiral_listZZ')
case default; call fatal_error('init_chiral','no such initZZ_chiral: '//trim(initZZ_chiral))
endselect
endif
!
! Interface for user's own initial condition
!
if (linitial_condition) call initial_condition_chiral(f)
!
endsubroutine init_chiral
!***********************************************************************
subroutine chiral_list(fact,f,i,filename)
!
! Read intial conditions from file
!
! 13-apr-20/axel: coded
!
integer :: nrow, irow, ix, iy, i, l, m, n=1
character (len=*) :: filename
real, dimension (mx,my,mz,mfarray) :: f
real :: ampl, fact
!
open(1,file=trim(filename)//'.dat')
read(1,*) nrow
do irow=1,nrow
read(1,*) ix,iy,ampl
if (ix/nx==ipx.and.iy/ny==ipy) then
l=l1+mod(ix,nx)
m=m1+mod(iy,ny)
n=n1
print*,'AXEL: ',ix,iy,ipx,ipy,l,m
f(l,m,n,i)=f(l,m,n,i)+fact*ampl
endif
enddo
close(1)
!
endsubroutine chiral_list
!***********************************************************************
subroutine pencil_criteria_chiral()
!
! All pencils that the Chiral module depends on are specified here.
!
! 21-nov-04/anders: coded
!
lpenc_requested(i_uu)=.true.
if (llorentzforceEP) lpenc_requested(i_rho1)=.true.
!
endsubroutine pencil_criteria_chiral
!***********************************************************************
subroutine pencil_interdep_chiral(lpencil_in)
!
! Interdependency among pencils provided by the Chiral module
! is specified here.
!
! 21-11-04/anders: coded
!
logical, dimension(npencils) :: lpencil_in
!
call keep_compiler_quiet(lpencil_in)
!
endsubroutine pencil_interdep_chiral
!***********************************************************************
subroutine calc_pencils_chiral(f,p)
!
! Calculate Chiral pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 21-11-04/anders: coded
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
intent(in) :: f,p
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(p)
!
endsubroutine calc_pencils_chiral
!***********************************************************************
subroutine dXY_chiral_dt(f,df,p)
!
! passive scalar evolution
! calculate chirality equations in reduced form; see q-bio/0401036
!
! 28-may-04/axel: adapted from pscalar
! 1-jul-09/axel: included gradX, gradY, and allowed for different reactions
!
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
intent(in) :: p
intent(inout) :: f
intent(inout) :: df
!
real, dimension (nx,3) :: bbEP, jjEP, jxbEP
real, dimension (nx) :: XX_chiral, YY_chiral, ZZ_chiral
real, dimension (nx) :: ugXX_chiral, dXX_chiral
real, dimension (nx) :: ugYY_chiral, dYY_chiral
real, dimension (nx) :: ugZZ_chiral, dZZ_chiral
real, dimension (nx) :: RRXX_chiral,XX2_chiral
real, dimension (nx) :: RRYY_chiral,YY2_chiral
real, dimension (nx) :: RR21_chiral
real, dimension (nx) :: diffus_chiral
real :: pp,qq
integer :: j
!
! identify module and boundary conditions
!
if (headtt.or.ldebug) print*,'SOLVE dXY_dt'
if (headtt) call identify_bcs('XX_chiral',iXX_chiral)
if (headtt) call identify_bcs('YY_chiral',iYY_chiral)
!
! gradient of passive scalar
!
call grad(f,iXX_chiral,gXX_chiral)
call grad(f,iYY_chiral,gYY_chiral)
!
! Add diffusion of imposed spatially constant gradient of X or Y.
! This makes sense mainly for periodic boundary conditions.
!
if (limposed_gradient) then
do j=1,3
gXX_chiral(:,j)=gXX_chiral(:,j)+gradX0(j)
gYY_chiral(:,j)=gYY_chiral(:,j)+gradY0(j)
enddo
endif
!
! Advection term.
!
if (lhydro) then
call u_dot_grad(f,iXX_chiral,gXX_chiral,p%uu,ugXX_chiral,UPWIND=lupw_chiral)
call u_dot_grad(f,iYY_chiral,gYY_chiral,p%uu,ugYY_chiral,UPWIND=lupw_chiral)
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)-ugXX_chiral
df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)-ugYY_chiral
endif
!
! Diffusion term.
!
call del2(f,iXX_chiral,del2XX_chiral)
call del2(f,iYY_chiral,del2YY_chiral)
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+chiral_diffXX*del2XX_chiral
df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+chiral_diff *del2YY_chiral
!
! Advection & diffusion for the containment in the SIR model.
!
if (lZZ_chiral) then
if (lhydro) then
call u_dot_grad(f,iZZ_chiral,gZZ_chiral,p%uu,ugZZ_chiral,UPWIND=lupw_chiral)
df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)-ugZZ_chiral
endif
call del2(f,iZZ_chiral,del2ZZ_chiral)
df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)+chiral_diffZZ*del2ZZ_chiral
endif
!
! For Euler Potentials, possibility to add Lorentz force
!
if (llorentzforceEP) then
if (lhydro) then
call cross(gXX_chiral,gYY_chiral,bbEP)
do j=1,3
jjEP(:,j)=gXX_chiral(:,j)*del2YY_chiral &
-gYY_chiral(:,j)*del2XX_chiral
enddo
call cross(jjEP,bbEP,jxbEP)
if (ldensity) call multsv_mn(p%rho1,jxbEP,jxbEP)
df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+jxbEP
endif
endif
!
! selection of different reaction terms
!
select case (chiral_reaction)
!
! BAHN model
!
case ('BAHN_model')
!
! reaction terms
! X^2/Rtilde^2 - X*R
! Y^2/Rtilde^2 - Y*R
!
! for finite crossinhibition (=kI/kS) and finite fidelity (=f) we have
! R --> RX=X+Y*kI/kS, R --> RY=Y+X*kI/kS, and
! X2tilde=X^2/2RX, Y2tilde=Y^2/2RY.
!
XX_chiral=f(l1:l2,m,n,iXX_chiral)
YY_chiral=f(l1:l2,m,n,iYY_chiral)
RRXX_chiral=XX_chiral+YY_chiral*chiral_crossinhibition
RRYY_chiral=YY_chiral+XX_chiral*chiral_crossinhibition
!
! abbreviations for quadratic quantities
!
XX2_chiral=.5*XX_chiral**2/max(RRXX_chiral, tini)
YY2_chiral=.5*YY_chiral**2/max(RRYY_chiral, tini)
RR21_chiral=1./max(XX2_chiral+YY2_chiral, tini)
!
! fidelity factor
!
pp=.5*(1.+chiral_fidelity)
qq=.5*(1.-chiral_fidelity)
!
! final reaction equation
!
dXX_chiral=(pp*XX2_chiral+qq*YY2_chiral)*RR21_chiral-XX_chiral*RRXX_chiral
dYY_chiral=(pp*YY2_chiral+qq*XX2_chiral)*RR21_chiral-YY_chiral*RRYY_chiral
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
case('fisher')
if (headtt) print*,"chiral_reaction='fishers equation'"
if (headtt) print*,"growth rate=", chiral_fishernu,"carrying capacity=",chiral_fisherK
XX_chiral=f(l1:l2,m,n,iXX_chiral)
YY_chiral=f(l1:l2,m,n,iYY_chiral)
dXX_chiral=(1.-XX_chiral/chiral_fisherK)*XX_chiral*chiral_fishernu
dYY_chiral=(1.-YY_chiral/chiral_fisherK)*YY_chiral*chiral_fishernu
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
case('SIR')
if (headtt) print*,"chiral_reaction='SIR equation'"
if (headtt) print*,"growth rate=", chiral_fishernu,"carrying capacity=",chiral_fisherK
!
XX_chiral=f(l1:l2,m,n,iXX_chiral)
YY_chiral=f(l1:l2,m,n,iYY_chiral)
!
dXX_chiral=-chiral_fishernu*XX_chiral*YY_chiral &
+chiral_fisherR2_tdep*(1.-(XX_chiral+YY_chiral))
dYY_chiral=+chiral_fishernu*XX_chiral*YY_chiral-chiral_fisherK*YY_chiral &
+chiral_fisherR*(1.-(XX_chiral+YY_chiral))
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
! Containment in the SIR model.
!
if (lZZ_chiral) then
ZZ_chiral=f(l1:l2,m,n,iZZ_chiral)
dXX_chiral=-chiral_fishermu*XX_chiral*ZZ_chiral
dZZ_chiral=+chiral_fishermu*XX_chiral*ZZ_chiral-chiral_fisherH*ZZ_chiral
df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)+dZZ_chiral
endif
!
! Alternative when nothing is set.
!
case ('nothing')
if (lroot.and.ip<=5) print*,"chiral_reaction='nothing'"
!
case default
call fatal_error('chiral_reaction','no such chiral_reaction: '//trim(chiral_reaction))
endselect
!
! For the timestep calculation, need maximum diffusion
!
if (lupdate_courant_dt) then
diffus_chiral=max(chiral_diffXX,chiral_diff)*dxyz_2
if (headtt.or.ldebug) print*,'dXY_chiral_dt: max(diffus_chiral) =',maxval(diffus_chiral)
maxdiffus=max(maxdiffus,diffus_chiral)
endif
call calc_diagnostics_chiral(f,p)
endsubroutine dXY_chiral_dt
!***********************************************************************
subroutine calc_diagnostics_chiral(f,p)
!
! diagnostics
!
! output for double and triple correlators (assume z-gradient of cc)
! <u_k u_j d_j c> = <u_k c uu.gradXX_chiral>
!
use Diagnostics
use Sub
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
real, dimension (nx) :: XX_chiral,YY_chiral,ZZ_chiral,QQ_chiral,QQ21QQ_chiral,bbEP2,jjEP2,jbEP
real, dimension (nx,3) :: bbEP, jjEP
real, dimension (nx,3,3) :: gXXij_chiral,gYYij_chiral
real :: lamchiral
integer :: j
if (ldiagnos) then
XX_chiral=f(l1:l2,m,n,iXX_chiral)
YY_chiral=f(l1:l2,m,n,iYY_chiral)
call max_mn_name(XX_chiral,idiag_XX_chiralmax)
call max_mn_name(YY_chiral,idiag_YY_chiralmax)
call sum_mn_name(XX_chiral,idiag_XX_chiralm)
call sum_mn_name(YY_chiral,idiag_YY_chiralm)
if (lZZ_chiral) then
ZZ_chiral=f(l1:l2,m,n,iZZ_chiral)
call sum_mn_name(ZZ_chiral,idiag_ZZ_chiralm)
call max_mn_name(ZZ_chiral,idiag_ZZ_chiralmax)
endif
call save_name(chiral_fisherR2_tdep,idiag_R2tdep)
!
if (idiag_QQm_chiral/=0 .or. idiag_QQ21m_chiral/=0 .or. idiag_QQ21QQm_chiral/=0) then
QQ_chiral=XX_chiral-YY_chiral
call sum_mn_name(QQ_chiral,idiag_QQm_chiral)
if (idiag_QQ21m_chiral/=0) call sum_mn_name(1.-QQ_chiral**2,idiag_QQ21m_chiral)
if (idiag_QQ21QQm_chiral/=0) then
lamchiral=2.*chiral_fidelity-1.
QQ21QQ_chiral=(lamchiral-QQ_chiral**2)/(1.+QQ_chiral**2)*QQ_chiral
call sum_mn_name(QQ21QQ_chiral,idiag_QQ21QQm_chiral)
endif
endif
!
! Calculate BB = gXX_chiral x gYY_chiral
!
if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0.or.idiag_jbmEP/=0) then
call cross(gXX_chiral,gYY_chiral,bbEP)
if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0) then
call dot2(bbEP,bbEP2)
call sum_mn_name(bbEP2,idiag_brmsEP,lsqrt=.true.)
call max_mn_name(bbEP2,idiag_bmaxEP,lsqrt=.true.)
endif
endif
!
! Calculate Ji = Xi*del2Y - Yi*del2X + Xij Yj - Yij Xj
! and then the max and rms values of that
!
if (idiag_jrmsEP/=0.or.idiag_jmaxEP/=0.or.idiag_jbmEP/=0) then
do j=1,3
jjEP(:,j)=gXX_chiral(:,j)*del2YY_chiral &
-gYY_chiral(:,j)*del2XX_chiral
enddo
call g2ij(f,iXX_chiral,gXXij_chiral)
call g2ij(f,iYY_chiral,gYYij_chiral)
call multmv_mn(+gXXij_chiral,gYY_chiral,jjEP,ladd=.true.)
call multmv_mn(-gYYij_chiral,gXX_chiral,jjEP,ladd=.true.)
call dot2(jjEP,jjEP2)
call sum_mn_name(jjEP2,idiag_jrmsEP,lsqrt=.true.)
call max_mn_name(jjEP2,idiag_jmaxEP,lsqrt=.true.)
!
! For J.B, we need B, but only if not already calculated.
!
if (idiag_jbmEP/=0) then
call dot(jjEP,bbEP,jbEP)
call sum_mn_name(jbEP,idiag_jbmEP)
endif
endif
endif
!
endsubroutine calc_diagnostics_chiral
!***********************************************************************
subroutine chiral_before_boundary(f)
!
! initialize aa from Euler potentials (EP), but only once when t=0
!
! 4-jul-09/axel: coded
!
use Cdata
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,3) :: gXX_chiral,gYY_chiral,gXX2_chiral,gYY2_chiral
integer :: j
!
intent(inout) :: f
!
! Possibility to initialize magnetic vector potential in terms
! of Euler potentials X and Y.
!
if (linitialize_aa_from_EP) then
if (t>=tinitialize_aa_from_EP) then
linitialize_aa_from_EP=.false.
if (lmagnetic) then
if (lroot) print*,'initialize_aa_from_EP at t=',t
do n=n1,n2
do m=m1,m2
call grad(f,iXX_chiral,gXX_chiral)
call grad(f,iYY_chiral,gYY_chiral)
if (initYY_chiral=='power_randomphase_higgs') then
call grad(f,iXX2_chiral,gXX2_chiral)
call grad(f,iYY2_chiral,gYY2_chiral)
endif
if (limposed_gradient) then
do j=1,3
gXX_chiral(:,j)=gXX_chiral(:,j)+gradX0(j)
gYY_chiral(:,j)=gYY_chiral(:,j)+gradY0(j)
enddo
endif
do j=1,3
if (linitialize_aa_from_EP_alpgrad) then
f(l1:l2,m,n,j+iaa-1)=f(l1:l2,m,n,iXX_chiral)*gYY_chiral(:,j)
elseif (linitialize_aa_from_EP_betgrad) then
f(l1:l2,m,n,j+iaa-1)=-f(l1:l2,m,n,iYY_chiral)*gXX_chiral(:,j)
else
f(l1:l2,m,n,j+iaa-1)=.5*( f(l1:l2,m,n,iXX_chiral)*gYY_chiral(:,j) &
-f(l1:l2,m,n,iYY_chiral)*gXX_chiral(:,j))
if (initYY_chiral=='power_randomphase_higgs') then
f(l1:l2,m,n,j+iaa-1)=f(l1:l2,m,n,j+iaa-1) &
+.5*( f(l1:l2,m,n,iXX2_chiral)*gYY2_chiral(:,j) &
-f(l1:l2,m,n,iYY2_chiral)*gXX2_chiral(:,j))
if (lroot) print*,'initialize_aa_from_Higgs at t=',t
endif
endif
enddo
enddo
enddo
endif
endif
endif
!
! time-dependent reinfection rate
!
if (chiral_fisherR2_tend==0.) then
chiral_fisherR2_tdep=chiral_fisherR2
else
chiral_fisherR2_tdep = chiral_fisherR2*(max(0.,1.-(max(0.,real(t)-chiral_fisherR2_tstart)/ &
(chiral_fisherR2_tend-chiral_fisherR2_tstart))**2)**2)
endif
endsubroutine chiral_before_boundary
!***********************************************************************
subroutine read_chiral_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=chiral_init_pars, IOSTAT=iostat)
!
endsubroutine read_chiral_init_pars
!***********************************************************************
subroutine write_chiral_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=chiral_init_pars)
!
endsubroutine write_chiral_init_pars
!***********************************************************************
subroutine read_chiral_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=chiral_run_pars, IOSTAT=iostat)
!
endsubroutine read_chiral_run_pars
!***********************************************************************
subroutine write_chiral_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=chiral_run_pars)
!
endsubroutine write_chiral_run_pars
!***********************************************************************
subroutine rprint_chiral(lreset,lwrite)
!
! reads and registers print parameters relevant for chirality
!
! 28-may-04/axel: adapted from pscalar
!
use Diagnostics
!
integer :: iname
logical :: lreset
logical, optional :: lwrite
!
! reset everything in case of reset
! (this needs to be consistent with what is defined above!)
!
if (lreset) then
idiag_XX_chiralmax=0; idiag_XX_chiralm=0
idiag_YY_chiralmax=0; idiag_YY_chiralm=0
idiag_ZZ_chiralmax=0; idiag_ZZ_chiralm=0
idiag_QQm_chiral=0; idiag_QQ21m_chiral=0; idiag_QQ21QQm_chiral=0
idiag_brmsEP=0; idiag_bmaxEP=0
idiag_jrmsEP=0; idiag_jmaxEP=0
idiag_jbmEP=0; idiag_R2tdep=0
endif
!
! check for those quantities that we want to evaluate online
!
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),&
'XXm',idiag_XX_chiralm)
call parse_name(iname,cname(iname),cform(iname),&
'YYm',idiag_YY_chiralm)
call parse_name(iname,cname(iname),cform(iname),&
'ZZm',idiag_ZZ_chiralm)
call parse_name(iname,cname(iname),cform(iname),&
'XXmax',idiag_XX_chiralmax)
call parse_name(iname,cname(iname),cform(iname),&
'YYmax',idiag_YY_chiralmax)
call parse_name(iname,cname(iname),cform(iname),&
'ZZmax',idiag_ZZ_chiralmax)
call parse_name(iname,cname(iname),cform(iname),&
'QQm',idiag_QQm_chiral)
call parse_name(iname,cname(iname),cform(iname),&
'QQ21m',idiag_QQ21m_chiral)
call parse_name(iname,cname(iname),cform(iname),&
'QQ21QQm',idiag_QQ21QQm_chiral)
call parse_name(iname,cname(iname),cform(iname),&
'brmsEP',idiag_brmsEP)
call parse_name(iname,cname(iname),cform(iname),&
'bmaxEP',idiag_bmaxEP)
call parse_name(iname,cname(iname),cform(iname),&
'jrmsEP',idiag_jrmsEP)
call parse_name(iname,cname(iname),cform(iname),&
'jmaxEP',idiag_jmaxEP)
call parse_name(iname,cname(iname),cform(iname),&
'jbmEP',idiag_jbmEP)
call parse_name(iname,cname(iname),cform(iname),&
'R2tdep',idiag_R2tdep)
enddo
!
! check for those quantities for which we want video slices
!
if (lwrite_slices) then
where(cnamev=='XX_chiral' .or. &
cnamev=='YY_chiral' .or. &
cnamev=='ZZ_chiral' .or. &
cnamev=='DQ_chiral') &
cformv='DEFINED'
endif
!
endsubroutine rprint_chiral
!***********************************************************************
subroutine get_slices_chiral(f,slices)
!
! Write slices for animation of Chiral variables.
!
! 26-jul-06/tony: coded
!
use Slices_methods, only: assign_slices_scal
real, dimension (mx,my,mz,mfarray) :: f
type (slice_data) :: slices
!
! Loop over slices
!
select case (trim(slices%name))
!
! Chirality fields: XX
!
case ('XX_chiral'); call assign_slices_scal(slices,f,iXX_chiral)
!
! Chirality fields: YY
!
case ('YY_chiral'); call assign_slices_scal(slices,f,iYY_chiral)
!
! Chirality fields: ZZ
!
case ('ZZ_chiral'); call assign_slices_scal(slices,f,iZZ_chiral)
!
! Chirality fields: DQ
!
case ('DQ_chiral')
if (lwrite_slice_yz ) &
call QQ_chiral(f(ix_loc,m1:m2,n1:n2,iXX_chiral)-f(ix_loc,m1:m2,n1:n2,iYY_chiral),slices%yz)
if (lwrite_slice_xz ) &
call QQ_chiral(f(l1:l2,iy_loc,n1:n2,iXX_chiral)-f(l1:l2,iy_loc,n1:n2,iYY_chiral),slices%xz)
if (lwrite_slice_xy ) &
call QQ_chiral(f(l1:l2,m1:m2,iz_loc,iXX_chiral)-f(l1:l2,m1:m2,iz_loc,iYY_chiral),slices%xy)
if (lwrite_slice_xy2) &
call QQ_chiral(f(l1:l2,m1:m2,iz2_loc,iXX_chiral)-f(l1:l2,m1:m2,iz2_loc,iYY_chiral),slices%xy2)
if (lwrite_slice_xy3) &
call QQ_chiral(f(l1:l2,m1:m2,iz3_loc,iXX_chiral)-f(l1:l2,m1:m2,iz3_loc,iYY_chiral),slices%xy3)
if (lwrite_slice_xy4) &
call QQ_chiral(f(l1:l2,m1:m2,iz4_loc,iXX_chiral)-f(l1:l2,m1:m2,iz4_loc,iYY_chiral),slices%xy4)
if (lwrite_slice_xz2) &
call QQ_chiral(f(l1:l2,iy2_loc,n1:n2,iXX_chiral)-f(l1:l2,iy2_loc,n1:n2,iYY_chiral),slices%xz2)
slices%ready=.true.
!
endselect
!
contains
subroutine QQ_chiral(source,dest)
!
real, dimension(:,:) :: source,dest
!
dest=source**2
dest=source*(1.-dest)/(1.+dest)
!
endsubroutine QQ_chiral
endsubroutine get_slices_chiral
!***********************************************************************
endmodule Chiral