This document describes how the problem dependent amrvacusr.t and the optional amrvacusrpar.t files should be written for user defined initial-, boundary conditions, source terms and special variables for IO.
There is additional functionality like controling AMR and (de-) activating parts of the mesh and we refer to the dummy subroutines or implemented tests for reference on the more advanced features.
This page: [General] [Creating] [Initial conditions] [Boundary conditions] [Source treatments] [Custom variables for IO]


The AMRVACUSR modules contain the problem dependent user-written subroutines.

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.

Creating a New Setup

To start fresh, in your designated simulation directory, copy the src/amrvacnul/amrvacusr.t file. It consists of a few include statements and three subroutines needed to specify the initial conditions of your problem. The arguments are declared and the purpose of the subroutines is described below. Comment out the INCLUDE:amrvacnul/XXX.t statement(s) for additional subroutine(s) that you intend to implement (copy in the defaults from src/amrvacnul/XXX.t to get started).

Initial conditions

Your start file should look something like this, which is actually the src/amrvacnul/amrvacusr.t default file:
  ! amrvacusr.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:
!      call b_from_vectorpotential(ixI^L,ixO^L,w,x)
!      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.

Boundary conditions

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).

Specialsource part

There are lots of possible physical source terms for the same basic equation. Rather than writing a new physics module for each, it is simpler to define a problem dependent source term in the AMRVACUSR module. The specialsource subroutine is called at least once in every time step by BHAC. The number of calls depends on the time integration scheme defined by typeadvance, the parameter sourcesplit and also on the number of dimensions if the parameter dimsplit=T.

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


  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.

Custiom variables for IO

The amrvacnul/speciallog.t file contains additional subroutines related to special I/O requests.

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})

      call div_staggered(ixO^L,s,w(ixO^S,nw+1))
      ! 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.