ProblemGridInit

This routine is responsible for initializing/modifying cell centered quantities with in each Info structure. Given an initial function f(x) that specifies the state of the gas at each position, this routine could be as simple as

DO i=1,Info%mX(1)
  DO j=1,Info%mX(2)
    DO k=1,Info%mX(3)
      pos=Info%xBounds(:,1)+((/i,j,k/)-.5)*levels(Info%level)%dx
      Info%q(i,j,k,:)=f(pos)
    END DO
  END DO
END DO

We don't need to initialize values in ghost regions of grids since internal overlaps and physical boundary conditions will supply the required values.

Often the initial conditions are defined piecewise over various physical regions. For example multiple clumps would each have a physical location in which they require the density to be modified. For this reason it is useful to first calculate the overlap bounds between an Info structure and a physical object. Because of periodic boundary conditions - it is possible for a physical object to have several overlap regions. Since this is a common thing to need - there is a routine for doing so. Let's say that there is a clump whose physical extents are limited by a cube whose vertices are defined in an array xbound(3,2).

 CALL CalcPhysicalOverlaps(Info, xBounds, mSs, nOverlaps, offsets)
 DO n=1,nOverlaps
    mS=mSs(n,:,:)
    offset=offsets(n,:)
    DO k=mS(3,1),mS(3,2)
       pos(3)=Info%xBounds(3,1)+offset(3)+(k-1)*dx
       DO j=mS(2,1),mS(2,2)
          pos(2)=Info%xBounds(2,1)+offset(2)+(j-1)*dx
          DO i=mS(1,1),mS(1,2)
             pos(1)=Info%xBounds(1,1)+offset(1)+(i-1)*dx
             Info%q(i,j,k,:)=f(pos)
          END DO
       END DO
    END DO
 END DO

Additionally since the cell centered value represents the volume averaged value, it is sometimes convenient to sub-sample the initial conditions within cells. This can be accomplished as follows:

 sample_res=1
 sample_res(1:nDim)=SubSample                 
 sample_fact=0
 sample_fact=1d0/REAL(sample_res)
 q_fact=product(sample_fact(1:nDim))
 CALL CalcPhysicalOverlaps(Info, xBounds, mSs, nOverlaps, offsets)
 DO n=1,nOverlaps
    mS=mSs(n,:,:)
    offset=offsets(n,:)
    DO k=mS(3,1),mS(3,2)
       xpos(3)=Info%xBounds(3,1)+offset(3)+(k-1)*dx
       DO j=mS(2,1),mS(2,2)
          xpos(2)=Info%xBounds(2,1)+offset(2)+(j-1)*dx
          DO i=mS(1,1),mS(1,2)
             xpos(1)=Info%xBounds(1,1)+offset(1)+(i-1)*dx
             q_Source=0
             DO kk=1,sample_res(3)
                pos(3)=xpos(3)+(REAL(kk, 8)-half)*dx*sample_fact(3)
                DO jj=1,sample_res(2)
                   pos(2)=xpos(2)+(REAL(jj, 8)-half)*dx*sample_fact(2)
                   DO ii=1,sample_res(1)
                      pos(1)=xpos(1)+(REAL(ii, 8)-half)*dx*sample_fact(1)
                      q_source=q_source+f(pos)-Info%q(i,j,k,1:NrHydroVars))
                   END DO
                END DO
             END DO
             Info%q(i,j,k,1:NrHydroVars)=Info%q(i,j,k,1:NrHydroVars)+q_source*q_fact
          END DO
       END DO
    END DO
 END DO

(Of course if you can numerically integrate the density function then sub-sampling would not be necessary - however this can be difficult especially for functions that are defined piecewise over spherical regions etc…)

If the source function is just a constant value over a geometric shape - then you can avoid writing all of this code and just use the UniformRegions object. If you create a UniformRegions object in your ProblemModuleInit routine and specify the shape and orientation as well as the constant value for Info%q then you don't have to every worry about the Info structures.

Initializing Aux Fields

In AstroBear, the B-fields are stored at cell faces (or edges in 2D). This makes it easier to ensure that a divergence free magnetic field remains divergence free. Of course this only works if the initial conditions are divergence free. To further complicate things - discretization errors can cause a divergence free magnetic field solution to have residual divergence due to discretization errors. To see how let us consider a divergence free B-field.

\mathbf{B}=[cos(y)cos(x), -sin(x)sin(y)]$

Now let us evaluate the numerical divergence of a cell centered at (y_0, x_0)

