Using USim to solve the Magnetohydrodynamic Equations

In Using USim to solve the Euler Equations we discussed the basic methods used by USim to solve the Euler equations. In this tutorial, we show how USim can be used to integrate the magnetohydrodynamic (MHD) equations for problems in one-dimension. This tutorial is based on the Shock Tube (shockTube.pre) example. The Shock Tube simulation is designed to set up a variety of tube simulations including those by Einfeldt, Sod, Liska & Wendroff, Brio & Wu, and Ryu & Jones. In this tutorial, we will look at the Brio & Wu shock Tube, described by:

Brio, M., & Wu, C.~C. (1988), Journal of Computational Physics, 75, 400

Note

While we will be following the Shock Tube (shockTube.pre) example, we will not reproduce the example file in its entirety. The shockTube example is designed to set up a variety of simulations, and by using “Flags” (if statements) in the shockTube.pre file, only certain parts of the file are used when simulating the Brio & Wu Shock Tube. Therefore, at the end of this tutorial, our input will not be an exact copy of the shockTube.pre file but will directly represent only the Brio & Wu Shock Tube portion of it.

In this example we will look at the use of equations for ideal compressible magnetohydrodynamics (the MHD equations) in one-dimension. It is appropriate to use these equations for transonic, supersonic and hypersonic flows (Mach numbers of 0.1 and above) where compressibility effects are important, at high Reynolds numbers where the effects of viscosity and thermal conductivity are relatively unimportant and at high magnetic Reynolds number where Ohmic resistivity is relatively unimportant.

We will follow the pattern established in Using USim to solve the Euler Equations in order that the user can see the differences between solving the MHD equations and the Euler equations. These are summarized in Notes at the end of each section of the tutorial.

The Magnetohydrodynamic Equations

The MHD equations for an adiabatic gas can be written as a hyperbolic conservation law, which has the form:

\[\notag \begin{align} \frac{\partial \mathbf{q}}{\partial t} + \nabla\cdot\left[ \mathcal{F} \left( \mathbf{w} \right) \right] = 0 \end{align}\]

where \(\mathbf{q}\) is a vector of conserved variables (e.g. density, momentum, total energy), \(\mathcal{F}\left( \mathbf{w} \right)\) is a non-linear flux tensor computed from a vector of primitive variables, (e.g. density, velocity, pressure), and \(\mathbf{w} = \mathbf{w}(\mathbf{q})\). Assuming an adiabatic equations of state, the MHD equations can be written in this form such that:

\[\notag \begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot\left[ \rho\,\mathbf{u} \right] = 0 \\ \frac{\partial \rho\,\mathbf{u}}{\partial t} + \nabla\cdot\left[ \rho\,\mathbf{u}\,\mathbf{u}^{T} - \mathbf{b}\,\mathbf{b}^{T} + \mathbb{I} \left(P+ \tfrac{1}{2}\left|\mathbf{b}\right|^{2}\right) \right] = 0 \\ \frac{\partial E}{\partial t} + \nabla\cdot\left[ \left(E + P \right) \mathbf{u} + \mathbf{e} \times \mathbf{b} \right] = 0 \\ \frac{\partial \mathbf{b^{\mathrm{plasma}}}}{\partial t} + \nabla\times\mathbf{e} + \nabla \psi= 0 \\ \frac{\partial \psi}{\partial t} + \nabla\cdot\left[ c^{2}_{\mathrm{fast}} \mathbf{b} \right] = 0 \end{align}\]

Here, \(\mathbb{I}\) is the identity matrix, \(P = \rho\epsilon(\gamma-1)\) is the pressure of an ideal gas, \(\epsilon\) is the specific internal energy and \(\gamma\) is the adiabatic index (ratio of specific heats). The quantity \(c_{\mathrm{fast}}\) corresponds to the fastest wave speed over the entire simulation domain; divergence errors are advected out of the domain with this speed.

The MHD equations are represented in USim by mhdDednerEqn. USim solves these equations by calculating an upwind approximation of the flux tensor, \(\mathcal{F}(\mathbf{w})\) (see classicMusclUpdater (1d, 2d, 3d)) and then using this approximation to advance the conserved state from time \(t\) to \(t+\Delta t\) (see multiUpdater (1d, 2d, 3d)).

Solving the MHD Equations in One Dimension

The first step in any USim simulation is to import the macros that are needed to define the simulation. For the ShockTube example, there are two macros that are needed:

$ import fluidsBase.mac
$ import idealmhd.mac

These macros provide basic capabilities for setting up a USim simulation of the ideal MHD Equations. We can now define global parameters that tell USim basic information about the simulation. This is done through the use of the initializeFluidSimulation macro:

initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

The parameters for the version of the macro used for the mhdDednerEqn are documented at Ideal MHD Macro. For completeness, we include them here, showing how the parameters specified here map onto the parameters used in the initializeFluidSimulation macro:

