The setup is represented by two files in your simulation-directory, amrvacusr.t and amrvacusrpar.t that can be adopted from specific examples in tests/rmhd/.../amrvacusr.t and tests/rmhd/.../amrvacusrpar.t.
!=============================================================================
! amrvacusr.t
!=============================================================================
INCLUDE:amrvacnul/speciallog.t
INCLUDE:amrvacnul/specialbound.t
INCLUDE:amrvacnul/specialsource.t
INCLUDE:amrvacnul/specialimpl.t
INCLUDE:amrvacnul/usrflags.t
INCLUDE:amrvacnul/correctaux_usr.t
!=============================================================================
subroutine initglobaldata_usr
include 'amrvacdef.f'
!-----------------------------------------------------------------------------
! Here you set global parameters, e.g. the adiabatic index:
! eqpar(gamma_) = 5.0d0/3.0d0
end subroutine initglobaldata_usr
!=============================================================================
subroutine initonegrid_usr(ixI^L,ixO^L,s)
! initialize one grid within ixO^L
include 'amrvacdef.f'
integer, intent(in) :: ixI^L, ixO^L
type(state) :: s
!-----------------------------------------------------------------------------
associate(x=>s%x%x,w=>s%w%w{#IFDEF STAGGERED ,ws=>s%ws%w})
! Set your initial conditions here
! e.g.:
! w(ixG^S,rho_) = 1.0d0
! w(ixO^S,pp_) = 1.0d0
! w(ixO^S,u1_) = 0.0d0
! w(ixO^S,u2_) = 0.0d0
! w(ixO^S,u3_) = 0.0d0
! For divB-free fields, obtain B-fields like this:
! {#IFNDEF STAGGERED
! call b_from_vectorpotential(ixI^L,ixO^L,w,x)
! }{#IFDEF STAGGERED
! call b_from_vectorpotential(s%ws%ixG^L,ixI^L,ixO^L,ws,x)
! call faces2centers(ixO^L,s)
! }
! Eventually convert to primitive variables:
! call conserve(ixI^L,ixO^L,w,x,patchfalse)
end associate
end subroutine initonegrid_usr
!=============================================================================
subroutine initvecpot_usr(ixI^L, ixC^L, xC, A, idir)
! initialize the vectorpotential on the corners
! used by b_from_vectorpotential()
include 'amrvacdef.f'
integer, intent(in) :: ixI^L, ixC^L, idir
double precision, intent(in) :: xC(ixI^S,1:ndim)
double precision, intent(out) :: A(ixI^S)
!-----------------------------------------------------------------------------
! Set your vectorpotential in direction idir here
A(ixC^S) = zero
end subroutine initvecpot_usr
!=============================================================================
! amrvacusr.t
!=============================================================================
Now you should edit these subroutines according to your needs: the idea is that in
initglobaldata_usr you must set the global equation parameter values (i.e. all eqpar(*) entries. In the subroutine initonegrid_usr, you have to make sure that at the end of this subroutine, all conservative variable values are provided on the full grid, i.e. you need to specify physically meaningfull w entries. You have the grid info available in the x variable.
You can write the subroutine(s) either in the dimension independent notation, described in source, or in Fortran 90 if the number of dimensions is fixed for your problem.
When the predefined boundary types provided by BHAC are not sufficient, the specialbound_usr subroutine can solve the problem. It is called for each boundary region for which at least one boundary type typeB='special'... is selected. The specialbound_usr subroutine is called after the variables with predefined boundary types have been updated in the ghost cells.
!=============================================================================
subroutine specialbound_usr(qt,ixI^L,ixO^L,iB,s)
! special boundary types, user defined
! user must assign conservative variables in bounderies
include 'amrvacdef.f'
integer, intent(in) :: ixI^L, ixO^L, iB
double precision, intent(in) :: qt
type(state), intent(inout) :: s
!----------------------------------------------------------------------------
associate(x=>s%x%x,w=>s%w%w{#IFDEF STAGGERED ,ws=>s%ws%w})
call mpistop("specialbound not defined")
! implement boundary conditions for a side, for example
! set the radial momentum at the inner radial boundary (assuming spherical-
! like coordinates) to zero:
! select case(iB)
! case(1)
! w(ixO^S,s1_) = 0.0d0
! end select
end associate
end subroutine specialbound_usr
!=============================================================================
You need to provide conserved variables (or call convert at the end of the subroutine, but then make sure all variables are actually in primitive form prior to the call).
In any case the subroutine should integrate w from time qt to qt+qdt for the variables listed in iw^LIM in the region ixO^S. The source terms should be evaluated for the wCT array, which corresponds to the physical time qtC. In case of explicit time dependence, qtC should be used as time. Only elements within ixI^S can be used from wCT.
The default subroutine looks like this:
!=============================================================================
subroutine specialsource(qdt,ixI^L,ixO^L,iw^LIM,qtC,wCT,qt,w,x)
! Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices
! iw=iwmin...iwmax. wCT is at time qCT
include 'amrvacdef.f'
integer, intent(in) :: ixI^L, ixO^L, iw^LIM
double precision, intent(in) :: qdt, qtC, qt, x(ixI^S,1:ndim)
double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
! integer :: iw
! double precision :: s(ixG^T)
!-----------------------------------------------------------------------------
! do iw= iw^LIM
! select case(iw)
! case(m1_)
! ! The source is based on the time centered wCT
! call getmyforce(wCT,ixO^L,s)
! w(ixO^S,m1_)=w(ixO^S,m1_) + qdt*s(ixO^S)
! case(e_)
! call getmyheating(wCT,ixO^L,s)
! w(ixO^S,e_) =w(ixO^S,e_) + qdt*s(ixO^S)
! end select
! end do
end subroutine specialsource
!=============================================================================
where in the commented out example 'force' and 'heating' depend on wCT as needed. Together with the special source terms, a timestep limitation can be implemented in getdt_special:
!=============================================================================
subroutine getdt_special(w,ixI^L,ixO^L,dtnew,dx^D,x)
! Limit "dt" further if necessary, e.g. due to the special source terms.
! The getdt_courant (CFL condition) and the getdt subroutine in the AMRVACPHYS
! module have already been called.
include 'amrvacdef.f'
integer, intent(in) :: ixI^L, ixO^L
double precision, intent(in) :: dx^D, x(ixI^S,1:ndim)
double precision, intent(inout) :: w(ixI^S,1:nw), dtnew
!-----------------------------------------------------------------------------
dtnew=bigdouble
end subroutine getdt_special
!=============================================================================
Any dtnew smaller than 'bigdouble' will be taken into account in the calculation of the timestep. This subroutine is called
after the CFL condition getdt_courant and the getdt
subroutine of the AMRVACPHYS module have been executed.
Two more subroutines from amrvacnul/specialsource.t should be pointed out:
The specialrefine_grid subroutine allows to add user controlled (de)refinement, by setting the integers refine,coarsen. You have all info available to do this refining (grid level, physical values in w, coordinates in x, time in qt). Similarly, the specialvarforerrest subroutine allows to compute a (local) new variable, to be stored in var, on which you can then base refinement as well. This is true for the lohner error estimator only.
The most important one is the specialvar_output subroutine which is extremely handy to compute variables from the actually computed conserved variables, that can then be visualized directly. It is only used in combination with the conversion subroutines. E.g., one may here compute current density components using the actual code discretizations for computing a curl, and then visualize those with any of the visualization tools applicable. You then also need to specify a label for this variable, in specialvarnames_output.
For example, compute the divergence of the magnetic field as extra (nwauxio) variable as follows:
!=============================================================================
subroutine specialvar_output(ixI^L,ixO^L,nwmax,w,s,normconv)
! this subroutine can be used in convert, to add auxiliary variables to the
! converted output file, for further analysis using tecplot, paraview, ....
! these auxiliary values need to be stored in the nw+1:nw+nwauxio slots
!
! the array normconv can be filled in the (nw+1:nw+nwauxio) range with
! corresponding normalization values (default value 1)
include 'amrvacdef.f'
integer, intent(in) :: ixI^L,ixO^L,nwmax
double precision :: w(ixI^S,nw+nwmax)
type(state) :: s
double precision :: normconv(0:nw+nwmax)
! .. local ..
!-----------------------------------------------------------------------------
associate(x=>s%x%x{#IFDEF STAGGERED ,ws=>s%ws%w})
{#IFDEF STAGGERED
call div_staggered(ixO^L,s,w(ixO^S,nw+1))
}{#IFNDEF STAGGERED
! Reduce output array size, +1 was added for eventual pointdata output
call get_divb(ixI^L,ixO^L^LSUB1,w(ixI^S,1:nw),w(ixI^S,nw+1))
}
end associate
end subroutine specialvar_output
!=============================================================================
subroutine specialvarnames_output
! newly added variables need to be concatenated with the varnames/primnames string
include 'amrvacdef.f'
!-----------------------------------------------------------------------------
primnames= TRIM(primnames)//' '//'divB'
wnames=TRIM(wnames)//' '//'divB'
end subroutine specialvarnames_output
!=============================================================================
To output this (one) extra variable, be sure to set nwauxio = 1 in the filelist section of the parameter file.