Units in Astrobear

Astrobear uses something like rationalized electromagnetic units (extra factor of \sqrt{4 \pi} in the electric and magnetic fields) - or Lorentz-Heaviside but scaling the E field by c and the charge density \rho (and current J ) by \frac{1}{c}

E = c E^{LH} = \frac{c }{\sqrt{4\pi}}E^{G}

\rho = \frac{1}{c}\rho^{LH} = \frac{\sqrt{4 \pi} }{c}\rho^{G}

J = \frac{1}{c} J^{LH} = \frac{\sqrt{4 \pi}}{c} J^{G}

B = B^{LH} = \frac{1}{\sqrt{4\pi}} B^{G}

Using the approach in the appendix of Jackson, we have

k_1 = \frac{c^2}{4 \pi}

k_2 = \frac{1}{4\pi}

k_3 = 1

\alpha = 1

\mu_0 = 1

\epsilon_0 = \frac{1}{c^2}

This allows us to write Maxwell's equations as

\nabla \cdot \mathbf{E} = c^2 \rho

\nabla \times \mathbf{B} = \mathbf{J} + \frac{1}{c^2} \frac{\partial \mathbf{E}}{\partial t}

\nabla \times \mathbf{E} = - \frac{\partial \mathbf{B}}{\partial t}

\nabla \cdot \mathbf{B} = 0

as well as

Lorentz Force Law	\mathbf{F} = q \left ( \mathbf{E} + \mathbf{v} \times \mathbf{B} \right)
Coulomb's Law	\mathbf{F} = -\frac{c^2}{4 \pi} \frac{q_1 q_2}{r^2}\hat{\mathbf{r}}

Most of the time we don't care about E, \rho, or J, however for the Hall MHD terms, we need to determine the electron charge in our system of units.

For our system of equations, the elementary charge is

e = 4.80320425 \times 10^{-10} \mbox{statC} = 1.70269157 \times 10^{-9} \mbox{hsu} = 5.67956774 \times 10^{-20} \left [\mbox{g}^{1/2} \mbox{cm}^{1/2} \right]

MHD equations

Using the system of units described above, we can write the following simplified equations

Induction equation	\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times \mathbf{E}
Ampere's equation (without Maxwell's correction)	\mathbf{J} = \nabla \times \mathbf{B}
Lorentz Force Law	\frac{\partial \rho \mathbf{v}}{\partial t}=\mathbf{J} \times \mathbf{B}
Poynting's theorem (electromagnetic energy)	\frac{\partial e}{\partial t} = - \nabla \cdot \mathbf{S} - \mathbf{J}\cdot \mathbf{E}
Poynting vector	\mathbf{S} = \mathbf{E} \times \mathbf{B}

Equations for momentum, magnetic energy, and magnetic field

Subsituting Ampere's equation and the definition of the Poynting vector into the equations for momentum, magnetic energy, and the magnetic field, we arrive at

\frac{\partial \rho \mathbf{v}}{\partial t} = \mathbf{J} \times \mathbf{B} = \left (\nabla \times \mathbf{B} \right) \times \mathbf{B} = \nabla \cdot \left [ \mathbf{BB}-\frac{B^2}{2} \mathbf{I} \right ]

\frac{\partial e}{\partial t} = -\nabla \cdot \left ( \mathbf{E} \times \mathbf{B} \right) - \left ( \nabla \times \mathbf{B} \right) \cdot \mathbf{E}

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right)

Ohm's Law

Not the momentum equation only involves the magnetic field \mathbf{B} and is the same regardless of the electric field \mathbf{E}. Note it is also already in a form involving a total divergence making it straightforward to implement in a conservative finite volume scheme. The energy and induction equation, however, will depend on our choice of \mathbf{E}.

The \mathbf{E} field can be calculated using a Generalized Ohm's law

\frac{m_e}{n e^2} \frac{\partial J}{\partial t} - \frac{1}{ne}\nabla \cdot \mathbf{P_e} = \mathbf{E} + \mathbf{v} \times \mathbf{B} - \eta\mathbf{J} - \frac{1}{ne} \mathbf{J} \times \mathbf{B}

Under under certain conditions, we can ignore various terms resulting in the following approximations