NDIM (dimensionality):
1,2,3. Number of dimensions for the simulation
0.0 (tStart):
Start time for simulation
END_TIME (tEnd):
End time for simulation
NUMDUMPS (numFrames):
Number of data outputs
CFL (cflNum):
Cfl limit, typically \(\Delta t = cflNum * \Delta x / V_{max}\)
GAS_GAMMA (gammaIn):
Adiabatic index for ideal gas eqn. of state. Pressure = (gammaIn - 1.0) * density * internal energy
MU0 (muIn):
Permeability of free space
WRITE_RESTART (writeRestartIn):
Output data required for simulation restart
DEBUG (debugIn):
Run simulation in debug mode

USim includes the ability to perform in place substitution of variables within the input file, as described in Symbol Definition. This means that we can define the options given above earlier in the input file and USim will replace them in the initializeFluidSimulation macro:

# X-extent of domain
$PAR_LENGTH = 1.0
# Y-extent of domain
$PERP_LENGTH = 1.0
# Zones parallel to shock normal direction
$PAR_ZONES = 256
# Zones perpendicular to shock normal direction
$PERP_ZONES = 256
# Shock tube to solve
$SHOCK_TUBE = "BRIOWU"
# adiabatic index
$GAS_GAMMA = 5.0/3.0
# permeability of free-space
$MU0 = 1.0
# Atmospheric Pressure
$REFERENCE_PRESSURE = 1.0 # [Pa]
# density of air
$REFERENCE_DENSITY = 1.0 # [kg/m^3]
# end time for simulation (units of # of sound wave crossings)
$TEND = 0.125
# number of frames
$NUMDUMPS = 10
# Whether to use a diffusive (but robust) integration scheme
$DIFFUSIVE = False
# Order in time
$TIME_ORDER = "second"
# Write data for restarting the simulation
$WRITE_RESTART = False
# Output info for debugging purposes
$DEBUG = False
# Default dimensionality
$NDIM = 1

Note that we have included a range of other variables that are defined in the Shock Tube (shockTube.pre) example. We will refer to these later in this tutorial as needed. In the above, the variable TEND is given in units of the the number of times a sound wave crosses the grid. In order to use this to specify the end time for the simulation (END_TIME), we convert this to have units of time through:

# nominal speed of sound
$ c0 = math.sqrt(GAS_GAMMA*REFERENCE_PRESSURE/REFERENCE_DENSITY) # [m/s]

# Set end time according to (user specified)
# number of times a wave crosses the box.
$END_TIME = $PERP_LENGTH*TEND/c0$

The specification of the CFL parameter is described in the next section. Our simulation starts off as:

# Import macros to setup simulation
$ import fluidsBase.mac
$ import idealmhd.mac

# X-extent of domain
$PAR_LENGTH = 1.0
# Y-extent of domain
$PERP_LENGTH = 1.0
# Zones parallel to shock normal direction
$PAR_ZONES = 256
# Zones perpendicular to shock normal direction
$PERP_ZONES = 256
# Shock tube to solve
$SHOCK_TUBE = "BRIOWU"
# adiabatic index
$GAS_GAMMA = 5.0/3.0
# permeability of free-space
$MU0 = 1.0
# Atmospheric Pressure
$REFERENCE_PRESSURE = 1.0 # [Pa]
# density of air
$REFERENCE_DENSITY = 1.0 # [kg/m^3]
# end time for simulation (units of # of sound wave crossings)
$TEND = 0.125
# number of frames
$NUMDUMPS = 10
# Whether to use a diffusive (but robust) integration scheme
$DIFFUSIVE = False
# Order in time
$TIME_ORDER = "second"
# Write data for restarting the simulation
$WRITE_RESTART = False
# Output info for debugging purposes
$DEBUG = False
# Default dimensionality
$NDIM = 1

# nominal speed of sound
$ c0 = math.sqrt(GAS_GAMMA*REFERENCE_PRESSURE/REFERENCE_DENSITY) # [m/s]

# Set end time according to (user specified)
# number of times a wave crosses the box.
$END_TIME = $PERP_LENGTH*TEND/c0$

# Initialize a USim simulation
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

Note

Compared to the Euler equations discussed in Using USim to solve the Euler Equations, there are three differences:

  1. We have replaced $ import euler.mac with $ import idealmhd.mac to change the system of equations from the Euler equations to the ideal MHD equations.

  2. We have defined the permeability of free-space, \(\mu_0\) through the parameter MU0 (here MU0 = 1.0).

  3. We have added \(\mu_0\) to the parameters used to initialized the simulation, replacing:

    initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,WRITE_RESTART,DEBUG)

    with:

    initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

Adding a Simulation Grid

The next step in setting up a USim simulation is to specify a simulation grid. The Shock Tube (shockTube.pre) uses a cart (1d, 2d, 3d) grid, which is added to the simulation through the addGrid macro:

addGrid(lowerBounds, upperBounds, numCells, periodicDirections)

