! $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