Diffusion Processes

There are several general categories of diffusion processes that Athena++ can model:

• Hydrodynamic diffusion such as fluid viscosity and thermal conduction
• Field diffusion such as Ohmic dissipation and ambipolar diffusion
• Passive Scalars diffusion of the dye concentration

Each specific diffusion process has been natively implemented in the code.

Configuration

Users do not need to manually enable the diffusion solver during the configuration stage. These capabilities are automatically invoked when the runtime parameters specify nonzero coefficients.

However, certain diffusion processes depend on specific physics and will be ignored if the code is configured without the requisite options. For example, in order to calculate Ohmic dissipation and ambipolar diffusion, the magnetic field has to be enabled during configuration with -b. Similarly, thermal conduction is calculated only when a non-barotropic equation of state is enabled when Configuring the code. And passive scalar diffusion will only be enabled if --nscalars=1 or greater is used when configuring the solver.

Input File

The additional microphysics are enrolled automatically through the following optional input parameters:

    <problem>
    nu_iso      = 0.01       # isotropic viscosity coefficient
    nu_aniso    = 0.0        # anisotropic viscosity coefficient
    kappa_iso   = 0.01       # isotropic thermal conduction coefficient
    kappa_aniso = 0.0        # anisotropic thermal conduction coefficient
    eta_ohm     = 0.01       # Ohmic resistivity coefficient
    eta_ad      = 0.0        # Ambipolar diffusion constant (R/T to diff coeff by factor |B|^2)
    nu_scalar_iso  = 0.01    # isotropic passive scalar diffusion coefficient

Positive coefficients invoke the calculation of the corresponding physical processes. Note, the anisotropic viscosity and thermal conduction are not currently implemented in Athena++. However, users may define such processes in the Problem Generators and enroll functions via the Mesh::EnrollViscosityCoefficient(), Mesh::EnrollConductionCoefficient(), Mesh::EnrollFieldDiffusivity() methods.

Algorithms

Currently all diffusion processes are calculated from the primitive variables and update the conservative quantities; the resulting diffusion fluxes are directly added to the total hydro and/or field fluxes in an unsplit fashion (no temporal operator splitting is used). In order to calculate the viscosity stress tensor, we make use of the covariant derivative of vectors, with the help of geometric scale factors and their derivatives. Please see Covariant Expressions for Vector and Tensor Operators for a summary of the general expressions and a table of the scale factors used in the code.

Units

Note, the kappa_iso and kappa_aniso coefficients correspond to diffusivities, not conductivities. Also note that the current implementation uses a dimensionless system of units in that the factor (m̅/k_B) is not included in calculating the temperature (instead, T=P/d is used). That is, the energy flux is set to - kappa*d d(P/d)/dx, and kappa must have dimensions L²/t. For an energy flux of the form -kappa dT/dx (e.g. with density-independent Spitzer conductivity), the value assigned to kappa must include the factor (k_B/m̅) and must include a factor 1/ d.

Ambipolar Diffusion Input Parameter

Note, the eta_ad constant does not correspond to a diffusion coefficient. At the time of writing, diffusion coefficient is automatically computed by Athena++ via the following code:

pfdif->etaB(FieldDiffusion::DiffProcess::ambipolar, k,j,i) = pfdif->eta_ad*SQR(bmag(k,j,i));

and is equal to the product of the input <problem>/eta_ad parameter and the magnitude of the magnetic field strength squared. The induction equation therefore takes the form:

$\frac{\partial \mathbf{B}}{\partial t} - \nabla \times \left[\left(\mathbf{v} \times \mathbf{B} \right) - \eta \mathbf{J} - \eta_\mathrm{AD} \mathbf{B} \times \left(\mathbf{J} \times \mathbf{B} \right) \right ] = 0$

where $D_{\eta, AD} = \eta_\mathrm{AD} \lvert \mathbf{B}^2 \rvert \equiv$ the ambipolar diffusion coefficient (with units [length]^2/[time]).

In the above, $\eta_\mathrm{AD}$ is the input parameter <problem>/eta_ad. See Issue #369 for more details and discussion. We note that the documentation for the predecessor Athena (C) code has incorrectly placed a negative sign in front of of eta_ad term in the induction equation. This mistake has propagated to some publications, e.g. Lam et al (2019).