This options for this macro are documented in Grid Macro. For completeness, we include them here:

lowerBounds:
Vector of coordinates for lower edge of grid, lowerBounds = [ XMIN YMIN ZMIN ]
upperBounds:
Vector of coordinates for upper edge of grid, upperBounds = [ XMAX YMAX ZMAX ]
numCells:
Vector of number of cells in grid, numCells = [ NX NY NZ ]
periodicDirections:

List of directions that are periodic

periodicDirections = [ 0 ] (x-direction periodic)
periodicDirections = [ 0 1 ] (x,y-directions periodic)
periodicDirections = [ 0 1 2 ] (x,y,z-directions periodic)

The options are setup in the ShockTube example according to the dimensionality of the simulation, as specified by the variable NDIM above:

$XMIN = -0.5*PERP_LENGTH
$XMAX =  0.5*PERP_LENGTH
$YMIN = -0.5*PAR_LENGTH
$YMAX =  0.5*PAR_LENGTH
$ZMIN = -0.5*PAR_LENGTH
$ZMAX =  0.5*PAR_LENGTH

$if NDIM==1
$ CFL = 0.5
$ numCells = [PERP_ZONES]
$ periodicDirections = []
$ lowerBounds = [XMIN]
$ upperBounds = [XMAX]
$ else
$if NDIM==2
$ CFL = 0.4
$ numCells = [PERP_ZONES, PAR_ZONES]
# Make the direction perpendicular to the shock
# normal periodic.
$ periodicDirections = [1]
$ lowerBounds = [XMIN,  YMIN]
$ upperBounds = [XMAX,  YMAX]
$ else
$ CFL = 0.32
$ numCells = [PERP_ZONES, PAR_ZONES, PAR_ZONES]
# Make the directions perpendicular to the shock
# normal periodic.
$ periodicDirections = [1,2]
$ lowerBounds = [XMIN,  YMIN,  ZMIN]
$ upperBounds = [XMAX,  YMAX,  ZMAX]
$ endif
$ endif

Note as well that the CFL condition that is used for the simulation changes according to the dimension of the simulation. This is because the scheme for solving the hyperbolic equations has different stability requirements in different dimensions:

\[\notag \begin{align} \mathrm{CFL} <= \frac{1}{\mathrm{NDIM}} \end{align}\]

So, our simulation now looks like:

# Import macros to setup simulation
$ import fluidsBase.mac
$ import idealmhd.mac

# X-extent of domain
$PAR_LENGTH = 1.0
# Y-extent of domain
$PERP_LENGTH = 1.0
# Zones parallel to shock normal direction
$PAR_ZONES = 256
# Zones perpendicular to shock normal direction
$PERP_ZONES = 256
# Shock tube to solve
$SHOCK_TUBE = "BRIOWU"
# adiabatic index
$GAS_GAMMA = 5.0/3.0
# permeability of free-space
$MU0 = 1.0
# Atmospheric Pressure
$REFERENCE_PRESSURE = 1.0 # [Pa]
# density of air
$REFERENCE_DENSITY = 1.0 # [kg/m^3]
# end time for simulation (units of # of sound wave crossings)
$TEND = 0.125
# number of frames
$NUMDUMPS = 10
# Whether to use a diffusive (but robust) integration scheme
$DIFFUSIVE = False
# Order in time
$TIME_ORDER = "second"
# Write data for restarting the simulation
$WRITE_RESTART = False
# Output info for debugging purposes
$DEBUG = False
# Default dimensionality
$NDIM = 1

# nominal speed of sound
$ c0 = math.sqrt(GAS_GAMMA*REFERENCE_PRESSURE/REFERENCE_DENSITY) # [m/s]

# Set end time according to (user specified)
# number of times a wave crosses the box.
$END_TIME = $PERP_LENGTH*TEND/c0$

$XMIN = -0.5*PERP_LENGTH
$XMAX =  0.5*PERP_LENGTH
$YMIN = -0.5*PAR_LENGTH
$YMAX =  0.5*PAR_LENGTH
$ZMIN = -0.5*PAR_LENGTH
$ZMAX =  0.5*PAR_LENGTH

$if NDIM==1
$ CFL = 0.5
$ numCells = [PERP_ZONES]
$ periodicDirections = []
$ lowerBounds = [XMIN]
$ upperBounds = [XMAX]
$ else
$if NDIM==2
$ CFL = 0.4
$ numCells = [PERP_ZONES, PAR_ZONES]
# Make the direction perpendicular to the shock
# normal periodic.
$ periodicDirections = [1]
$ lowerBounds = [XMIN,  YMIN]
$ upperBounds = [XMAX,  YMAX]
$ else
$ CFL = 0.32
$ numCells = [PERP_ZONES, PAR_ZONES, PAR_ZONES]
# Make the directions perpendicular to the shock
# normal periodic.
$ periodicDirections = [1,2]
$ lowerBounds = [XMIN,  YMIN,  ZMIN]
$ upperBounds = [XMAX,  YMAX,  ZMAX]
$ endif
$ endif

