\section{GIZMO~: variables and basic equations}

This section is left unchanged and as the first in this document such that the comments in the code still correspond to the correct equations.


\subsection{Variables}
6 standard variables~: $\vec{x}$, $\vec{v}$\\  \\
5 conserved variables ($C$)~: $m$, $\vec{p}$, $E$\\ \\
5 primitive variables ($X$)~: $\rho$, $\vec{w}$, $P$\\
It could be that $\vec{w} = \vec{v}$, but this is not generally true for all methods...\\ \\
6 matrix elements + the volume (can be computed from $m$ and $\rho$ if this is
accurate enough)\\ \\
15 gradients (3 for every primitive variable)\\

In total: 37 quantities per particle (without the volume), 34 if the velocity is equal to the fluid velocity

\subsection{Equations by Neighbour Loops}

\subsubsection{Loop 1: volumes (=densities) + matrix elements}
Volumes are given by (combination of eqns. (7) and (27) from Hopkins)~:
\begin{equation}
V_i = \frac{1}{\sum_j W(|\vec{x_i}-\vec{x_j}|, h_i)}
\end{equation}

The 6 elements of the (symmetric) matrix $E_i$ by (eqn. (14) in Hopkins without the normalization, because $E_i$ is only used in combination with $\psi_i$ (eqn. (12)), which contains the same normalization and cancels it out again)~:
\begin{equation}
E_i^{\alpha \beta} = \sum_j (\vec{x_j}-\vec{x_i})^\alpha (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i)
\end{equation}

The kernel functions in SWIFT do not include the normalization factor $\frac{1}{h_i^3}$, so this is added after the loop in the ghost.

\subsubsection{Ghost 1: primitive variables + matrix inversion}
Use the volume to convert conserved variables to primitive variables (eqn. (5), (7) and (9) in the AREPO paper)~:
\begin{align}
\rho &= \frac{m}{V}\\
\vec{v} &= \frac{\vec{p}}{m}\\
P &= \frac{\gamma - 1}{V}\left(E - \frac{1}{2} \frac{|\vec{p}|^2}{m}\right)
\end{align}

Invert the matrix and call the inverse matrix $B_i$ (eqn. (13) in Hopkins).

\subsubsection{Loop 2: gradients}
Calculate gradients for the primitive variables (eqn. (12) in Hopkins with eqn. (6) inserted and the normalization constant eliminated)~:
\begin{equation}
	\left(\vec{\nabla} X_i\right)^\alpha = \sum_j \left(X_j - X_i\right) \sum_\beta B_i^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i)
\end{equation}
If you want to use a slope limiter of some sorts, then this would also be the time to collect the necessary data.

\subsubsection{Ghost 2: nothing?}
Perform the slope limiting if wanted.

\subsubsection{Loop 3: hydro}
For every neighbour $j$, calculate an interface area (combination of eqn. (12) in Hopkins and the definition of $\vec{A_{ij}}$, given in between eqns. (18) and (19))
\begin{equation}
	\left(\vec{A_{ij}}\right)^\alpha = V_i \sum_\beta B_i^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i) + V_j \sum_\beta B_j^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_j)
\end{equation}

Calculate a position for the interface (eqn. (20) in Hopkins)~:
\begin{equation}
	\vec{x_{ij}} = \vec{x_i} + \frac{h_i}{h_i+h_j} (\vec{x_j} - \vec{x_i})
\end{equation}
and a velocity (eqn. (21) in Hopkins)~:
\begin{equation}
	\vec{v_{ij}} = \vec{v_i} + (\vec{v_j}-\vec{v_i})\left[ \frac{(\vec{x_{ij}}-\vec{x_i}).(\vec{x_j}-\vec{x_i})}{|\vec{x_i}-\vec{x_j}|^2} \right]
\end{equation}

Use the velocity to boost the primitive variables $X_k$ ($k=i,j$) to the rest frame of the interface ($X'_k$). This comes down to applying a correction to the fluid velocities.

Use the gradients to interpolate the primitive variables from positions $\vec{x_i}$ and $\vec{x_j}$ to $\vec{x_{ij}}$ and predict them forward in time
by half a time step to obtain second order accuracy in time (eqn. (A3) in Hopkins)~:
\begin{equation}
	X''_k = X'_k + \left(\vec{\nabla}X\right)_k.(\vec{x_{ij}}-\vec{x_k}) + \frac{\Delta t}{2} \frac{\partial X'_k}{\partial t}.