Resistive + Hall MHD	\mathbf{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{J} + \frac{1}{ne} \mathbf{J} \times \mathbf{B}
Resistive MHD	\mathbf{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{J}
Ideal MHD	\mathbf{E} = - \mathbf{v}\times \mathbf{B}

Ideal MHD

For ideal MHD we substitute

\mathbf{E} = -\mathbf{v} \times \mathbf{B}

into the equations above and after some vector calculus (including a very messy vector quadruple product with a del operator), manage to write the energy term as a total divergence (and the field as a curl)

\frac{\partial e}{\partial t} = \nabla \cdot \left ( \left (\mathbf{v} \times \mathbf{B} \right ) \times \mathbf{B} \right) + \left ( \nabla \times \mathbf{B} \right) \cdot \left ( \mathbf{v} \times \mathbf{B} \right) = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \right]

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right)

Resistive MHD

For resistive MHD, it will be convenient to define the resistive component of the electric field \color{green}{\mathbf{E}_{\scriptscriptstyle R} = \eta \mathbf{J} = \eta \nabla \times \mathbf{B}}

Then our total electric field is given by

\mathbf{E} = -\mathbf{v} \times \mathbf{B} \color{green}{+ \mathbf{E}_{\scriptscriptstyle R}}

We then need to add this contribution to the induction equation \color{green}{-\nabla \times \mathbf{E}_{\scriptscriptstyle R}} and the energy equation where it enters both the Poynting flux -\color{green}{\nabla \cdot \left ( \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B} \right)} and the work done on the the charge distribution \color{green}{ - \mathbf{J} \cdot \mathbf{E}_{\scriptscriptstyle R}}.

However the work done on the charge distribution would need to also be added to the thermal energy cancelling any changes in the total energy due to the \color{green}{\mathbf{J} \cdot \mathbf{E}_{\scriptscriptstyle R}} term. This leaves us with.

\frac{\partial e}{\partial t} = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \color{green}{- \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B}} \right]

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{ - \nabla \times \mathbf{E}_{\scriptscriptstyle R}} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right )}

Note the resistive part of the induction equation

\frac{\partial \mathbf{B}}{\partial t} = \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right )}

can be expanded as

\frac{\partial \mathbf{B}}{\partial t} = \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right ) = \eta \nabla^2 \mathbf{B} - \vec{\nabla} \mathbf{B} \cdot \nabla \eta}

Which hints to the diffusive nature of the resistive term. It also involves a second derivative so the time stepping criteria will be \Delta t \propto \Delta x^2 and an implicit solve might be worth considering.

Hall+Resistive MHD

As in resistive MHD, it will be convenient to define the Hall component of the electric field \color{red}{\mathbf{E}_{\scriptscriptstyle H} =} \color{red}{ \frac{1}{ne} \mathbf{J} \times \mathbf{B} = \frac{1}{ne} \left ( \nabla \times \mathbf{B} \right ) \times \mathbf{B}}

Then our total electric field is given by

\mathbf{E} = -\mathbf{v} \times \mathbf{B} \color{green}{+ \mathbf{E}_{\scriptscriptstyle R}} \color{red}{+\mathbf{E}_{\scriptscriptstyle H}}

Again we need to add the additional electric field component to the energy equation (in both the Poynting Flux and the work on the charge distribution) as well as the induction equation. Also note in this case the electric field is perpendicular to the current - so it does no work on the charge distribution.

\mathbf{J} \cdot \color{red}{\mathbf{E}_{\scriptscriptstyle H}} = \mathbf{J} \cdot \color{red}{\frac{1}{ne} \mathbf{J} \times \mathbf{B}} = 0

So as in the Resistive case, only the Poynting Flux term needs to be added to the energy equation.

\frac{\partial \rho \mathbf{v}}{\partial t} = \left (\nabla \times \mathbf{B} \right) \times \mathbf{B} = \nabla \cdot \left [ \mathbf{BB}-\frac{B^2}{2} \mathbf{I} \right ]

\frac{\partial e}{\partial t} = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \color{green}{ - \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B}} \color{red}{ -\mathbf{E}_{\scriptscriptstyle H} \times \mathbf{B}} \right]

\frac{\partial B}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{- \nabla \times \mathbf{E}_{\scriptscriptstyle R}} \color{red}{- \nabla \times \mathbf{E}_{\scriptscriptstyle H}}

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