# Initialize a USim simulation
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

# Setup the grid
addGrid(lowerBounds, upperBounds, numCells, periodicDirections)

Note

The process of adding the simulation grid is identical for both the Euler equations and the MHD equations. This is true for all USim simulations.

Creating a Fluid Simulation

The next step in setting up a USim simulation is to create the basic set of variables need to simulate the fluid. This is accomplished through the createFluidSimulation macro, documented at Ideal MHD Macro:

createFluidSimulation()

Now that these variables have been created, it is possible to specify the distribution of fluid on the grid. This is a three-step process:

  1. Use addVariable to add variables that are independent of grid position.
  2. Use addPreExpression to add quantities that are functions of grid position, variables and any previously defined PreExpression in this block. Evaluated before expressions and the result is not accessible outside of this block. Any number of PreExpressions can be added.
  3. Use addExpression to define each initial condition for the fluid. There is one expression for density, each component of momentum and the total energy. The order of the exprssions correspond to the order in the state vector and there can only be one expression per entry in the state vector.

For the Brio & Wu variant of the ShockTube example, the variables added in Step 1 in the above process are:

addVariable(pi_value,math.pi)
addVariable(gas_gamma,GAS_GAMMA)
addVariable(densityL,DENSITY_L)
addVariable(densityR,DENSITY_R)
addVariable(pressureL,PRESSURE_L)
addVariable(pressureR,PRESSURE_R)
addVariable(normalVelocityL,NORMAL_VELOCITY_L)
addVariable(normalVelocityR,NORMAL_VELOCITY_R)
addVariable(perpendicularVelocityL,PERPENDICULAR_VELOCITY_L)
addVariable(perpendicularVelocityR,PERPENDICULAR_VELOCITY_R)
addVariable(tangentialVelocityL,TANGENTIAL_VELOCITY_L)
addVariable(tangentialVelocityR,TANGENTIAL_VELOCITY_R)
addVariable(mu0,MU0)
addVariable(normalFieldL,NORMAL_FIELD_L)
addVariable(normalFieldR,NORMAL_FIELD_R)
addVariable(perpendicularFieldL,PERPENDICULAR_FIELD_L)
addVariable(perpendicularFieldR,PERPENDICULAR_FIELD_R)
addVariable(tangentialFieldL,TANGENTIAL_FIELD_L)
addVariable(tangentialFieldR,TANGENTIAL_FIELD_R)

The variables are then used in the PreExpressions that are defined in Step 2:

addPreExpression(rho = if (x>0.0, densityR,               densityL))
addPreExpression(pr  = if (x>0.0, pressureR,              pressureL))
addPreExpression(vx  = if (x>0.0, normalVelocityR,        normalVelocityL))
addPreExpression(vy  = if (x>0.0, perpendicularVelocityR, perpendicularVelocityL))
addPreExpression(vz  = if (x>0.0, tangentialVelocityR,    tangentialVelocityL))
addPreExpression(bx  = if (x>0.0, normalFieldR,           normalFieldL))
addPreExpression(by  = if (x>0.0, perpendicularFieldR,    perpendicularFieldL))
addPreExpression(bz  = if (x>0.0, tangentialFieldR,       tangentialFieldL))
addPreExpression(psi = 0.0)

Finally, these PreExpressions are used to specify the initial conditions defined in Step 3:

addExpression(rho)
addExpression(rho*vx)
addExpression(rho*vy)
addExpression(rho*vz)
addExpression(pr/(gas_gamma-1.0) + 0.5*rho*(vx*vx+vy*vy+vz*vz)+0.5*((bx*bx+by*by+bz*bz)/mu0))
addExpression(bx)
addExpression(by)
addExpression(bz)
addExpression(psi)

The ShockTube example includes many different initial conditions that are based on a range of cases proposed in the literature. The Brio & Wu shock tube is specified through $SHOCK_TUBE = “BRIOWU” above. This corresponds to the following variable definitions:

$ DENSITY_L = $1.0*REFERENCE_DENSITY$
$ DENSITY_R = $0.125*REFERENCE_DENSITY$
$ PRESSURE_L = $1.0*REFERENCE_PRESSURE$
$ PRESSURE_R =  $0.1*REFERENCE_PRESSURE$
$ NORMAL_VELOCITY_L = 0.0
$ NORMAL_VELOCITY_R = 0.0
$ PERPENDICULAR_VELOCITY_L = 0.0
$ PERPENDICULAR_VELOCITY_R = 0.0
$ TANGENTIAL_VELOCITY_L = 0.0
$ TANGENTIAL_VELOCITY_R = 0.0
$ NORMAL_FIELD_L = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ NORMAL_FIELD_R = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_L = $1.0*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_R = $-1.0*math.sqrt(REFERENCE_DENSITY)$
$ TANGENTIAL_FIELD_L = 0.0
$ TANGENTIAL_FIELD_R = 0.0

