64 | | There are two major data arrays in {{{InfoDef}}}: the {{{q}}} array and the {{{aux}}} array. {{{q}}} holds the cell-centered data and is used by all AstroBEAR simulations. 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 |
| 64 | 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 f(x) defined on the interval [0:h]. The average of f(x) over the interval is |
| 65 | |
| 66 | f(0)+d,,1,,f(0)h/2+d,,2,,f(0)h^2^/3+... |
| 67 | |
| 68 | where the cell centered value is |
| 69 | |
| 70 | f(0)+d,,1,,f(0)h/2+d,,2,,f(0)h^2^/6+... |
| 71 | |
| 72 | 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: |
| 73 | * Introduce smoothing to the physical model to remove large 2nd and higher derivatives |
| 74 | * Better approximate the volume average either by |
| 75 | * Analytical integration (often non-trivial if possible) |
| 76 | * Numerical integration (ie sub-sampling) |
| 77 | |
| 78 | 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 |
100 | | |
101 | | The ''core'' region of the {{{Info%q}}} array (which does not include ghost zones) is {{{Info%q(1:Info%mx(1),1:Info%mx(2),1:Info%mx(3),1:NrHydroVars)}}} where {{{Info%mx(1)}}}, {{{Info%mx(2)}}}, and {{{Info%mx(3)}}} denote the number of cells in the x, y, and z directions, respectively and {{{NrHydroVars}}} represents the number of fluid variables including tracers. If running with fewer than 3 dimensions, the unused dimensions have an extent of 1. |
102 | | |
103 | | 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. |
104 | | |
105 | | {{{Info%q}}} is cell-centered, so the values are assumed to be taken from the midpoint of the cell. The lower bounds of the info patch for the x, y, and z direction are given by the components of the array {{{Info%xBounds(1:3,1)}}}. And the cell size is given by {{{levels(Info%level)%dx}}}. Thus, the index-to-space conversion is: |
106 | | {{{ |
107 | | x=Info%xBounds(1,1) + (REAL(i,xPrec)-half) * levels(Info%level)%dx) |
108 | | y=(ylower + (REAL(j,xPrec)-half) * dx) |
109 | | z=(zlower + (REAL(k,xPrec)-half) * dx) |
110 | | |
111 | | }}} |
112 | | |
113 | | Note that {{{i, j, k, x, y, z, xlower, ylower, zlower,}}} and {{{dx}}} are user defined variables and are explained in further detail below. The variable {{{half}}} is used since the data is in the center of the cell, and {{{xPrec}}} is a type of precision. {{{half}}} and {{{xPrec}}} are already defined elsewhere and can be used if the following statement is at the beginning of {{{problem.f90}}}: |
| 114 | The number of cells in the x, y, & z direction for the '''core''' region of each Info structure is stored in the array |
| 115 | {{{ |
| 116 | Info%mX(1:3) |
| 117 | }}} |
| 118 | and often one will declare local variables {{{mx, my, & mz}}} to avoid repeatedly having to type Info%mx(d). |
| 119 | {{{ |
| 120 | mx=Info%mX(1) |
| 121 | my=Info%mX(2) |
| 122 | mz=Info%mX(3) |
| 123 | }}} |
| 124 | |
| 125 | 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. |
| 126 | |
| 127 | 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. |
117 | | |
118 | | [[BR]] |
119 | | The {{{Info%aux}}} array is a little different. The {{{aux}}} array holds magnetic flux values, which are face-centered. This means that every cell-centered value in {{{Info%q}}} is bracked 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: |
| 131 | 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 |
| 132 | {{{ |
| 133 | xlower=Info%xBounds(1,1) |
| 134 | dx=levels(Info%level)%dx |
| 135 | x=xlower+(REAL(i)-.5)*dx |
| 136 | }}} |
| 137 | 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 |
| 138 | {{{ |
| 139 | xlower=Info%xBounds(1,1) |
| 140 | dx=levels(Info%level)%dx |
| 141 | x=xlower + (REAL(i)-.5)*dx |
| 142 | y=ylower + (REAL(j)-.5) * dx |
| 143 | z=zlower + (REAL(k)-.5) * dx |
| 144 | IF (nDim < 2) y=ylower |
| 145 | IF (nDim < 3) z=zlower |
| 146 | }}} |
| 147 | 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 |
| 148 | {{{ |
| 149 | pos=Info%xBounds(:,1)+merge((REAL((/i,j,k/))-.5)*dx, (/0d0,0d0,0d0/), nDim < (/1,2,3/)) |
| 150 | }}} |
| 151 | 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. |
| 152 | {{{ |
| 153 | pos=Info%xBounds(:,1)+merge((REAL((/i,j,k/),KIND=xPREC)-half)*dx, (/0d0,0d0,0d0/), nDim < (/1,2,3/)) |
| 154 | }}} |
| 155 | Note that the variable {{{half}}} is a parameter equal to {{{REAL(.5, KIND=xPREC)}}} declared in GlobalDeclarations |
| 156 | |
| 157 | Finally there is a function already called !CellPos that does the same calculation which makes life much easier. |
| 158 | {{{ |
| 159 | pos=CellPos(Info, i, j, k) |
| 160 | }}} |
| 161 | |
| 162 | |
| 163 | [[BR]] |
| 164 | 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: |
127 | | The additional cells (the ones in the "upper-front right" corner of the {{{aux}}} array) are not used. |
128 | | |
| 172 | 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. |
| 173 | {{{ |
| 174 | x_pos=CellPos(Info, i, j, k)-(/half,0d0,0d0/) |
| 175 | }}} |
| 176 | and for By and Bz we could use |
| 177 | {{{ |
| 178 | y_pos=CellPos(Info, i, j, k)-(/0d0,half,0d0/) |
| 179 | z_pos=CellPos(Info, i, j, k)-(/0d0,0d0,half/) |
| 180 | }}} |
133 | | Astrophysical problems involve many different physical units and constants with a wide range of scales. To avoid the loss of precision that comes when computers try to work with, say, a 10^-8^ variable and a 10^24^ constant in the same expression, we scale our units into ''computational units'' before storing them in the data arrays. |
134 | | |
135 | | Usually, the physical scales are defined in the [PhysicsDataExplained physics.data] file--you simply enter the scales for density, temperature, velocity, etc in that file, and AstroBEAR will read them in. More complicated scaling would be defined in the {{{ProblemModuleInit()}}} routine (see above). |
| 185 | Astrophysical problems involve many different physical units and constants with a wide range of scales. To avoid overflow or underflow - we scale our units into ''computational units'' before storing them in the data arrays. Note with double precision this would be quite rare - but it still convenient to work within physical units appropriate to the problem. |
| 186 | |
| 187 | Usually, the physical scales are defined in the [PhysicsDataExplained physics.data] file --you simply enter the scales for density, temperature, velocity, etc in cgs units, and AstroBEAR will read them in. ''Note that nScale is in cm^{-3}^ and !TempScale is in Kelvin'' |