Summary of the Compressible MHD Code Implementation

A possible implementation of the algorithm is:

• Step 1: Read in initial conditions ${\bf U}({\bf x},0)$, evaluate the fields at the element quadrature points. Set n = 0.
• Step 2: Calculate the upwinded fluxes ${\bf \hat{F}_{\bf n}}$ at the Gauss quadrature points on the element interfaces. At domain boundaries use the prescribed boundary conditions for the exterior values of the fields.
• Step 3: Calculate the Euler flux terms: ${\bf F}^{Euler}_x$, ${\bf F}^{Euler}_y$ and ${\bf F}^{Euler}_z$ at the element quadrature points.
• Step 4: For each component of the state vector $U_k = ({\bf U}({\bf x},t^n))_k$ calculate $(\frac{\partial {\bf F}^{Euler}_x}{\partial x}+\frac{\partial {\bf F}^{Euler}_y}{\partial y}+\frac{\partial {\bf F}^{Euler}_z}{\partial z})_k$.
• Step 5: Interpolate the upwinded flux to the correct quadrature for the element interface. Scale the fluxes with the edge Jacobian divided by the volume Jacobians. Add this to the divergence of the Euler fluxes calculated in Step 4.
• Step 6: Take inner product of the fluxes with the orthogonal basis. Evaluate the resulting polynomials at the quadrature points.
• Step 7: Calculate the spatial derivatives of the primitive fields. For example, we will need $\frac{\partial u}{\partial x}$. We evaluate this as described in section 7.1.1 using the DGM to compensate for the derivative across element interfaces, and the average of u at the element interface.
• Step 8: Use the derivatives of the primitive fields to construct the Viscous flux terms: ${\bf F}^{Visc}_x$, ${\bf F}^{Visc}_y$ and ${\bf F}^{Visc}_z$.
• Step 9: Take the divergence of the viscous flux terms.
• Step 10: Subtract the divergence of the Euler fluxes obtained in Step 6 from the divergence of viscous fluxes obtained in Step 9.
• Step 11: Take inner product of the result from Step 10 with the orthogonal basis. Evaluate the resulting polynomials at the quadrature points and place it in ${\bf Uf}({\bf x},t^{n-q})$
• Step 12: Update the state vector ${\bf U}({\bf x},t^{n+1}) = {\bf U}({\bf x},t^{n}) + \Delta t \sum_q \beta_q {\bf Uf}({\bf x},t^{n-q})$ using an Adams-Bashforth integration scheme.
• Step 13: Increase n by one. If tⁿ is less than the termination time return to Step 2.
• Step 14: Output final values of the state vector ${\bf U}({\bf x},t^{end})$.