Note that this initialization procedure sets up up the following initial condition:

\[\begin{eqnarray} & \rho &= 1.0 \rho_0 \quad x\le0.0\\ & \rho &= 0.125 \rho_0 \quad x>0.0\\ & P_g &= 1.0 P_0 \quad x\le0.0 \\ & P_g &= 0.1 P_0 \quad x>0.0 \\ & b_y &= \sqrt{\rho_0} \quad x\le0.0 \\ & b_y &=-\sqrt{\rho_0} \quad x>0.0 \\ & b_x &= 0.75 \sqrt{\rho_0} \end{eqnarray}\]

where \(\rho_0,P_0\) are the reference density and pressure respectively. Our simulation now looks like:

# Import macros to setup simulation
$ import fluidsBase.mac
$ import idealmhd.mac

# X-extent of domain
$PAR_LENGTH = 1.0
# Y-extent of domain
$PERP_LENGTH = 1.0
# Zones parallel to shock normal direction
$PAR_ZONES = 256
# Zones perpendicular to shock normal direction
$PERP_ZONES = 256
# Shock tube to solve
$SHOCK_TUBE = "BRIOWU"
# adiabatic index
$GAS_GAMMA = 5.0/3.0
# permeability of free-space
$MU0 = 1.0
# Atmospheric Pressure
$REFERENCE_PRESSURE = 1.0 # [Pa]
# density of air
$REFERENCE_DENSITY = 1.0 # [kg/m^3]
# end time for simulation (units of # of sound wave crossings)
$TEND = 0.125
# number of frames
$NUMDUMPS = 10
# Whether to use a diffusive (but robust) integration scheme
$DIFFUSIVE = False
# Order in time
$TIME_ORDER = "second"
# Write data for restarting the simulation
$WRITE_RESTART = False
# Output info for debugging purposes
$DEBUG = False
# Default dimensionality
$NDIM = 1

# nominal speed of sound
$ c0 = math.sqrt(GAS_GAMMA*REFERENCE_PRESSURE/REFERENCE_DENSITY) # [m/s]

# Set end time according to (user specified)
# number of times a wave crosses the box.
$END_TIME = $PERP_LENGTH*TEND/c0$

$XMIN = -0.5*PERP_LENGTH
$XMAX =  0.5*PERP_LENGTH
$YMIN = -0.5*PAR_LENGTH
$YMAX =  0.5*PAR_LENGTH
$ZMIN = -0.5*PAR_LENGTH
$ZMAX =  0.5*PAR_LENGTH

$if NDIM==1
$ CFL = 0.5
$ numCells = [PERP_ZONES]
$ periodicDirections = []
$ lowerBounds = [XMIN]
$ upperBounds = [XMAX]
$ else
$if NDIM==2
$ CFL = 0.4
$ numCells = [PERP_ZONES, PAR_ZONES]
# Make the direction perpendicular to the shock
# normal periodic.
$ periodicDirections = [1]
$ lowerBounds = [XMIN,  YMIN]
$ upperBounds = [XMAX,  YMAX]
$ else
$ CFL = 0.32
$ numCells = [PERP_ZONES, PAR_ZONES, PAR_ZONES]
# Make the directions perpendicular to the shock
# normal periodic.
$ periodicDirections = [1,2]
$ lowerBounds = [XMIN,  YMIN,  ZMIN]
$ upperBounds = [XMAX,  YMAX,  ZMAX]
$ endif
$ endif

# Parameters to specify the fluid state at t=0.0
$ DENSITY_L = $1.0*REFERENCE_DENSITY$
$ DENSITY_R = $0.125*REFERENCE_DENSITY$
$ PRESSURE_L = $1.0*REFERENCE_PRESSURE$
$ PRESSURE_R =  $0.1*REFERENCE_PRESSURE$
$ NORMAL_VELOCITY_L = 0.0
$ NORMAL_VELOCITY_R = 0.0
$ PERPENDICULAR_VELOCITY_L = 0.0
$ PERPENDICULAR_VELOCITY_R = 0.0
$ TANGENTIAL_VELOCITY_L = 0.0
$ TANGENTIAL_VELOCITY_R = 0.0
$ NORMAL_FIELD_L = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ NORMAL_FIELD_R = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_L = $1.0*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_R = $-1.0*math.sqrt(REFERENCE_DENSITY)$
$ TANGENTIAL_FIELD_L = 0.0
$ TANGENTIAL_FIELD_R = 0.0

# Initialize a USim simulation
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

# Setup the grid
addGrid(lowerBounds, upperBounds, numCells, periodicDirections)

# Create data structures needed for the simulation
createFluidSimulation()