begin{equation} \nabla_c \cdot \mathbf{B} = \frac{1}{\Delta} \left [ cos(y_0)[cos(x_0+\Delta/2)-cos(x_0-\Delta/2)] - sin(x_0)[sin(y_0+\Delta/2)] - sin(y_0-\Delta/2)] \right ] \end{equation} \begin{equation} =\frac{1}{\Delta} \left [ cos(y_0)[-2sin(x_0)sin(\Delta/2)] - sin(x_0)[2cos(y_0)sin(\Delta/2)] \right ]=-\frac{4}{\Delta}cos(y_0)sin(x_0)sin(\Delta/2) \end{equation}

If \Delta is small compared to the wavelength, then this will be small - but still non-zero.

The reason is that the face centered fields should store the face-averaged value of the B-field. For functions that vary across a cell face, the value at the center will not equal the average value over the cell face. Note however, that if we had used a simpler form of the potential \mathbf{B}=[f(y), g(x)], then the value of the B-field at face centers would equal the value averaged over the cell face, and therefore the numerical divergence would also be zero by Stoke's theorem.

Alternatives in 2D

In 2D, the magnetic fields are stored on cell edges and the edge-centered value equal the average field along the cell edge. This requires integration - either numerically or algebraically. Numerical integration is essentially sub-sampling thereby reducing the effective \Delta. Algebraic integration is identical to solving for the vector potential, and then discretizing it at cell corners before numerically taking the curl. \nabla_c \cdot ( \nabla_c \times A)=0.

 IF (nDim == 2) THEN
   ALLOCATE(emf(1:Info%mX(1)+1, 1:Info%mX(2)+1))
   emf=0
   DO i=1,Info%mX(1)+1
     xpos(1)=Info%xBounds(1,1)+offset(1)+(i-1)*dx
     DO j=Info%mX(2)+1
       xpos(2)=Info%xBounds(2,1)+offset(2)+(j-1)*dx
       emf(i,j)=A_z(xpos)
     END DO
   END DO
   CALL AddCurl(Info%aux(1:Info%mX(1)+1, 1:Info%mX(2)+1,1,1:2),emf(:,:),dx)
   DEALLOCATE(emf)
 END IF

Alternatives in 3D

In 3D, the magnetic fields are stored on cell faces and the face-centered value should equal the average value over the entire face. This requires a 2D integration. Calculating the vector potential would take care of one of these integrations and assuming that it is done along cell edges can be used to calculate a numerically divergence free field - however it will not necessarily equal the true face-averaged magnetic field. This may only cause problems at coarse fine boundaries, where the discretization errors due to grid sizes may lead to different values for the magnetic field (even though the magnetic fields are divergence free within a given AMR level). To ensure that both coarse and fine grids have the same discretization error - it is only necessary to sub-sample the coarse grids along cell edges at the same resolution as the finest grids.

  IF (nDim == 3) THEN
   sample_res=2**(MaxLevel-Info%level)
   sample_fact=1d0/REAL(sample_res)
   ALLOCATE(emf(mS(1,1):mS(1,2)+1, mS(2,1):mS(2,2)+1,mS(3,1):mS(3,2)+1,3), emf_source(3))
   emf=0
   DO i=mS(1,1),mS(1,2)
     xpos(1)=Info%xBounds(1,1)+offset(1)+(i-1)*dx
     DO j=mS(2,1),mS(2,2)
       xpos(2)=Info%xBounds(2,1)+offset(2)+(j-1)*dx
       DO k=mS(3,1),mS(3,2)
         xpos(3)=Info%xBounds(3,1)+offset(3)+(k-1)*dx
         DO m=1,3
           pos=xpos
           DO ii=1,sample_res
             pos(m)=xpos(m)+(REAL(ii, 8)-half)*dx*sample_fact
             emf_source=emf(pos)
             emf(i,j,k,m)=emf(i,j,k,m)+emf_source(m)
           END DO
           emf(i,j,k,m)=emf(i,j,k,m)*sample_fact
         END DO
       END DO
     END DO
   END DO
   CALL AddCurl(Info%aux(1:Info%mX(1)+1, 1:Info%mX(2)+1, 1:Info%mX(3)+1,1:3),emf,dx)
   DEALLOCATE(emf, emf_source)
 END IF

Finally, if an emf is applied to only part of the domain - it must physically be zero everywhere it is not applied. For example if there is a localized B-field which can be expressed as the curl of a potential - the potential must go to zero at the surface. Since you can add a gradient to a potential without effecting it's curl, this can always (and often is required) to be done.

After initializing the aux fields it is often necessary to calculate the cell-centered magnetic energy in order to update the cell-averaged total energy. During the simulation, the cell-centered B-field is always calculated as the average of the face-centered B-fields - so it is a good idea to skip the calculation of the the volume averaged B-field and instead replace it with the average of the face-averaged B-fields. After doing this you can then update the total energy with the volume-averaged thermal and kinetic energies plus the magnetic energy from the face-averaged fields.