The recommended place to set-up initial data for BHAC is the routine initonegrid_usr() in the user file (amrvacusr.t). However, in some cases, it may be preferable to read the initial data for the simulation directly from a table. This may be the case for data-driven simulations, or for situations when initial data is generated by another program and the effort to implement natively-generated data is not justified.
To facilitate the use of numerical initial data, BHAC possess the module mod_oneblock, which contains routines to read initial data in ASCII format and perform interpolations. The initial data must be of the same dimension as the problem set-up; however, there are ways to circumvent this.
It is not necessary that the tabulated data is uniformly spaced, but it needs to be structured in a grid. It is important to mention that it does not need to have the same resolution as the BHAC simulation, although it is recommended that the resolution is approximately similar or higher in order to avoid interpolation artifacts. The domain size does not need to coincide either, but in case the domain of tha BHAC simulation is larger, flat interpolation will be applied.
Input format
The ASCII file containing the data needs to be in a pre-defined format, which is the same as that of BHAC’s ‘oneblock’ output. This is:
list of coordinates and variable names, separated by spaces
Ntotal N1 [N2] [N3]
time
coord1_value1 [coord2_value2] [coord3_value3] variable_values
...Where N1, N2, and N3 are the numbers of cells along dimensions 1, 2 and 3, and Ntotal = N1xN2xN3. The elements in square brackets are optional. Here is an example of a one-dimensional initial data file:
r rho u1 u2 u3 p b1 b2 b3 lfac xi
938 938
0.000000
1.023826E-04 1.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 2.004809E-04 0.000000E+00 0.000000E+00 0.000000E+00 1.000000E+00 1.000000E+00
1.289948E-04 9.999608E-01 0.000000E+00 0.000000E+00 0.000000E+00 2.004701E-04 0.000000E+00 0.000000E+00 0.000000E+00 1.000000E+00 1.000000E+00
...We can call this file, for instance, idata.blk.
Reading initial data into BHAC
To use the tabulated initial data, it is first necessary to load the file. This can be done by adding a call to read_oneblock(), which accepts as argument the name of the initial data file. A good place to do it is at the end of initglobaldata_usr(), in amrvacusr.t. For example,
subroutine initglobaldata_usr
use mod_oneblock
! Other modules and variable declarations
...
!-----------------------------------
! Subroutine contents
...
call read_oneblock('idata.blk')
end subroutine initglobaldata_usrOnce the data has been read, we can use the routine interpolate_oneblock() to fill the initial data in our problem set-up where required. The recommended place to do this is the routine initonegrid_usr() in amrvacusr.t. For instance, one can assign the initial data by interpolating point by point using a loop.
One can also modify the initial data within the same routine as desired. In the following example, density and pressure are read from the tabulated initial data, while magnetic field and velocity components are set to zero. Later, regions of low pressure and density are replaced by floors. The calls to the interpolation routine are preceded by a coordinate transformation (not shown), to highlight the possibility that the tabulated initial data could be in coordinates different from those of the simulation. It is important to recall that, the variables must exit initonegrid_usr() in a conserved state, hence we explicitly show the call to conserve().
subroutine initonegrid_usr(...)
use mod_oneblock
! Other modules and variable declarations
...
!----------------------------------------
...
! Coordinate transformation from code coordinates 'x'
! to Boyer-Lindquist coordinate 'xBL
...
do r_idx=ixmin1,ixmax1
call interpolate_oneblock(xBL_point,rho_,w(r_idx,rho_))
call interpolate_oneblock(xBL_point,pp_,w(r_idx,pp_))
end do
w(ixG^S,u1_) = 0.0d0
w(ixG^S,u2_) = 0.0d0
w(ixG^S,u3_) = 0.0d0
w(ixG^S,b1_) = 0.0d0
w(ixG^S,b2_) = 0.0d0
w(ixG^S,b3_) = 0.0d0
where(w(ix^S,rho_) .lt. eqpar(rhomin_))
w(ix^S,rho_)=eqpar(rhomin_)
end where
where(w(ix^S,pp_) .lt. eqpar(pmin_))
w(ix^S,pp_)=eqpar(pmin_)
end where
call conserve(ixG^L,ix^L,w,x,patchfalse)
...
end subroutine initonegrid_usr
With these small additions, we can compile BHAC and it will be able to start a simulation from our tabulated initial data.