# Step 1: Add Variables
addVariable(pi_value,math.pi)
addVariable(gas_gamma,GAS_GAMMA)
addVariable(densityL,DENSITY_L)
addVariable(densityR,DENSITY_R)
addVariable(pressureL,PRESSURE_L)
addVariable(pressureR,PRESSURE_R)
addVariable(normalVelocityL,NORMAL_VELOCITY_L)
addVariable(normalVelocityR,NORMAL_VELOCITY_R)
addVariable(perpendicularVelocityL,PERPENDICULAR_VELOCITY_L)
addVariable(perpendicularVelocityR,PERPENDICULAR_VELOCITY_R)
addVariable(tangentialVelocityL,TANGENTIAL_VELOCITY_L)
addVariable(tangentialVelocityR,TANGENTIAL_VELOCITY_R)
addVariable(mu0,MU0)
addVariable(normalFieldL,NORMAL_FIELD_L)
addVariable(normalFieldR,NORMAL_FIELD_R)
addVariable(perpendicularFieldL,PERPENDICULAR_FIELD_L)
addVariable(perpendicularFieldR,PERPENDICULAR_FIELD_R)
addVariable(tangentialFieldL,TANGENTIAL_FIELD_L)
addVariable(tangentialFieldR,TANGENTIAL_FIELD_R)

# Step 2: Add Pre-Expressions
addPreExpression(rho = if (x>0.0, densityR,               densityL))
addPreExpression(pr  = if (x>0.0, pressureR,              pressureL))
addPreExpression(vx  = if (x>0.0, normalVelocityR,        normalVelocityL))
addPreExpression(vy  = if (x>0.0, perpendicularVelocityR, perpendicularVelocityL))
addPreExpression(vz  = if (x>0.0, tangentialVelocityR,    tangentialVelocityL))
addPreExpression(bx  = if (x>0.0, normalFieldR,           normalFieldL))
addPreExpression(by  = if (x>0.0, perpendicularFieldR,    perpendicularFieldL))
addPreExpression(bz  = if (x>0.0, tangentialFieldR,       tangentialFieldL))
addPreExpression(psi = 0.0)

# Step 3: Add expressions specifying initial condition on density,
# momentum, total energy, magnetic field and correction potential
addExpression(rho)
addExpression(rho*vx)
addExpression(rho*vy)
addExpression(rho*vz)
addExpression(pr/(gas_gamma-1.0) + 0.5*rho*(vx*vx+vy*vy+vz*vz)+0.5*((bx*bx+by*by+bz*bz)/mu0))
addExpression(bx)
addExpression(by)
addExpression(bz)
addExpression(psi)

Note

Compared to the Euler equations discussed in Using USim to solve the Euler Equations, there are three important differences in specifying the initial conditions:

  1. We have to specify the three components of the magnetic field, here these are denoted \(b_{\{x,y,z\}}\).
  2. We have to specify the correction potential, here denoted, psi. Typically, \(\psi = 0\) at \(t=0\).
  3. We have to change the definition of the total energy from \(E = \frac{P}{\Gamma-1} + \frac{\rho |\mathbf{u}|^2}{2}\) to \(E = \frac{P}{\Gamma-1} + \frac{\rho |\mathbf{u}|^2}{2} + \frac{|\mathbf{b}|^2}{2}\); that is we include the contribution of the magnetic field to the total energy

Evolving the Fluid

USim implements the well-known MUSCL scheme to advance the conserved variables in time. There is a simple macro that can be called to implement this scheme as shown below:

finiteVolumeScheme(DIFFUSIVE)

This macro is documented at Ideal MHD Macro. It’s purpose is to compute the numerical flux for the hyperbolic system:

\[\notag \begin{align} \nabla\cdot\left[ \mathcal{F} \left( \mathbf{w} \right) \right] - \mathcal{S} \left( \mathbf{w} \right) \end{align}\]

The next part of evolving the fluid is to apply boundary conditions at the left and right of the domain to ensure that at the next time step, physically-valid data is used to update the conserved state. Without this, the simulation will fail. It is possible to specify arbitrary boundary conditions in USim. For the Brio & Wu shock tube example considered here, appropriate boundary conditions are outflow (“open”) boundary conditions at both ends of the domain

boundaryCondition(copy,left)
boundaryCondition(copy,right)

The copy boundary condition block is described in copy (1d, 2d, 3d). This boundary condition updater copies the values on the layer next to the ghost cells into the ghost cells - this is equivalent to a zero derivative boundary condition. If we are evolving in more than one-dimension, we have to specify boundary conditions on the rest of the domain boundaries. For this example, this is done using periodic boundary conditions (periodicCartBc (1d, 2d, 3d)):

boundaryCondition(periodic)

The final element of advancing the conserved quantities from time \(t\) to \(t+\Delta t\) is to apply a time integration scheme, specified through:

timeAdvance(TIME_ORDER)

This applies an explicit Runge-Kutta time-integration scheme with the order of accuracy determined by the TIME_ORDER parameter. TIME_ORDER can be one of first, second, third or fourth according to the desired order of accuracy.

Note

The method for evolving the fluid is identical for the Euler equations and the MHD equations in one-dimension.

