| Version 1 (modified by , 8 years ago) ( diff ) |
|---|
Simulation Data
All AstroBEAR modules have at least one thing in common: initializing the problem domain. Within our code, the problem domain's data is held in InfoDef structures, which is why so many module subroutines take an InfoDef structure as a parameter. To make use of the InfoDef structure, the following statement is required at the beginning of problem.f90:
USE DataDeclarations
Anything that is defined in the InfoDef type is now available. For example, q need not be defined…just reference it by Info%q.
There are two major data arrays in InfoDef: the q array and the aux array. q holds the volume averaged data and is used by all AstroBEAR simulations while aux holds face averaged data and is used only for MHD when nDim > 1. Note that volume averaged data and cell-centered data are often used interchangeably, but there is an important distinction. Take for instance a simple function defined on the interval . Applying a Taylor expansion to and finding the average of gives
where the cell centered value is
so the cell centered value is second order accurate for the volume average and usually is a quick way to estimate the volume average. However if the function has large 2nd derivatives (or higher) this can lead to large errors in the volume average. This is often apparent when modeling discontinuities along curved boundaries. There are two ways to handle this problem:
- Introduce smoothing to the physical model to remove large 2nd and higher derivatives
- Better approximate the volume average either by
- Analytical integration (often non-trivial if possible)
- Numerical integration (ie SubSampling)
The q array takes the form q(x,y,z,variable) where variable is an index that refers to the various physical quantities such as density, momentum, energy, etc. in each cell. The order of the quantities in the variable array is dependent on the equation of state, whether or not magnetic fields are being tracked, etc… For 2D hydro (non MHD) the order of the fields is (rho, px, py, E). So if we wanted to set the energy of the cell at integer location i,j,k we would use
Info%q(i,j,k,4) = 1.0
However if we were to change the number of dimensions from 2 to 3, then the order of the fields would be rho, px, py, pz, E and the above statement would not set the energy, but the z momentum to 1.0 and leave the energy unchanged. The solution is to avoid using integer constants for the 4th array index and instead use integer variables that are adjusted based on the equations of state, number of dimensions, etc… These variables are declared in PhysicsDeclarations so we need to also add
USE PhysicsDeclarations
to the top of our module. Then we can set the energy of cell i,j,k regardless of the specifics of our problem by using
Info%q(i,j,k,iE) = 1.0
Also, if we happen to be using an isothermal equation of state, then the energy is no longer stored within the q array and the value of iE is set to 0 to indicate this. So it is generally a good idea to check the value of iE as follows
IF (iE /= 0) Info%q(i,j,k,iE)=1.0
Additional variables used to store slots are:
- irho - density (always non-zero)
- ivx - x momentum (always non-zero)
- ivy - y momentum (non-zero unless 2D, 3D, MHD)
- ivz - z momentum (non-zero unless 3D, MHD)
- iE - Energy (non-zero unless isothermal EOS)
- iBx - x magnetic field (non-zero unless MHD)
- iBy - y magnetic field (non-zero unless MHD)
- iBz - z magnetic field (non-zero unless MHD)
There are also two arrays that are sometimes useful as well
- iB(1:3) = (/iBx, iBy, iBz/)
- imom(1:3) = (/ivx, ivy, ivz/)
The aux array holds face-centered data, and is only used in MHD problems. If you are running a strictly hydrodynamic problem or a hydrodynamic + elliptic problem, then you will not need aux.
Dimensions
The number of cells in the x, y, & z direction for the core region of each Info structure is stored in the array
Info%mX(1:3)
and often one will declare local variables mx, my, & mz to avoid repeatedly having to type Info%mx(d).
mx=Info%mX(1) my=Info%mX(2) mz=Info%mX(3)
The data within this core region (which does not include ghost zones) is stored in Info%q(1:mx,1:my,1:mz,1:NrHydroVars) where NrHydroVars represents the number of fluid variables including tracers. If running with fewer than 3 dimensions, the unused dimensions have an extent of 1.
Before we can initialize a cell we must calculate it's physical location and extent. To do so we need to know the cell size for the Info's AMR level. The properties of each level are stored in the levels(:) array. To access this data we must use the GlobalDeclarations module by adding the following to our module at the top.
USE GlobalDeclarations
Then to access properties of level n - for example the current time that level has advanced to we would use levels(n)%tnow. If we wanted to now the current time step for level n we could use levels(n)%dt. And to access the cell size for level n we could use levels(n)%dx. Since the level a given info structure resides on is stored in Info%level, the cell size is given by levels(Info%level)%dx. So to get the x-position of the center of a cell with x-index i we could use
xlower=Info%xBounds(1,1) dx=levels(Info%level)%dx x=xlower+(REAL(i)-.5)*dx
Note we subtract 0.5 from the index before multiplying by the spacing since we are calculating the cell center. And that the cell actually goes from x-.5*dx to x+.5*dx. Also note that we convert the integer to a real before subtracting .5. And if we want to calculate x,y,z we could use
xlower=Info%xBounds(1,1) dx=levels(Info%level)%dx x=xlower + (REAL(i)-.5)*dx y=ylower + (REAL(j)-.5) * dx z=zlower + (REAL(k)-.5) * dx IF (nDim < 2) y=ylower IF (nDim < 3) z=zlower
The last two lines are necessary since we don't want to add .5 to the y or z dimensions if we are only in 1D or 2D. We could also streamline this using the Fortran MERGE function and storing (/x,y,z/) in an array pos(:) using
pos=Info%xBounds(:,1)+merge((REAL((/i,j,k/))-.5)*dx, (/0d0,0d0,0d0/), nDim < (/1,2,3/))
Finally since the precision of the various info fields related to spatial position is a parameter xPrec (could be single or double), some compilers will complain unless you convert (/i,j,k/ as well as .5 to the right kind of REAL.
pos=Info%xBounds(:,1)+merge((REAL((/i,j,k/),KIND=xPREC)-half)*dx, (/0d0,0d0,0d0/), nDim < (/1,2,3/))
Note that the variable half is a parameter equal to REAL(.5, KIND=xPREC) declared in GlobalDeclarations
Finally there is a function already called CellPos that does the same calculation which makes life much easier.
pos=CellPos(Info, i, j, k)
The Info%aux array is a little different. The aux array holds magnetic flux values, which are face-averaged. This means that every volume averaged value in Info%q is bracketed in each dimension by two Info%aux values. To accommodate the extra values, Info%aux is a 1:mx+1 by 1:my+1 by 1:mz+1 box, but the aux dimensions are actually different for each variable:
Bx = Info%aux(1:mx+1, 1:my, 1:mz, 1) By = Info%aux(1:mx, 1:my+1, 1:mz, 2) Bz = Info%aux(1:mx, 1:my, 1:mz+1, 3)
The additional cells (the ones in the "upper-front right" corner of the aux array) are not used. To locate the center of the face for the Bx fields, we would subtract half*dx from the cell center.
x_pos=CellPos(Info, i, j, k)-(/half,0d0,0d0/)
and for By and Bz we could use
y_pos=CellPos(Info, i, j, k)-(/0d0,half,0d0/) z_pos=CellPos(Info, i, j, k)-(/0d0,0d0,half/)
