Sub Sampling

Cell centered quantities

• AstroBEAR stores most fluid quantities as volume averages in the array Info%q(i,j,k,:)
• Each cell is a cube with a side dx=levels(Info%level)%dx
• The 'lower' corner of the cube is located at Info%xbounds(:,1) + ((/i,j,k/)-1)*dx
• The 'upper' corner of the cube is located at Info%xbounds(:,1) + ((/i,j,k/))*dx
• The 'center' of the cube is located at Info%xbounds(:,1) + (REAL(/i,j,k/)-.5)*dx

Sub-sampling

• If we break each cube into NxNxN sub-pieces, and set ddx=dx/N then the center of the [ii,jj,kk] sub piece will be at Info%xBounds(:,1)+(/i-1,j-1,k-1/)*dx + (REAL(/ii,jj,kk/)-.5)*ddx
• So we could use the following routine to initialize cell centered quantities in q. Here is a routine for density.

  dx=levels(Info%level)%dx
  ddx=dx/N
  DO i=1,Info%mX(1)
    DO j=1,Info%mx(2)
      DO k=1,Info%mx(3)
        new_rho=0d0
        DO ii=1,N
          DO jj=1,N
            DO kk=1,N
              pos=Info%xBounds(:,1)+(/i-1,j-1,k-1/)*dx + (REAL(/ii,jj,kk/)-.5)*ddx
              new_rho=new_rho+rho_func(pos)
            END DO
          END DO
        END DO
        Info%q(i,j,k,1)=new_rho/N**nDim
      END DO
    END DO
  END DO
  

or we could try to speed it up by reducing some repeated computations.

dx=levels(Info%level)%dx
ddx=dx/N
DO i=1,Info%mX(1)
  pos(1)=Info%xBbounds(1,1)+(i-1)*dx
  DO j=1,Info%mx(2)
    pos(2)=Info%xBbounds(1,1)+(i-1)*dx
    DO k=1,Info%mx(3)
      pos(3)=Info%xBbounds(1,1)+(i-1)*dx
      new_rho=0d0
      DO ii=1,N
        ppos(1)=pos(1)+(REAL(ii)-.5)*ddx
        DO jj=1,N
          ppos(2)=pos(2)+(REAL(jj)-.5)*ddx
          DO kk=1,N
            ppos(3)=pos(3)+(REAL(kk)-.5)*ddx
            new_rho=new_rho+rho_func(pos)
          END DO
        END DO
      END DO
      Info%q(i,j,k,1)=new_rho/N**nDim
    END DO
  END DO
END DO

and sometimes we may be implementing a function that is only defined for a certain region (like within a certain radius etc…) where the value for q outside of a given radius is unknown and should be left unchanged. In that case, we just want to add up changes to q for subcells within the region.

dx=levels(Info%level)%dx
ddx=dx/N
DO i=1,Info%mX(1)
  pos(1)=Info%xBbounds(1,1)+(i-1)*dx
  DO j=1,Info%mx(2)
    pos(2)=Info%xBbounds(1,1)+(i-1)*dx
    DO k=1,Info%mx(3)
      pos(3)=Info%xBbounds(1,1)+(i-1)*dx
      drho=0d0
      DO ii=1,N
        ppos(1)=pos(1)+(REAL(ii)-.5)*ddx
        DO jj=1,N
          ppos(2)=pos(2)+(REAL(jj)-.5)*ddx
          DO kk=1,N
            ppos(3)=pos(3)+(REAL(kk)-.5)*ddx
            IF (sum(ppos(:)**2) > r2) THEN
              drho=drho+(rho_func(ppos)-Info%q(i,j,k,1))
            END IF
          END DO
        END DO
      END DO
      Info%q(i,j,k,1)=Info%q(i,j,k,1)+drho/N**nDim
    END DO
  END DO
END DO

Face centered quantities

• AstroBEAR stores magnetic fields as face averages in the array Info%aux(i,j,k,:)
• Each cell is a face perpendicular to the component of the B-field with sides of length dx=levels(Info%level)%dx
• What is more important is the location of the edges since those will store in temporary Vector Potentials which can then be used to update the aux fields.
• The 'lower' end point of every component of the vector potential A(i,j,k,m) is at Info%xbounds(:,1) + ((/i,j,k/)-1)*dx
• The 'upper' end point of each component of the vector potential A(i,j,k,m) is at Info%xbounds(:,1) + ((/i,j,k/)+I(m,:))*dx where I(:,:) is the identity matrix so it effectively adds 1 along the 'm' coordinate.

Sub-sampling

For face centered fields we could also subsample, but there is the divergence criterion to consider for B-fields. We could calculate the potential on a subgrid, and then take a bunch of curls on each face and add them up, but Stoke's theorem lets us just do the integral around the outside. So we really just need to subsample along the appropriate edge for each component of the vector potential. So for example to calculate the x component of the potential along a given edge we would

dx=levels(Info%level)%dx
ddx=dx/N
DO i=1,Info%mX(1)
  pos(1)=Info%xBbounds(1,1)+(i-1)*dx
  DO j=1,Info%mx(2)+1
    pos(2)=Info%xBbounds(1,1)+(i-1)*dx
    DO k=1,Info%mx(3)+1
      pos(3)=Info%xBbounds(1,1)+(i-1)*dx      
      A(i,j,k,1)=0d0
      DO ii=1,N
        ppos(1)=pos(1)+(REAL(ii)-.5)*ddx
        A(i,j,k,1)=Ax(ppos)
      END DO
      A(i,j,k,1)=A(i,j,k,1)/N
    END DO
  END DO
END DO

or we could do all three components of the potential like this

dx=levels(Info%level)%dx
ddx=dx/N
DO m=1,3
  p=1
  p(m)=0
  DO i=1,Info%mX(1)+p(1)
    pos(1)=Info%xBbounds(1,1)+(i-1)*dx
    DO j=1,Info%mx(2)+p(2)
      pos(2)=Info%xBbounds(1,1)+(i-1)*dx
      DO k=1,Info%mx(3)+p(3)
        pos(3)=Info%xBbounds(1,1)+(i-1)*dx
        A(i,j,k,m)=0d0
        DO ii=1,N
          ppos(m)=pos(m)+(REAL(ii)-.5)*ddx
          temp(:)=VectorPotential(ppos) !This returns a 3-component vector for simplicity
          A(i,j,k,m)=A(i,j,k,m)+temp(m) !but we only retain the 'm' component.
        END DO
        A(i,j,k,m)=A(i,j,k,m)/N
      END DO
    END DO
  END DO
END DO

There are many examples of loops like these in the various objects (ie clumps, outflows, disks, etc..)

One could imagine writing functions like 'volume average', or 'edge average' that could replace those inner loops - that would require a user function to be passed.