Putting it all Together

The final step in the USim simulation is to add:

runFluidSimulation()

This tells USim that we’re done specifying the simulation and that it can be run. So, our simulation now looks like:

# Import macros to setup simulation
$ import fluidsBase.mac
$ import idealmhd.mac

# X-extent of domain
$PAR_LENGTH = 1.0
# Y-extent of domain
$PERP_LENGTH = 1.0
# Zones parallel to shock normal direction
$PAR_ZONES = 256
# Zones perpendicular to shock normal direction
$PERP_ZONES = 256
# Shock tube to solve
$SHOCK_TUBE = "BRIOWU"
# adiabatic index
$GAS_GAMMA = 5.0/3.0
# permeability of free-space
$MU0 = 1.0
# Atmospheric Pressure
$REFERENCE_PRESSURE = 1.0 # [Pa]
# density of air
$REFERENCE_DENSITY = 1.0 # [kg/m^3]
# end time for simulation (units of # of sound wave crossings)
$TEND = 0.125
# number of frames
$NUMDUMPS = 10
# Whether to use a diffusive (but robust) integration scheme
$DIFFUSIVE = False
# Order in time
$TIME_ORDER = "second"
# Write data for restarting the simulation
$WRITE_RESTART = False
# Output info for debugging purposes
$DEBUG = False
# Default dimensionality
$NDIM = 1

# nominal speed of sound
$ c0 = math.sqrt(GAS_GAMMA*REFERENCE_PRESSURE/REFERENCE_DENSITY) # [m/s]

# Set end time according to (user specified)
# number of times a wave crosses the box.
$END_TIME = $PERP_LENGTH*TEND/c0$

$XMIN = -0.5*PERP_LENGTH
$XMAX =  0.5*PERP_LENGTH
$YMIN = -0.5*PAR_LENGTH
$YMAX =  0.5*PAR_LENGTH
$ZMIN = -0.5*PAR_LENGTH
$ZMAX =  0.5*PAR_LENGTH

$if NDIM==1
$ CFL = 0.5
$ numCells = [PERP_ZONES]
$ periodicDirections = []
$ lowerBounds = [XMIN]
$ upperBounds = [XMAX]
$ else
$if NDIM==2
$ CFL = 0.4
$ numCells = [PERP_ZONES, PAR_ZONES]
# Make the direction perpendicular to the shock
# normal periodic.
$ periodicDirections = [1]
$ lowerBounds = [XMIN,  YMIN]
$ upperBounds = [XMAX,  YMAX]
$ else
$ CFL = 0.32
$ numCells = [PERP_ZONES, PAR_ZONES, PAR_ZONES]
# Make the directions perpendicular to the shock
# normal periodic.
$ periodicDirections = [1,2]
$ lowerBounds = [XMIN,  YMIN,  ZMIN]
$ upperBounds = [XMAX,  YMAX,  ZMAX]
$ endif
$ endif

# Parameters to specify the fluid state at t=0.0
$ DENSITY_L = $1.0*REFERENCE_DENSITY$
$ DENSITY_R = $0.125*REFERENCE_DENSITY$
$ PRESSURE_L = $1.0*REFERENCE_PRESSURE$
$ PRESSURE_R =  $0.1*REFERENCE_PRESSURE$
$ NORMAL_VELOCITY_L = 0.0
$ NORMAL_VELOCITY_R = 0.0
$ PERPENDICULAR_VELOCITY_L = 0.0
$ PERPENDICULAR_VELOCITY_R = 0.0
$ TANGENTIAL_VELOCITY_L = 0.0
$ TANGENTIAL_VELOCITY_R = 0.0
$ NORMAL_FIELD_L = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ NORMAL_FIELD_R = $0.75*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_L = $1.0*math.sqrt(REFERENCE_DENSITY)$
$ PERPENDICULAR_FIELD_R = $-1.0*math.sqrt(REFERENCE_DENSITY)$
$ TANGENTIAL_FIELD_L = 0.0
$ TANGENTIAL_FIELD_R = 0.0

# Initialize a USim simulation
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)

# Setup the grid
addGrid(lowerBounds, upperBounds, numCells, periodicDirections)

# Create data structures needed for the simulation
createFluidSimulation()

# Step 1: Add Variables
addVariable(pi_value,math.pi)
addVariable(gas_gamma,GAS_GAMMA)
addVariable(densityL,DENSITY_L)
addVariable(densityR,DENSITY_R)
addVariable(pressureL,PRESSURE_L)
addVariable(pressureR,PRESSURE_R)
addVariable(normalVelocityL,NORMAL_VELOCITY_L)
addVariable(normalVelocityR,NORMAL_VELOCITY_R)
addVariable(perpendicularVelocityL,PERPENDICULAR_VELOCITY_L)
addVariable(perpendicularVelocityR,PERPENDICULAR_VELOCITY_R)
addVariable(tangentialVelocityL,TANGENTIAL_VELOCITY_L)
addVariable(tangentialVelocityR,TANGENTIAL_VELOCITY_R)
addVariable(mu0,MU0)
addVariable(normalFieldL,NORMAL_FIELD_L)
addVariable(normalFieldR,NORMAL_FIELD_R)
addVariable(perpendicularFieldL,PERPENDICULAR_FIELD_L)
addVariable(perpendicularFieldR,PERPENDICULAR_FIELD_R)
addVariable(tangentialFieldL,TANGENTIAL_FIELD_L)
addVariable(tangentialFieldR,TANGENTIAL_FIELD_R)

