Cold Relativistic Fluid Model

An alternate to the representation of the plasma as macroparticles, is the cold relativistic fluid model. The cold relativistic fluid model is typically faster to compute than particles and useful to replace particles when kinetic effects are negligible. This can be done for certain parts of the plasma, such as the background plasma, while still representing particle beams as a group of macroparticles. The two models then couple through Maxwell’s equations.

In the cold limit (zero internal temperature and pressure) of a relativistic plasma, the Maxwell-Fluid equations govern the plasma evolution. The fluid equations per species, s, are given by,

\[\begin{split}\frac{\partial N_s}{\partial t} + \nabla \cdot (N_s\mathbf{V}_s) &= 0 \\ \frac{\partial (N\mathbf{U})_s}{\partial t} + \nabla \cdot ((N\mathbf{U})_s\mathbf{V}_s) &= \frac{q_sN_s}{m_s}(\mathbf{E}_s + \mathbf{V}_s \times \mathbf{B}_s).\end{split}\]

Where the fields are updated via Maxwell’s equations,

\[\begin{split}\nabla \cdot \mathbf{E} &= \frac{\rho}{\varepsilon_0} \\ \nabla \cdot \mathbf{B} &= 0 \\ \nabla \times \mathbf{E} &= -\frac{\partial \mathbf{B}}{\partial t} \\ \nabla \times \mathbf{B} &= \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \frac{\partial \mathbf{E}}{\partial t}.\end{split}\]

The fluids are coupled to the fields through,

\[\begin{split}\rho &= \rho_{ptcl}+\sum_s q_sN_s \\ \mathbf{J} &= \mathbf{J}_{ptcl}+\sum_s q_sN_s\mathbf{V}_s \\ \mathbf{V}_s &= \frac{ \mathbf{U}_s }{ \sqrt{ 1 + \mathbf{U}_s^2/c^2} } \\ (N\mathbf{U})_s &= N_s\mathbf{U}_s\end{split}\]

where the particle quantities are calculated by the PIC algorithm.

Implementation details

Fig. 41 Fluid Loop embedded within the overall PIC loop.

The fluid timeloop is embedded inside the standard PIC timeloop and consists of the following steps: 1. Higuera and Cary push of the momentum 2. Non-inertial (momentum source) terms (only in cylindrical geometry) 3. boundary conditions and MPI Communications 4. MUSCL scheme for advection terms 5. Current and Charge Deposition. Fig. 41 gives a visual representation of these steps, and we describe each of these in more detail.

Step 0: Preparation

Before the fluid loop begins, it is assumed that the program is in the state where fields \(\mathbf{E}\) and \(\mathbf{B}\) are available integer timestep. The fluids themselves are represented by arrays of fluid quantities (density and momentum density, \(\mathbf{Q} \equiv \{ N, NU_x, NU_y, NU_z \}\)) known on a nodal grid and at half-integer timestep.

Step 1: Higuera and Cary Push

The time staggering of the fields is used by the momentum source term, which is solved with a Higuera and Cary push [1]. We do not adopt spatial grid staggering, all discretized fluid quantities exist on the nodal grid. External fields can be included at this step.

Step 2: Non-inertial Terms

In RZ, the divergence of the flux terms has additional non-zero elements outside of the derivatives. These terms are Strang split and are time integrated via equation 2.18 from Shu and Osher [2], which is the SSP-RK3 integrator.

Step 3: Boundary Conditions and Communications

At this point, the code applies boundary conditions (assuming Neumann boundary conditions for the fluid quantities) and exchanges guard cells between MPI ranks in preparation of derivative terms in the next step.

Step 4: Advective Push

For the advective term, a MUSCL scheme with a low-diffusion minmod slope limiting is used. We further simplify the conservative equations in terms of primitive variables, \(\{ N, U_x, U_y, U_z \}\). Which we found to be more stable than conservative variables for the MUSCL reconstruction. Details of the scheme can be found in Van Leer [3].

Step 5: Current and Charge Deposition

Once this series of steps is complete and the fluids have been evolved by an entire timestep, the current and charge is deposited onto the grid and added to the total current and charge densities.

Note

The algorithm is safe with zero fluid density.

It also implements a positivity limiter on the density to prevent negative density regions from forming.

There is currently no ability to perform azimuthal mode decomposition in RZ.

Mesh refinement is not supported for the fluids.

The implemented MUSCL scheme has a simplified slope averaging, see the extended writeup for details.

More details on the precise implementation are available here, WarpX_Cold_Rel_Fluids.pdf.

Warning

If using the fluid model with the Kinetic-Fluid Hybrid model or the electrostatic solver, there is a known issue that the fluids deposit at a half-timestep offset in the charge-density.

[1]

A. V. Higuera and J. R. Cary. Structure-preserving second-order integration of relativistic charged particle trajectories in electromagnetic fields. Physics of Plasmas, 24(5):052104, 04 2017. URL: https://doi.org/10.1063/1.4979989, arXiv:https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4979989/15988441/052104\_1\_online.pdf, doi:10.1063/1.4979989.

[2]

C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2):439–471, 1988. URL: https://www.sciencedirect.com/science/article/pii/0021999188901775, doi:https://doi.org/10.1016/0021-9991(88)90177-5.

[3]

B. Van Leer. On The Relation Between The Upwind-Differencing Schemes Of Godunov, Engquist—Osher and Roe, pages 33–52. Springer Berlin Heidelberg, 1997. URL: https://doi.org/10.1007/978-3-642-60543-7_3, doi:10.1007/978-3-642-60543-7_3.