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.
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 ...
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,
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
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_usr
Once 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
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
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.