\end{equation}
The time derivatives are given by (eqn. (A4) in Hopkins)
\begin{equation}
	\frac{\partial X'_k}{\partial t} = -\begin{pmatrix}
	\vec{w'} & \rho' & 0\\
	0 & \vec{w'} & 1/\rho'\\
	0 & \gamma P' & \vec{w'}
	\end{pmatrix} \vec{\nabla}X_k, \label{eq:euler_primitive}
\end{equation}
or, for example~:
\begin{equation}
	\frac{\rho_k^{n+1/2}-\rho_k^n}{\Delta t/2} = -\left( \vec{w_k}-\vec{v_{ij}}^n \right).\vec{\nabla}\rho_k^n - \rho_k^n \vec{\nabla}.\vec{w_k}^n
\end{equation}

We then feed the $X''_k$ to a 1D Riemann solver, that gives us a vector of fluxes $\vec{F''_{ij}}$, which we deboost to a static frame of reference ($\vec{F_{ij}}$). To account for the arbitrary orientation of the interface, we also feed the normal vector to the interface to the Riemann solver. The Riemann solver internally uses the fluid velocity along the interface normal to solve an effective 1D Riemann problem. The velocity solution along the interface normal is also internally added to the velocities (see GIZMO source code). Mathematically, this is the same as first rotating the velocities to a frame aligned with the interface and then rotating the solution back to the original frame (e.g. AREPO paper), but it is computationally cheaper and causes less round off error in the solution.

Given the solution vector $X_\text{half}$ returned by the Riemann problem, the fluxes are given by (eqns. (1)-(3) in Hopkins)
\begin{align}
	\vec{F}_\rho &= 
		\rho_\text{half} \left( \vec{w}_\text{half} - \vec{v}_{ij}\right)\\
	\vec{F}_{w_k} &= 
		\rho_\text{half} w_{k, \text{half}} \left( \vec{w}_\text{half} - \vec{v}_{ij} \right) + P_\text{half} \vec{n}_k\\
	\vec{F}_P &= 
		\rho_\text{half} e_\text{half} \left( \vec{w}_\text{half} - \vec{v}_{ij}\right) + P_\text{half} \vec{w}_\text{half},
\end{align}
with $\vec{n}_k$ the unit vector along the coordinate axes ($k=x,y,z$) and $e = \frac{P_\text{half}}{\left(\gamma -1\right)\rho_\text{half}} + \frac{1}{2}\vec{w}_\text{half}^2$.

Finally, we use the fluxes to update the conserved variables (eqn. (23) in Hopkins)~:
\begin{equation}
	\Delta C_k = -\Delta t \sum_j \vec{F_{ij}}.\vec{A_{ij}}.
\end{equation}

\subsection{Other equations that are used}

In order to couple the GIZMO schemes to subgrid physics, we need to provide a way to obtain the time 
derivative of the internal energy, which is not a basic quantity in these schemes. This time derivative can be 
obtained from the Euler equations, using the gradients. From \eqref{eq:euler_primitive}, we get

\begin{equation}
\frac{\partial{}P}{\partial{}t} = -\gamma{}P \vec{\nabla{}}.\vec{w} - \vec{w}.\vec{\nabla{}}P,
\end{equation}
while we also have
\begin{equation}
\frac{\partial{}P}{\partial{}t} = \frac{\partial{}}{\partial{}t} \left(
  (\gamma{}-1) \rho{} u
\right) = (\gamma{}-1)\rho{}\frac{\partial{}u}{\partial{}t} + (\gamma{}-1)u\frac{\partial{}\rho{}}{\partial{}t}
\end{equation}
or hence
\begin{equation}
\frac{\partial{}u}{\partial{}t} = \frac{1}{(\gamma{}-1)\rho{}} \frac{\partial{}P}{\partial{}t} -
\frac{u}{\rho{}} \left( -\vec{w}.\vec{\nabla{}}\rho{} - \rho{}\vec{\nabla{}}.\vec{w} \right),
\end{equation}
where we again used \eqref{eq:euler_primitive} in the last term.

Combining these equations, we get
\begin{align}
\frac{\partial{}u}{\partial{}t} &= \frac{1}{(\gamma{}-1)\rho{}} \left(
-\gamma{}P \vec{\nabla{}}.\vec{w} - \vec{w}.\vec{\nabla{}}P
\right) -
\frac{u}{\rho{}} \left( -\vec{w}.\vec{\nabla{}}\rho{} - \rho{}\vec{\nabla{}}.\vec{w} \right) \\
&= -\gamma{}u \vec{\nabla{}}.\vec{w} - \frac{1}{(\gamma{}-1)\rho{}} \vec{w}.\vec{\nabla{}}P
+ \frac{u}{\rho{}}\vec{w}.\vec{\nabla{}}\rho{} + u \vec{\nabla{}}.\vec{w} \\
&= -(\gamma{}-1)u \vec{\nabla{}}.\vec{w} - \frac{1}{(\gamma{}-1)\rho{}} \vec{w}.\vec{\nabla{}}P
+ \frac{u}{\rho{}}\vec{w}.\vec{\nabla{}}\rho{}.
\end{align}

For convenience, we can eliminate $u$ from this equation again:
\begin{equation}
\frac{\partial{}u}{\partial{}t} = -\frac{P}{\rho{}}\vec{\nabla{}}.\vec{w}
- \frac{1}{(\gamma{}-1)\rho{}} \vec{w} . \left(
\vec{\nabla{}}P - \frac{P}{\rho{}} \vec{\nabla{}}\rho{}
\right),
\end{equation}
where we recognise a factor $\vec{\nabla{}}u$ in the second term.