Cosmic‐ray diffusion with Multigrid

This module solves anisotropic cosmic-ray diffusion with a fully-implicit time integrator. The code solves $$\frac{\partial e_{\rm CR}}{\partial t}+\nabla\cdot\left(e_{\rm CR}\mathbf{v}\right)=-p_{\rm CR}\nabla\cdot\mathbf{v}+\nabla\cdot\left(\mathbb{D}\nabla e_{\rm CR}\right)+Q_{\rm CR},$$ $$p_{\rm CR} = \frac{1}{3}e_{\rm CR},$$ $$Q_{\rm CR} = S_{\rm CR}-\Lambda_{\rm coll}n_{\rm H}e_{\rm CR},$$ $$\mathbb{D} = D_{ab}= D_\perp\delta_{ab}+\left(D_\parallel-D_\perp\right)n_an_b,$$ $$n_{a} = \frac{B_a}{|\mathbf{B}|}.$$ This code is not very general yet. Currently it supports only single energy, and constant diffusion and absorption coefficients. Also, the velocity-related terms are not yet implemented; the algorithm assumes such terms are insignificant. This module is compatible with static and adaptive mesh refinement.

In many cases, the CR distribution quickly approaches a steady state. The code solves the steady-state equation when <crdiffusion> steady = true is set in an input file.

This is a demonstration of extending the Multigrid solver to various physical processes.

Configuration and compilation

To enable the Multigrid CR diffusion sovler, simply configure the code with the -crdiff. src/pgen/cr_diffusion_mg.cpp (with inputs/cosmic_ray/athinput.cr_diffusion_mg) is a simple test problem.

> python configure.py -b -crdiff --prob cr_diffusion_mg
> make clean ; make

Optionally you can enable MPI parallelization with -mpi.

Input File

The parameters for the CR diffusion module must be set in the <crdiffusion> block.

<crdiffusion>
mgmode          = FMG      # FMG (Full Multigrid) or MGI (Multigrid Iteration)
fas             = true     # Use FAS method
npresmooth      = 2        # Number of pre-smoothing (0,1,2)
npostsmooth     = 2        # Number of post-smoothing (0,1,2)
omega           = 1.0      # Acceleration parameter
niteration      = 30       # Number of iteration
threshold       = 0.0      # Convergence threshold
output_defect   = true     # Output defect
show_defect     = false    # Show defect at every cycle
steady          = false    # Steady-state approximation
ix1_bc          = user     # inner-x1 boundary condition
ox1_bc          = user     # outer-x1 boundary condition
ix2_bc          = user     # inner-x2 boundary condition
ox2_bc          = user     # outer-x2 boundary condition
ix3_bc          = user     # inner-x3 boundary condition
ox3_bc          = user     # outer-x3 boundary condition

Dpara           = 1000     # Diffusion coefficient along B-field
Dperp           = 10       # Diffusion coefficient across B-field
Lambda          = 100      # Absorption coefficient

Similar to self-gravity, you can change the Multigrid parameters. However, the default values shown above are safe and reasonably good choice. For convergence criterion, you need to specify either niteration or threshold (when both are specified, only threshold is used). With the above configuration, typically you can obtain a convergence factor of ~ 0.7, but we strongly recommend users to test it for yourself and set an appropriate convergence criterion. See Notes below and also Self Gravity with Multigrid.

Initialization and boundary conditions

Unless the steady-state mode is in use, the initial CR energy density must be set in the ProblemGenerator function.

  for(int k=ks; k<=ke; ++k) {
    for(int j=js; j<=je; ++j) {
      for(int i=is; i<=ie; ++i)
        pcrdiff->ecr(k,j,i) = 1.0;
    }
  }

Optionally, the source term distribution $S_{\rm CR}$ can be set in pcrdiff->source(k,j,i). Also, appropriate boundary conditions must be set. Please consult the example problem generator file about the boundary conditions.

Using results in another module

One particular use case of the CR diffusion is calculating the ionization rates for non-ideal MHD effects or chemical reactions. In the case of the non-ideal MHD, one can define a function like below and enroll it using Mesh::EnrollFieldDiffusivity() in your problem generator. Note that the CR diffusion module is called at the beginning of each time step (not substep), so the CR energy density at the beginning of the time step is stored in pcrdiff->ecr.

void CalcDiffusivity(FieldDiffusion *pfdif, MeshBlock* pmb, const AthenaArray<Real> &w,
     const AthenaArray<Real> &bmag, int is, int ie, int js, int je, int ks, int ke) {
  Real factor = 1.0; // constant coefficient
  for (int k = ks; k <= ke; ++k) {
    for (int j = js; j <= je; ++j) {
      for (int i = is; i <= ie; ++i)
        Real zeta = factor * w(IDN, k, j, i) * pmb->pcrdiff->ecr(k, j, i);
        // calculate your diffusion rates using zeta
        // ...
    }
  }

  return;
}

If you want to output zeta, you can store them in user-defined output variables.

Notes

The code uses the V(2,2) cycle with the red-black Jacobi smoother by default, because the discretized stencil for the anisotropic diffusion equation is 19-points (center + 6 faces + 12 edges), the red-black Gauss-Seidel solver is not compatible. Also, because of this complex stencil, we simply adopt the trilinear prolongation and volume-weighted restriction methods at level boundaries instead of the flux-conserving discretization used in the Poisson solver for the self-gravity.

The Multigrid solver is quite robust and efficient, but strongly anisotropic problems are difficult. As far as we tested, Dpara:Dperp = 100:1 is feasible, with a convergence factor of ~ 0.7 using omega = 1.0. When it is a purely diffusion problem without the source term, using a larger acceleration parameter, e.g. omega = 1.2, can accelerate convergence. However, when the source term is not negligible, the acceleration parameter causes numerical instabilities and poor convergence. Therefore, we recommend to use omega = 1.0 unless there is a clear advantage to use other value.

Cosmic-ray transport module

Compared to the Cosmic Ray Transport module which solves up to the first moment with the reduced-speed-of-light technique, this module adopts the diffusion approximation solving the zeroth moment only, which is less accurate. On the other hand, as this is a fully-implicit method, the time step is not constrained by the stability condition and a very long time step (e.g. the time step for the MHD part) can be used. Also, it does not use the reduced speed of light technique, which may be inappropriate in some situations. As they behave quite differently, we recommend users to learn the difference and consider which module to use depending on the problem of interest.