MHD

To enable MHD:

• In physics.data, set lMHD to T.

If the code is using constrained transport, then the logical flag MaintainAuxArrays will be .True.. In that case you will need to initialize the face averaged magnetic field values in the aux array in your module file. The following diagram displays the various aux variables and what positions they correspond to - as well as the positions of the EMF as well as the vector potential.

Here is an example for an initial magnetic field that is \mathbf {B}=Bx(x,y,z)\mathbf{e_x} + By(x,y,z) \mathbf{e_y} + Bz(x,y,z) \mathbf{e_z}

dx = (/(merge(1d0,0d0,nDim >= i),i=1,3)/)*levels(Info%level)%dx

DO i=1, Info%mX(1)
  x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
  DO j=1, Info%mX(2)
    y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
    DO k=1, Info%mX(3)
      z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
      Info%q(i,j,k,iBx)=Bx(x,y,z)
      Info%q(i,j,k,iBy)=By(x,y,z)
      Info%q(i,j,k,iBz)=Bz(x,y,z)
    END DO
  END DO
END DO

IF (MaintainAuxArrays) THEN
  IF (nDim >= 1) THEN
    !Note the +1 in the i loop DO statement and the -1d0 in 
    !the x coordinate calculation instead of -.5d0
    DO i=1, Info%mX(1)+1
      x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
      DO j=1, Info%mX(2)
        y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
        DO k=1, Info%mX(3)
          z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
          Info%aux(i,j,k,1)=Bx(x,y,z)
        END DO
      END DO
    END DO

    !Now we can update the cell centered value with the adjacent face average
    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBx)=.5d0 * (&
      Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1) + &
      Info%aux(2:Info%mX(1)+1,1:Info%mX(2),1:Info%mX(3),1) )


    IF (nDim >= 2) THEN
      !Note the +1 in the j loop DO statement and the -1d0 in 
      !the y coordinate calculation instead of -.5d0
      DO i=1, Info%mX(1)
        x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
        DO j=1, Info%mX(2)+1
          y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
          DO k=1, Info%mX(3)
            z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
            Info%aux(i,j,k,2)=By(x,y,z)
          END DO
        END DO
      END DO

      !Now we can update the cell centered value with the adjacent face average  
      Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBy)=.5d0 * (&
        Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),2) + &
        Info%aux(1:Info%mX(1),2:Info%mX(2)+1,1:Info%mX(3),2) )

      IF (nDim >= 3) THEN
       !Note the +1 in the k loop DO statement and the -1d0 in 
       !the z coordinate calculation instead of -.5d0
       DO i=1, Info%mX(1)
         x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
         DO j=1, Info%mX(2)
           y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
           DO k=1, Info%mX(3)+1
             z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
             Info%aux(i,j,k,3)=Bz(x,y,z)
           END DO
         END DO
       END DO

      !Now we can update the cell centered value with the adjacent face average
       Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBz)=.5d0 * (&
         Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),3) + &
         Info%aux(1:Info%mX(1),1:Info%mX(2),2:Info%mX(3)+1,3) )
      END IF
    END IF
  END IF
END IF

Of course the above example works for any field - but does not impose the divergence free constraint. Presumably the functions Bx, By, and Bz would take care of that to the degree they are accurate face averages - but a better approach is to instead approximate the vector potential at cell edges and then discretely difference it to get the face average fields.

IF (nDim == 2) THEN
  !In 2D the only component of the vector potential that matters is Az
  ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1,1))
  DO i=1, Info%mX(1)+1
    x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    DO j=1, Info%mX(2)+1
      y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
      A(i,j,1,1)=Az(x,y,z)
    END DO
  END DO

  !Then take the curl of A to get B
  FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2))
    Info%aux(i,j,1,1)=+(A(i,j+1,1,1)-A(i,j,1,1))/dx(2)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1)
    Info%aux(i,j,1,2)=-(A(i+1,j,1,1)-A(i,j,1,1))/dx(1)
  END FORALL

ELSEIF (nDim == 3) THEN
  !Now the vector potential is a 3 component vector
  ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1:Info%mX(3)+1,3))

  !We could do this with three separate loops to avoid calculating 
  !A(:,:,Info%mX(3)+1,3),A(:,Info%mX(2)+1,:,2), and A(Info%mX(1)+1,:,:,1)
  !Since these do not get used in taking the curl...

  DO i=1, Info%mX(1)+1
    x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    DO j=1, Info%mX(2)+1
      y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
      DO k=1, Info%mX(3)+1
        z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
        A(i,j,k,1)=Ax(x+half*dx(1),y,z)
        A(i,j,k,2)=Ay(x,y+half*dx(2),z)
        A(i,j,k,3)=Az(x,y,z+half*dx(3))
      END DO
    END DO
  END DO

  !Then take the curl of A to get B
  FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2),k=1:Info%mX(3))
    Info%aux(i,j,k,1)=+(A(i,j+1,k,3)-A(i,j,k,3))/dx(2)-(A(i,j,k+1,2)-A(i,j,k,2))/dx(3)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1,k=1:Info%mX(3))
    Info%aux(i,j,k,2)=+(A(i,j,k+1,1)-A(i,j,k,1))/dx(3)-(A(i+1,j,k,3)-A(i,j,k,3))/dx(1)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2),k=1:Info%mX(3)+1)
    Info%aux(i,j,k,3)=+(A(i+1,j,k,2)-A(i,j,k,2))/dx(1)-(A(i,j+1,k,1)-A(i,j,k,1))/dx(2)
  END FORALL
END IF