# Step 2: Add Pre-Expressions
 addPreExpression(rho = if (x>0.0, densityR,               densityL))
addPreExpression(pr  = if (x>0.0, pressureR,              pressureL))
addPreExpression(vx  = if (x>0.0, normalVelocityR,        normalVelocityL))
addPreExpression(vy  = if (x>0.0, perpendicularVelocityR, perpendicularVelocityL))
addPreExpression(vz  = if (x>0.0, tangentialVelocityR,    tangentialVelocityL))
addPreExpression(bx  = if (x>0.0, normalFieldR,           normalFieldL))
addPreExpression(by  = if (x>0.0, perpendicularFieldR,    perpendicularFieldL))
addPreExpression(bz  = if (x>0.0, tangentialFieldR,       tangentialFieldL))
addPreExpression(psi = 0.0)

# Step 3: Add expressions specifying initial condition on density,
# momentum, total energy, magnetic field and correction potential
addExpression(rho)
addExpression(rho*vx)
addExpression(rho*vy)
addExpression(rho*vz)
addExpression(pr/(gas_gamma-1.0) + 0.5*rho*(vx*vx+vy*vy+vz*vz)+0.5*((bx*bx+by*by+bz*bz)/mu0))
addExpression(bx)
addExpression(by)
addExpression(bz)
addExpression(psi)

# Add the spatial discretization of the fluxes
finiteVolumeScheme(DIFFUSIVE)

# Boundary conditions
boundaryCondition(copy,left)
boundaryCondition(copy,right)
boundaryCondition(periodic)

# Time integration
timeAdvance(TIME_ORDER)

# Run the simulation!
runFluidSimulation()

Note

For more depth, you can view the actual input blocks to Ulixes in the Setup window by choosing Save And Process Setup and then clicking on the shockTube.in file. In the .in file all macros are expanded to produce input blocks.

An Example Simulation

The input file for the Shock Tube (shockTube.pre) example for the USimBase package with SHOCK_TUBE = “BRIOWU” demonstrates each of the concepts described above to evolve the classic Brio & Wu shock tube problem in one-dimensional magnetohydrodynamics. Executing this input file within USimComposer and switching to the Visualize tab yields the plots shown in Fig. 3.

image 1

Figure 3: Visualization tab in USimComposer after executing the Brio & Wu shock tube input file for the tutorial.

Combining Euler and MHD Equations in the Same Input File

In Using USim to solve the Euler Equations we saw how most USimBase simulations for the Euler equations have an underlying pattern. Based on what we have discussed for the MHD equations, we can now extend this pattern to easily switch between the Euler and MHD equations:

# Are we solving the MHD equations?
$ MHD = True

# Import macros to setup simulation
$ import fluidsBase.mac
# if MHD
$ import idealmhd.mac
$ else
$ import euler.mac
$ endif

# Specify parameters for the specific physics problem
$ PARAM_1 = <value>
$ PARAM_2 = <value>
$ PARAM_N = <value>

# Initialize a USim simulation
$ if MHD
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,MU0,WRITE_RESTART,DEBUG)
$ else
initializeFluidSimulation(NDIM,0.0,END_TIME,NUMDUMPS,CFL,GAS_GAMMA,WRITE_RESTART,DEBUG)
$ endif

# Setup the grid
addGrid(lowerBounds, upperBounds, numCells, periodicDirections)

# Create data structures needed for the simulation
createFluidSimulation()

# Specify initial condition
# Step 1: Add Variables
addVariable(NAME,<value>)

# Step 2: Add Pre-Expressions
addPreExpression(<PreExpression>)

# Step 3: a) Add expressions specifying initial condition on density,
# momentum
addExpression(<expression>)
$ if MHD
# Step 3: b) Add expression specifying initial conditions on total
# energy, magnetic field, correction potential
addExpression(<expression>)
$ else
# Step 3: b) Add expression specifying initial conditions on total
# energy
addExpression(<expression>)
$ endif

# Add the spatial discretization of the fluxes
finiteVolumeScheme(DIFFUSIVE)

# Boundary conditions
boundaryCondition(<boundaryCondition,entity>)

# Time integration
timeAdvance(TIME_ORDER)

# Run the simulation!
runFluidSimulation()

You can see this pattern demonstrated in the Shock Tube (shockTube.pre) example.