PLUTO: Computation of the Right Hand Side

The right hand side is constructed by adding, in succession, 1-D contributions coming from different sweeps:

$\vec{U}^{n+1}_{ijk} - \vec{U}^n_{ijk} = \vec{\cal R}_{\{x_1\},i,j,k}^n + \vec{\cal R}_{\{x_2\},i,j,k}^n + \vec{\cal R}_{\{x_3\},i,j,k}^n$

It is called once for every direction at every integration stage. The single contribution, say for the $ x_1 $ direction, is computed in a quasi-conservative way as

$\vec{\cal R}_{\{x_1\},i,j,k} = -\frac{\Delta t}{\dvol_i} \left( {\cal A}_{i+\HALF}\vec{\cal F}_{i+\HALF} - {\cal A}_{i-\HALF}\vec{\cal F}_{i-\HALF}\right) + \Delta t\left(\vec{\cal S}^{geo}_i + \vec{\cal S}^{grav}_i + \vec{\cal S}^{8w}_i\right)$

where the indices $ j $ and $ k $ stay fixed, $\dvol_d$ is the cell volume, $\vec{\cal F}^{\cal A}$ is the Riemann solver flux (not including pressure terms), ${\cal A}$ is the interface area. The source terms $\vec{\cal S}^{geo}_i, \vec{\cal S}^{grav}_i, \vec{\cal S}^{8w}_d$ are, respectively, the geometrical, gravitational and Powell's source terms (only for the MHD and RMHD modules).

The pressure term is separately discretized as a gradient operator and therefore it does not contain area factors nor it appears in the geometrical source terms. For instance, while incrementing the component of $\vec{\cal R}$ corresponding to the normal momentum $m_{x_1}$ we add

${\cal R}_{\{x_1\},i,j,k}^{[m_{x_1}]} \quad -= \frac{\Delta t}{\Delta l_i} \left(p_{i+\HALF} - p_{i-\HALF}\right)$

where $\Delta l_i$ is a line element.

The computation of the right hand side closely reflects the nature of the divergence or gradient operators involved in the original equations. Indeed, for scalar quantities such as density or energy, one simply has

