Using USim to Solve an Anisotropic Diffusion Problem

In this tutorial we show how to use USim to solve a problem with anisotropic diffusion using the updater diffusion (1d, 2d, 3d) with option anisotropicDiffusion. The problem used an example is based off of the problem described in

Parrish, Ian J., and James M. Stone. "Nonlinear evolution of the magnetothermal instability in two dimensions." The Astrophysical Journal 633.1 (2005): 334.

but solved on an unstructured grid.

Required DataStructs

This tutorial follows many of the same concepts as prior tutorials

First of all we initialize a vector (in this case B) that defines the field that is used to construct the conductivity tensor

<DataStruct B>
  kind = nodalArray
  onGrid = domain
  numComponents = 3
  writeOut = 1
</DataStruct>

In addition we define the parallel conductivity (parallel to the vector field) kParallel

<DataStruct kParallel>
  kind = nodalArray
  onGrid = domain
  numComponents = 1
  writeOut = 1
</DataStruct>

And the conductivity perpendicular to the vector field kPerpendicular

<DataStruct kPerpendicular>
  kind = nodalArray
  onGrid = domain
  numComponents = 1
  writeOut = 1
</DataStruct>

Finally we define the conductivity tensor conductivityTensor. The conductivity tensor always has 9 components even if the simulation is 2 dimensional

<DataStruct conductivityTensor>
  kind = nodalArray
  onGrid = domain
  numComponents = 9
  writeOut = 1
</DataStruct>

Computing a Conductivity Tensor

In general the conductivityTensor could be filled up manually using a combiner (1d, 2d, 3d), but thermal diffusion in a magnetic field is best split into parallel and perpendicular conductivities. A built in source computes this automatically given the parallel and perpendicular conductivities. The conductivityTensor takes a vector field B a parallel conductivity kParallel (scalar) and a perpendicular conductivity kPerpendicular (scalar). The use of a Conductivity Tensor is given below

<Updater initConductivityTensor>
  kind = equation2d
  onGrid = domain

  in = [B, kParallel, kPerpendicular]

  out = [conductivityTensor]

  <Equation a>
    kind = conductivityTensor
  </Equation>

</Updater>

Solving problems using Derivatives with option anisotropicDiffusion

Now that conductivityTensor is defined we can compute the diffusion. We use the diffusion (1d, 2d, 3d) which supplies second order derivatives. In this case the we use anisotropicDiffusion which expects a scalar input (in this case temperature) and a 9 component conductivity tensor (in this case conductivityTensor)

<Updater computeDiffusion>
  kind = diffusion2d
  derivative = anisotropicDiffusion
  onGrid = domain
  numScalars = 3
  coefficient = 1.0

  orderAccuracy = 1
  numberOfInterpolationPoints = 8

  in = [temperature, conductivityTensor]
  out = [temperatureNew]
</Updater>

Note that the complete list of options available in diffusion (1d, 2d, 3d) are diffusion, anisotropicDiffusion and gradientOfDivergence.

Computing the time step for the diffusion operator

Time integration is performed using super time stepping. Super time stepping is a variable stage Runge-Kutta approach that is much faster (by the number of stages) than standard Runge-Kutta methods for solving diffusion problems. The approach requires two time steps. The desired (actual) time step is computed using a timeStepRestrictionUpdater (1d, 2d, 3d). The key here is that CFLSTEP=100.0 so it is much higher than the explicitly stable time step for a diffisive system

 <Updater timeStepRestriction>

  kind = timeStepRestrictionUpdater2d
  in = [maxConductivity]

  onGrid = domain
  restrictions = [quadratic]

  <TimeStepRestriction quadratic>
    kind = quadratic
    cfl = CFLSTEP
  </TimeStepRestriction>

</Updater>

Next the time step for the super time stepping method is computed where an explicitly stable CFL=0.1 is used. The time step is stored in diffDT1

 <Updater getDiffDT1>
  kind = getTimeStepUpdater2d
  in = [maxConductivity]
  out = [diffDT1]
  onGrid = domain
  restrictions = [quadratic]

  <TimeStepRestriction quadratic>
    kind = quadratic
    cfl = CFL
  </TimeStepRestriction>
</Updater>

The Super Time Stepping integrator then knows to take the ratio of the desired time step and the explicitly stable time step to compute the number of stages used in the STS Updater.

An Example Simulation

The input file for the problem Anisotropic Diffusion in the USimHEDP package demonstrates each of the concepts described above.