! $Id$
!
! This module provide a way for users to specify custom
! (i.e. not in the standard Pencil Code) physics, diagnostics etc.
!
! The module provides a set of standard hooks into the Pencil-Code and
! currently allows the following customizations:
!
! Description | Relevant function call
! ---------------------------------------------------------------------------
! Special variable registration | register_special
! (pre parameter read) |
! Special variable initialization | initialize_special
! (post parameter read) |
! Special variable finalization | finalize_special
! (deallocation, etc.) |
! |
! Special initial condition | init_special
! this is called last so may be used to modify |
! the mvar variables declared by this module |
! or optionally modify any of the other f array |
! variables. The latter, however, should be |
! avoided where ever possible. |
! |
! Special term in the mass (density) equation | special_calc_density
! Special term in the momentum (hydro) equation | special_calc_hydro
! Special term in the energy equation | special_calc_energy
! Special term in the induction (magnetic) | special_calc_magnetic
! equation |
! |
! Special equation | dspecial_dt
! NOT IMPLEMENTED FULLY YET - HOOKS NOT PLACED INTO THE PENCIL-CODE
!
!** 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 :: lspecial = .true.
!
! MVAR CONTRIBUTION 12
! MAUX CONTRIBUTION 3
!
!***************************************************************
!
! HOW TO USE THIS FILE
! --------------------
!
! Change the line above to
! lspecial = .true.
! to enable use of special hooks.
!
! The rest of this file may be used as a template for your own
! special module. Lines which are double commented are intended
! as examples of code. Simply fill out the prototypes for the
! features you want to use.
!
! Save the file with a meaningful name, eg. geo_kws.f90 and place
! it in the $PENCIL_HOME/src/special directory. This path has
! been created to allow users ot optionally check their contributions
! in to the Pencil-Code SVN repository. This may be useful if you
! are working on/using the additional physics with somebodyelse or
! may require some assistance from one of the main Pencil-Code team.
!
! To use your additional physics code edit the Makefile.local in
! the src directory under the run directory in which you wish to
! use your additional physics. Add a line with all the module
! selections to say something like:
!
! SPECIAL=special/geo_kws
!
! Where geo_kws it replaced by the filename of your new module
! upto and not including the .f90
!
module Special
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages, only: svn_id, fatal_error
!
implicit none
!
include '../special.h'
!
! Declare index of new variables in f array (if any).
!
integer :: ihij,igij
!
!! Diagnostic variables (needs to be consistent with reset list below).
!
integer :: idiag_g22pt=0 ! DIAG_DOC: $g_{22}(x_1,y_1,z_1,t)$
!
contains
!***********************************************************************
subroutine register_special()
!
! Set up indices for variables in special modules.
!
! 6-oct-03/tony: coded
!
use Sub, only: register_report_aux
use FArrayManager
!
if (lroot) call svn_id( &
"$Id$")
!
call farray_register_pde('hij',ihij,vector=6)
call farray_register_pde('gij',igij,vector=6)
!
! Set indices for auxiliary variables.
!
!call farray_register_auxiliary('bb',ibb)
call register_report_aux('bb', ibb, ibx, iby, ibz)
print*,'AXEL1: registered, ibb, ibx, iby, ibz=',ibb, ibx, iby, ibz
!
if (lroot) call svn_id( &
"$Id$")
!
endsubroutine register_special
!***********************************************************************
subroutine initialize_special(f)
!
! Called after reading parameters, but before the time loop.
!
! 06-oct-03/tony: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
call keep_compiler_quiet(f)
!
if (lfargo_advection) then
print*,''
print*,'Switch '
print*,' SPECIAL = special/fargo'
print*,'in src/Makefile.local if you want to use the fargo algorithm'
print*,''
call fatal_error('nospecial','initialize_special()')
endif
!
endsubroutine initialize_special
!***********************************************************************
subroutine finalize_special(f)
!
! Called right before exiting.
!
! 14-aug-2011/Bourdin.KIS: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine finalize_special
!***********************************************************************
subroutine init_special(f)
!
! initialise special condition; called from start.f90
! 06-oct-2003/tony: coded
!
real, dimension (mx,my,mz,mfarray) :: f
!
intent(inout) :: f
!!
!! SAMPLE IMPLEMENTATION
!!
!! select case (initspecial)
!! case ('nothing'); if (lroot) print*,'init_special: nothing'
!! case ('zero', '0'); f(:,:,:,iSPECIAL_VARIABLE_INDEX) = 0.
!! case default
!! call fatal_error("init_special: No such value for initspecial:" &
!! ,trim(initspecial))
!! endselect
!
call keep_compiler_quiet(f)
!
endsubroutine init_special
!***********************************************************************
subroutine pencil_criteria_special()
!
! All pencils that this special module depends on are specified here.
!
! 18-07-06/tony: coded
!
lpenc_requested(i_bb)=.true.
lpenc_requested(i_b2)=.true.
!
endsubroutine pencil_criteria_special
!***********************************************************************
subroutine pencil_interdep_special(lpencil_in)
!
! Interdependency among pencils provided by this module are specified here.
!
! 18-07-06/tony: coded
!
logical, dimension(npencils), intent(inout) :: lpencil_in
!
call keep_compiler_quiet(lpencil_in)
!
endsubroutine pencil_interdep_special
!***********************************************************************
subroutine calc_pencils_special(f,p)
!
! Calculate Special pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 24-nov-04/tony: coded
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
intent(in) :: f
intent(inout) :: p
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(p)
!
endsubroutine calc_pencils_special
!***********************************************************************
subroutine dspecial_dt(f,df,p)
!
! calculate right hand side of ONE OR MORE extra coupled PDEs
! along the 'current' Pencil, i.e. f(l1:l2,m,n) where
! m,n are global variables looped over in equ.f90
!
! Due to the multi-step Runge Kutta timestepping used one MUST always
! add to the present contents of the df array. NEVER reset it to zero.
!
! Several precalculated Pencils of information are passed for
! efficiency.
!
! 06-oct-03/tony: coded
!
use Diagnostics
use Sub, only: del2v
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (nx,3) :: del2hii,del2hij
type (pencil_case) :: p
!
integer :: j,jhij,jgij
!
intent(in) :: f,p
intent(inout) :: df
!
! Identify module and boundary conditions.
!
if (headtt.or.ldebug) print*,'dspecial_dt: SOLVE dspecial_dt'
!! if (headtt) call identify_bcs('special',ispecial)
!
if (lmagnetic) then
! print*,'AXEL:',p%bb(:,2)
!
! g11=1, g22=2, g33=3, g12=4, g13=5, g23=6,
! g11=0, g22=1, g33=2, g12=3, g13=4, g23=5,
!
do j=1,6
jhij=ihij-1+j
jgij=igij-1+j
df(l1:l2,m,n,jhij)=df(l1:l2,m,n,jhij)+f(l1:l2,m,n,jgij)
enddo
call del2v(f,ihij ,del2hii)
call del2v(f,ihij+3,del2hij)
df(l1:l2,m,n,igij+0)=df(l1:l2,m,n,igij+0)+del2hii(:,1)+ &
p%bb(:,1)**2-onethird*p%b2
df(l1:l2,m,n,igij+1)=df(l1:l2,m,n,igij+1)+del2hii(:,2)+ &
p%bb(:,2)**2-onethird*p%b2
df(l1:l2,m,n,igij+2)=df(l1:l2,m,n,igij+2)+del2hii(:,3)+ &
p%bb(:,3)**2-onethird*p%b2
df(l1:l2,m,n,igij+3)=df(l1:l2,m,n,igij+3)+del2hij(:,1)+p%bb(:,1)*p%bb(:,2)
df(l1:l2,m,n,igij+4)=df(l1:l2,m,n,igij+4)+del2hij(:,2)+p%bb(:,1)*p%bb(:,3)
df(l1:l2,m,n,igij+5)=df(l1:l2,m,n,igij+5)+del2hij(:,3)+p%bb(:,2)*p%bb(:,3)
else
call fatal_error("dspecial_dt","need magnetic field")
endif
!
! diagnostics
!
if (ldiagnos) then
if (lroot.and.m==mpoint.and.n==npoint) then
if (idiag_g22pt/=0) call save_name(p%bb(lpoint-nghost,2),idiag_g22pt)
endif
endif
!
call keep_compiler_quiet(f,df)
call keep_compiler_quiet(p)
!
endsubroutine dspecial_dt
!***********************************************************************
subroutine calc_lspecial_pars(f)
!
! dummy routine
!
! 15-jan-08/axel: coded
!
use Fourier, only: fourier_transform
!
real, dimension (:,:,:,:), allocatable :: B_re, B_im, v_re, v_im
real, dimension (:,:,:), allocatable :: k2, r
real, dimension (:), allocatable :: kx, ky, kz
real, dimension (mx,my,mz,mfarray) :: f
integer :: i,ikx,iky,ikz,stat
logical :: lscale_tobox1=.true.
real :: scale_factor
intent(inout) :: f
!
! Allocate memory for arrays.
!
allocate(k2(nx,ny,nz),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for k2')
allocate(r(nx,ny,nz),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for r')
!
allocate(B_re(nx,ny,nz,3),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for B_re')
allocate(B_im(nx,ny,nz,3),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for B_im')
!
allocate(v_re(nx,ny,nz,3),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for v_re')
allocate(v_im(nx,ny,nz,3),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars','Could not allocate memory for v_im')
!
allocate(kx(nxgrid),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars', &
'Could not allocate memory for kx')
allocate(ky(nygrid),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars', &
'Could not allocate memory for ky')
allocate(kz(nzgrid),stat=stat)
if (stat>0) call fatal_error('calc_lspecial_pars', &
'Could not allocate memory for kz')
!
! calculate k^2
!
scale_factor=1
if (lscale_tobox1) scale_factor=2*pi/Lx
kx=cshift((/(i-(nxgrid+1)/2,i=0,nxgrid-1)/),+(nxgrid+1)/2)*scale_factor
!
scale_factor=1
if (lscale_tobox1) scale_factor=2*pi/Ly
ky=cshift((/(i-(nygrid+1)/2,i=0,nygrid-1)/),+(nygrid+1)/2)*scale_factor
!
scale_factor=1
if (lscale_tobox1) scale_factor=2*pi/Lz
kz=cshift((/(i-(nzgrid+1)/2,i=0,nzgrid-1)/),+(nzgrid+1)/2)*scale_factor
!
! Set k^2 array. Note that in Fourier space, kz is the fastest index and has
! the full nx extent (which, currently, must be equal to nxgrid).
!
if (lroot .AND. ip<10) &
print*,'calc_lspecial_pars:fft ...'
do iky=1,nz
do ikx=1,ny
do ikz=1,nx
k2(ikz,ikx,iky)=kx(ikx+ipy*ny)**2+ky(iky+ipz*nz)**2+kz(ikz+ipx*nx)**2
enddo
enddo
enddo
if (lroot) k2(1,1,1) = 1. ! Avoid division by zero
!
! Compute Mij(x,t)=Bi(x,t)*Bj(x,t) in real space
!
! Find bb if as communicated auxiliary.
!
call zero_ghosts(f, iax, iaz)
call update_ghosts(f, iax, iaz)
mn_loop: do imn = 1, ny * nz
m = mm(imn)
n = nn(imn)
call gij(f, iaa, aij, 1)
call curl_mn(aij, bb, f(l1:l2,m,n,iax:iaz))
!
! Add imposed field, if any
!
bext: if (lB_ext_in_comaux) then
call get_bext(b_ext)
forall(j = 1:3, b_ext(j) /= 0.0) bb(:,j) = bb(:,j) + b_ext(j)
if (headtt .and. imn == 1) print *, 'magnetic_before_boundary: B_ext
= ', b_ext
endif bext
f(l1:l2,m,n,ibx:ibz) = bb
enddo mn_loop
endif getbb
!
! Go into Fourier space
!
v_im=0.0
v_re=f(l1:l2,m1:m2,n1:n2,iax:iaz)
do i=1,3
call fourier_transform(v_re(:,:,:,i),v_im(:,:,:,i))
enddo !i
!
! Compute the magnetic field in Fourier space from vector potential.
! In this definition of the Fourier transform, nabla = +ik.
!
do ikz=1,nz; do iky=1,ny; do ikx=1,nx
!
! (vx, vy, vz) -> Bx
!
B_re(ikz,ikx,iky,1)=+( &
-ky(iky+ipz*nz)*v_im(ikz,ikx,iky,3) &
+kz(ikz+ipx*nx)*v_im(ikz,ikx,iky,2) )
B_im(ikz,ikx,iky,1)=+( &
+ky(iky+ipz*nz)*v_re(ikz,ikx,iky,3) &
-kz(ikz+ipx*nx)*v_re(ikz,ikx,iky,2) )
!
! (vx, vy, vz) -> By
!
B_re(ikz,ikx,iky,2)=+( &
-kz(ikz+ipx*nx)*v_im(ikz,ikx,iky,1) &
+kx(ikx+ipy*ny)*v_im(ikz,ikx,iky,3) )
B_im(ikz,ikx,iky,2)=+( &
+kz(ikz+ipx*nx)*v_re(ikz,ikx,iky,1) &
-kx(ikx+ipy*ny)*v_re(ikz,ikx,iky,3) )
!
! (vx, vy, vz) -> Bz
!
B_re(ikz,ikx,iky,3)=+( &
-kx(ikx+ipy*ny)*v_im(ikz,ikx,iky,2) &
+ky(iky+ipz*nz)*v_im(ikz,ikx,iky,1) )
B_im(ikz,ikx,iky,3)=+( &
+kx(ikx+ipy*ny)*v_re(ikz,ikx,iky,2) &
-ky(iky+ipz*nz)*v_re(ikz,ikx,iky,1) )
!
enddo; enddo; enddo
!
! back to real space
!
print*,'AXEL2: registered, ibb, ibx, iby, ibz=',ibb, ibx, iby, ibz
do i=1,3
call fourier_transform(B_re(:,:,:,i),B_im(:,:,:,i),linv=.true.)
f(l1:l2,m1:m2,n1:n2,ibb+i-1)=B_re(:,:,:,i)
enddo !i
!
endsubroutine calc_lspecial_pars
!***********************************************************************
subroutine rprint_special(lreset,lwrite)
!
! Reads and registers print parameters relevant to special.
!
! 06-oct-03/tony: coded
!
use Diagnostics
!
integer :: iname
logical :: lreset,lwr
logical, optional :: lwrite
!
lwr = .false.
if (present(lwrite)) lwr=lwrite
!!!
!!! reset everything in case of reset
!!! (this needs to be consistent with what is defined above!)
!!!
if (lreset) then
idiag_g22pt=0
endif
!
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'g22pt',idiag_g22pt)
enddo
!!
!!! write column where which magnetic variable is stored
!! if (lwr) then
!! write(3,*) 'i_SPECIAL_DIAGNOSTIC=',i_SPECIAL_DIAGNOSTIC
!! endif
!!
endsubroutine rprint_special
!***********************************************************************
!************ DO NOT DELETE THE FOLLOWING **************
!********************************************************************
!** This is an automatically generated include file that creates **
!** copies dummy routines from nospecial.f90 for any Special **
!** routines not implemented in this file **
!** **
include '../special_dummies.inc'
!********************************************************************
endmodule Special