$\pd{q}{t} + \nabla\cdot\vec{F} = 0\quad\Longrightarrow\quad \left\{\begin{array}{c@{+}c@{+}c@{+}c@{=}cl} \DS \pd{q}{t} & \DS \pd{F_x}{x} & \DS \pd{F_y}{y} & \DS \pd{F_z}{z} & \DS\; 0 &\qquad\textrm{Cartesian} \\ \noalign{\medskip} \DS \pd{q}{t} & \DS \frac{1}{r}\pd{(rF_{r})}{r} & \DS \frac{1}{r}\pd{F_{\phi}}{\phi} & \DS \pd{F_z}{z} & \DS \;0 & \qquad\textrm{Polar} \\ \noalign{\medskip} \DS \pd{q}{t} & \DS \frac{1}{r^2}\pd{(r^2F_{r})}{r} & \DS \frac{1}{rs}\pd{(s F_{\theta})}{\theta} & \DS \frac{1}{rs}\pd{F_\phi}{\phi} & \DS\; 0 &\qquad\textrm{Spherical} \\ \noalign{\medskip} \end{array}\right.$

For vector quantities such as momentum and magnetic field, we exploit the symmetric or antisymmetric properties of the corresponding flux tensor.

In Cartesian coordinates, vector components are treated as scalar quantities.
In polar or cylindrical coordinates, for the symmetric tensor $M_{ij} = M_{ji}$ (as in the momentum equation), we use the angular momentum-conserving form:
$\pd{\vec{m}}{t} + \nabla\cdot\vec{M} = 0 \quad \Longrightarrow \quad \left\{ \begin{array}{c@{+}c@{+}c@{+}c@{=}c} \DS \pd{m_{R}}{t} & \DS \frac{1}{R}\pd{(RM_{RR})}{R} & \DS \frac{1}{R}\pd{M_{\phi R}}{\phi} & \DS \pd{ M_{zR}}{z} & \DS \frac{M_{\phi\phi}}{R} \\ \noalign{\medskip} \DS \pd{m_{\phi}}{t} & \DS \frac{1}{R^2}\pd{(R^2M_{R \phi})}{R} & \DS \frac{1}{R} \pd{ M_{\phi \phi}}{\phi} & \DS \pd{ M_{z \phi}}{z} & \DS 0 \\ \noalign{\medskip} \DS \pd{m_z}{t} & \DS \frac{1}{R}\pd{(R M_{Rz})}{R} & \DS \frac{1}{R}\pd{M_{\phi z}}{\phi} & \DS \pd{M_{zz}}{z} & \DS 0 \end{array}\right.$
If the tensor is antisymmetric, $\Omega_{ij} = -\Omega_{ji}$ (as for the induction equation) the conservation law takes the form
$\pd{\vec{B}}{t} + \nabla\cdot\vec{\Omega} = 0 \quad \Longrightarrow \quad \left\{ \begin{array}{c@{+}c@{+}c@{+}c@{=}c} \DS \pd{B_{R}}{t} & \DS & \DS \frac{1}{R}\pd{\Omega_{\phi R}}{\phi} & \DS \pd{\Omega_{z r}}{z} & \DS\; 0 \\ \noalign{\medskip} \DS \pd{B_{\phi}}{t} & \DS \pd{\Omega_{R\phi}}{R} & \DS & \DS \pd{\Omega_{z\phi}}{z} & \DS\; 0 \\ \noalign{\medskip} \DS \pd{B_z}{t} & \DS \frac{1}{R}\pd{(R \Omega_{Rz})}{R} & \DS \frac{1}{R}\pd{\Omega_{\phi z}}{\phi} & \DS & \DS\; 0 \end{array}\right.$
where $\vec{\Omega} = \vec{v}\vec{B} - \vec{B}\vec{v}$ .

In spherical coordinates, for a symmetric tensor we use

$\pd{\vec{m}}{t} + \nabla\cdot\vec{M} = 0 \quad \Longrightarrow \quad \left\{ \begin{array}{c@{+}c@{+}c@{+}c@{=}c} \DS \pd{m_{r}}{t} & \DS \frac{1}{r^2} \pd{(r^2M_{rr})}{r} & \DS \frac{1}{rs} \pd{(s M_{\theta r})}{\theta} & \DS \frac{1}{rs} \pd{ M_{\phi r}}{\phi} & \DS \DS\frac{M_{\theta\theta} + M_{\phi\phi}}{r} \\ \noalign{\medskip} \DS \pd{m_{\theta}}{t} & \DS \frac{1}{r^2}\pd{(r^2M_{r\theta})}{r} & \DS \frac{1}{rs} \pd{(s M_{\theta\theta})}{\theta} & \DS \frac{1}{rs} \pd{ M_{\phi\theta}}{\phi} & \DS \frac{M_{\phi\phi}\cot\theta - M_{\theta r}}{r} \\ \noalign{\medskip} \DS \pd{m_{\phi}}{t} & \DS \frac{1}{r^3} \pd{(r^3 M_{r\phi})}{r} & \DS \frac{1}{rs^2}\pd{(s^2 M_{\theta\phi})}{\theta} & \DS \frac{1}{rs} \pd{ M_{\phi\phi}}{\phi} & \DS 0 \end{array}\right.$

(here $s\equiv\sin\theta$ ) whereas, if the tensor is antisymmetric, $\Omega_{ij} = -\Omega_{ji}$ , one has

$\pd{\vec{B}}{t} + \nabla\cdot\vec{\Omega} = 0 \; \Longrightarrow \; \left\{\begin{array}{c@{+}c@{+}c@{+}c@{=}c} \DS \pd{B_{r}}{t} & \DS & \DS \frac{1}{rs} \pd{(s \Omega_{\theta r})}{\theta} & \DS \frac{1}{rs} \pd{\Omega_{\phi r}}{\phi} & \DS 0 \\ \noalign{\medskip} \DS \pd{B_{\theta}}{t} & \DS \frac{1}{r} \pd{(r\Omega_{r\theta})}{r} & \DS & \DS \frac{1}{rs}\pd{\Omega_{\phi\theta}}{\phi} & \DS 0 \\ \noalign{\medskip} \DS \pd{B_{\phi}}{t} & \DS \frac{1}{r} \pd{(r\Omega_{r\phi})}{r} & \DS \frac{1}{r} \pd{ \Omega_{\theta\phi}}{\theta} & \DS & \DS 0 \\ \noalign{\medskip} \end{array}\right.$

The geometrical source terms (right hand sides in the previous equations) are evaluated by averaging interface values at the cell-centered.