Commit 2843b8e3 authored by Mladen Ivkovic's avatar Mladen Ivkovic
Browse files

added gizmo tex stuff

parent 1d3dc690
......@@ -185,6 +185,7 @@ theory/Multipoles/potential_short.pdf
theory/Multipoles/force_short.pdf
theory/Cosmology/cosmology.pdf
theory/Cooling/eagle_cooling.pdf
theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.pdf
m4/libtool.m4
m4/ltoptions.m4
......
\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,
\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}
\ No newline at end of file
%=========================================
\section{Explicit Computations}
%=========================================
Main purpose of this section is to have the explicit computations written down clearly somewhere.
Implicitly swallowing indices in the formulae can make programming life a nightmare, so here we go.
%=========================================
\subsection{Normalization}
%=========================================
To compute the normalisations \ref{omega}, we need to sum over all neighbouring particles and sum the kernels correctly.
To evaluate the kernels, we use the \verb|kernel_deval(xij, wij, wij_dx)| function in SWIFT.
If a kernel is defined as
\begin{align*}
W_i(\x) = W(\x - \x_i, h(\x)) = \frac{1}{h(\x)^\nu} w\left(\frac{| \x - \x_i |}{h(\x)} \right)
\end{align*}
then \verb|kernel_deval| computes
\begin{align*}
\mathtt{wij} &= w(\mathtt{xij}) \\
\text{and \quad} \mathtt{wij\_dx} &= \frac{\del w(r)}{\del r} \big{|} _{r = \mathtt{xij}} \\
\text{with \quad} \mathtt{xij} &= \frac{| \x_i - \x_j |}{h(\x_i)}
\end{align*}
So for a specific particle position $i$, we need to compute
\begin{align*}
\omega(\x_i) &= \sum_j W(\x_i - \x_j, h(\x_i)) \\
&= \sum_j \frac{1}{h(\x_i)^\nu} w\left( \frac{| \x_i - \x_j |}{h(\x_i)} \right) \\
&= \sum_j \frac{1}{h_i^\nu} \ \mathtt{ wij }
\end{align*}
with $h_i = h(\x_i)$.
%================================================
\subsection{Analytical gradients of $\psi(\x)$}
%================================================
For the \cite{ivanovaCommonEnvelopeEvolution2013} expression of the effective surfaces \Aij, we need analytical gradients in Cartesian coordinates of $\psi_i(\x_j)$.
From eq. \ref{psi} we have that
\begin{align*}
\psi_j(\x_i) &= \frac{ W(\x_i - \x_j, h(\x_i))}{\omega(\x_i)}\\
\end{align*}
Let $r_{ij} \equiv |\x_i - \x_j|$ and $q_{ij} \equiv \frac{r_{ij}}{h_i}$. Then
\begin{align}
\deldx \psi_j(\x_i) &= \deldx \frac{ W(\x_i - \x_j, h(\x_i))}{\omega(\x_i)}
= \deldx \frac{ W(r_{ij}, h_i)}{\omega(\x_i)} \nonumber \\
%
&= \frac{
\DELDX{W}(r_{ij}, h_i) \ \omega(\x_i) -
W(r_{ij}, h_i) \DELDX{\omega}(\x_i)
}{\omega(\x_i)^2} \nonumber \\
%
&= \frac{1}{\omega(\x_i)} \DELDX{W}(r_{ij}, h_i) -
\frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \DELDX{}\sum_k W(r_{ik}, h_i) \nonumber \\
%
&= \frac{1}{\omega(\x_i)} \DELDX{W}(r_{ij}, h_i) -
\frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \sum_k \DELDX{W}(r_{ik}, h_i) \label{grad_psi}
\end{align}
If a kernel\footnote{
Helpful way to think of the indices: $W_j(\x_i) = W(|\x_j - \x_i|, h_i)$ is the kernel value of particle $i$ at position $\x_i$ evaluated at the position $\x_j$.
Personally I would've used the indices the other way around for simplified thinking, but I'm going to stick to this notation because Hopkins also uses it.
}
is defined as
\begin{align*}
W_j(\x_i) = W(\x_i - \x_j, h_i))
= \frac{1}{h_i^\nu} w \left(\frac{| \x_i - \x_j |}{h_i} \right)
= \frac{1}{h_i^\nu} w ( q_{ij} )
\end{align*}
and we assume that the smoothing length $h_i$ is treated as constant at this point, then the gradient of the kernel is given by
\begin{align}
\deldx W_j (\x_i) &= \deldx \left( \frac{1}{h_i^\nu} w ( q_{ij} ) \right) \nonumber \\
&= \frac{1}{h_i^\nu} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\del q_{ij}(r_{ij})}{\del r_{ij}} \frac{\del r_{ij}}{\del x} \label{grad_kernel}
\end{align}
We now use
\begin{align}
\frac{\del q_{ij}(r_{ij})}{\del r_{ij}} &= \frac{\del}{\del r_{ij}} \frac{r_{ij}}{h_i} = \frac{1}{h_i} \label{dqdr} \\
\DELDX{r_{ij}} &= \deldx \sqrt{(\x_i - \x_j)^2}
= \frac{1}{2} \frac{1}{\sqrt{(\x_i - \x_j)^2}} \cdot 2 (\x_i - \x_j) \nonumber \\
&= \frac{\x_i - \x_j}{r_{ij}} \label{drdx}
\end{align}
To be perfectly clear, we should in fact write
\begin{align*}
r_j(\x) & \equiv | \x - \x_j | \\
\end{align*}
which again leads to
\begin{align*}
\DELDX{r_{ij}} &= \DELDX{r_j(\x_i)} = \DELDX{r_j(\x)} \bigg{|}_{\x = \x_i} \\
&= \deldx \sqrt{(\x - \x_j)^2} \big{|}_{\x = \x_i}
= \frac{1}{2} \frac{1}{\sqrt{(\x - \x_j)^2}} \cdot 2 (\x - \x_j) \big{|}_{\x = \x_i} \\
&= \frac{\x - \x_j}{r_{j}(\x)} \big{|}_{\x = \x_i}
= \frac{\x_i - \x_j}{r_{ij}}
\end{align*}
Inserting expressions \ref{dqdr} and \ref{drdx} in \ref{grad_kernel}, we obtain
\begin{equation}
\deldx W_j (\x_i) = \frac{1}{h_i^{\nu + 1}} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\x_i - \x_j}{r_{ij}} \label{grad_kernel_final}
\end{equation}
$\frac{\del w(q_{ij})}{\del q_{ij}}$ is given by \verb|wij_dx| of \verb|kernel_deval|.
Finally, inserting \ref{grad_kernel_final} in \ref{grad_psi} we get
\begin{align}
\deldx \psi_j (\x_i) &=
\frac{1}{\omega(\x_i)} \frac{1}{h_i^{\nu + 1}} \deldr{w}(r_{ij}, h_i) \frac{\x_i - \x_j}{r_{ij}} -
\frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \sum_k \frac{1}{h_i^{\nu + 1}} \deldr{w}(r_{ik}, h_i) \frac{\x_i - \x_k}{r_{ik}} \nonumber \\
&=
\frac{1}{\omega(\x_i)} \frac{1}{h_i^{\nu + 1}} \mathtt{wij\_dx} \frac{\x_i - \x_j}{r_{ij}} -
\frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \frac{1}{h_i^{\nu + 1}} \sum_k \mathtt{wik\_dx} \frac{\x_i - \x_k}{r_{ik}}
\end{align}
The definition of $r_{ij}$ requires a bit more discussion.
Since kernels used in hydrodynamics (at least in those methods currently implemented in SWIFT) are usually taken to be spherically symmetric, we might as well have defined
\begin{align*}
r'_{ij} = |\x_j - \x_i |
\end{align*}
which would leave the evaluation of the kernels invariant [$r'_{ij} = r_{ij}$], but the gradients would have the opposite direction:
\begin{align*}
\DELDX{r'_{ij}}
= \frac{\x_j - \x_i}{r'_{ij}}
= - \frac{\x_i - \x_j}{r_{ij}}
= - \DELDX{r_{ij}}
\end{align*}
So which definition should we take? \\
Consider a one-dimensional case where we choose two particles $i$ and $j$ such that $x_j > x_i$ and $q_{ij} = | x_j - x_i | / h_i < H$, where $H$ is the compact support radius of the kernel of choice.
Because we're considering a one-dimensional case with $x_j > x_i$, we can now perform a simple translation such that particle $i$ is at the origin, i.e. $x'_i = 0$ and $x'_j = x_j - x_i = | x_j - x_i | = r_{ij}$.
In this scenario, the gradient in Cartesian coordinates and in spherical coordinates must be the same:
\begin{align}
\frac{\del}{\del x'} W (|x_i' - x'|, h_i) \big{|}_{x' = x_j' } &= \frac{\del}{\del r_{ij}} W (r_{ij}, h_i) \nonumber\\
%
\Rightarrow \quad \frac{1}{h_i^\nu} \frac{\del w(q'_{ij})}{\del q'_{ij}} \frac{\del q'_{ij}(r'_{ij})}{\del r'_{ij}} \frac{\del r'_{i}(x')}{\del x'} \big{|}_{x' = x_j' } &=
\frac{1}{h_i^\nu} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\del q_{ij}(r_{ij})}{\del r_{ij}} \label{spher_cart}
\end{align}
We have the trivial case where
\begin{align*}
r'_{ij} &= |x_i' - x_j'| = | x_i - x_j | = r_{ij} \\
q'_{ij} &= r'_{ij}/h_i = q_{ij} \\
\Rightarrow \quad \frac{\del w(q'_{ij})}{\del q'_{ij}} &= \frac{\del w(q_{ij})}{\del q_{ij}}, \quad \frac{\del q'_{ij}(r'_{ij})}{\del r'_{ij}} = \frac{\del q_{ij}(r_{ij})}{\del r_{ij}}
\end{align*}
giving us the condition from \ref{spher_cart}:
\begin{equation}
\frac{\del r'_{i}(x')}{\del x'} \big{|}_{x' = x_j' } = 1 = \frac{\del r_{i}(x)}{\del x}
\end{equation}
this is satisfied for
\begin{align*}
r_j(\x) &= | \x - \x_j | \\
\Rightarrow r_{ij} &= | \x_i - \x_j |,\ \text{ not } \ r_{ij} = | \x_j - \x_i |
\end{align*}
\section{Equations}
\subsection{Partition of Unity and Related Quantities}
The partition of unity is defined as:
\begin{align}
\psi_i(\x) &= \frac{1}{\omega(\x)} W(\x - \x_i, h(\x)) \label{psi} \\
\omega(\x) &= \sum_j W(\x - \x_j, h(\x)) \label{omega}
\end{align}
where $h(\x)$ is some ``kernel size'' and $\omega(\x)$ is used to normalise the volume partition at any point $\x$.
We use the compact support radius $H$ of the kernels as $h$.
It can be shown that (provided $W(\x)$ is normalized, i.e. $\int_V W(\x) \D V = 1$):
\begin{align}
V_i &= \int_V \psi_i(\x) \de V = \frac{1}{\omega(\x_i)} \\ \label{psi_volume_integral}
V &= \sum_i V_i \\
\int_V f(\x) \D V &= \sum_i f(\x_i) V_i + \mathcal{O}(h^2) \label{vol_integral_fx}
\end{align}
Eq. (\ref{vol_integral_fx}) can be derived by Taylor-expanding $f(\x)$ around $\x_i$ and integrating the first order term by parts to show that it is zero.
\subsection{Meshless Hydrodynamics \`a la Hopkins}
Following \cite{hopkinsGIZMONewClass2015}, we arrive at the equation
\begin{equation}
\frac{\D}{\D t} (V_i \U_{k,i}) + \sum_j \F_{k,ij} \cdot \Aijm = 0
\end{equation}
with
\begin{equation}
\Aijm^\alpha = V_i \psitilde_j^\alpha (\x_i) - V_j\psitilde_i^\alpha (\x_j) \label{Hopkins}
\end{equation}
for every component $k$ of the Euler equations and every gradient component $\alpha$
The $\psitilde(\x)$ come from the $\mathcal{O}(h^2)$ accurate discrete gradient expression from \cite{lansonRenormalizedMeshfreeSchemes2008}:
\begin{align}
\frac{\del}{\del x_{\alpha}} f(\x) \big{|}_{\x_i} &=
\sum_j \left( f(\x_j) - f(\x_i) \right) \psitilde_j^\alpha (\x_i) \\ \label{gradient}
\psitilde_j^\alpha (\x_i) &= \sum_{\beta = 1}^{\beta=\nu} \mathbf{B}_i^{\alpha \beta}
(\x_j - \x_i)^\beta \psi_j(\x_i) \\
\mathbf{B}_i &= \mathbf{E_i} ^ {-1} \\
\mathbf{E}_i^{\alpha \beta} &= \sum_j (\x_j - \x_i)^\alpha (\x_j - \x_i)^\beta \psi_j(\x_i)
\end{align}
where $\alpha$ and $\beta$ again represent the coordinate components for $\nu$ dimensions.
\subsection{Meshless Hydrodynamics \`a la Ivanova}
Following \cite{ivanovaCommonEnvelopeEvolution2013}, we arrive at the equation
\begin{equation}
\frac{\D}{\D t} (V_i \U_{k,i}) + \sum_j \F_{k,ij} \cdot \Aijm = 0
\end{equation}
In the paper, no discrete expression for \Aij is given, instead:
\begin{equation}
\Aijm = \int_V \left[ \psi_j(\x) \nabla \psi_i(\x) - \psi_i(\x) \nabla \psi_j(\x) \right] \D V
\end{equation}
By Taylor-expanding the gradients of $\psi$ the volume integral can be discretised as:
\begin{equation}
\Aijm^\alpha = V_i \nabla^\alpha \psi_j (\x_i) - V_j \nabla^\alpha \psi_i (\x_j) + \mathcal{O}(h^2) \label{Ivanova}
\end{equation}
This time we get analytical gradients of $\psi$ instead of $\psitilde$.
Note that the indices $i$ and $j$ in these equations are switched compared to the one given in \cite{ivanovaCommonEnvelopeEvolution2013} so that the evolution equations for Hopkins and Ivanova versions are identical.
\include{./header}
%------------------------------------------
%: Metadata
%------------------------------------------
\title{GIZMO Implementation Details}
%\subtitle{Subtitle}
\author{Bert Vandenbroucke (bert.vandenbroucke@ugent.be), Mladen Ivkovic (mladen.ivkovic@epfl.ch)}
\date{}
%------------------------------------------
\newcommand{\Aij}{$\mathbf{A}_{ij}$}
\newcommand{\Aijm}{\mathbf{A}_{ij}} % A_ij math
\newcommand{\U}{\mathbf{U}}
\newcommand{\F}{\mathbf{F}}
\newcommand{\psitilde}{\tilde{\boldsymbol{\psi}}}
\begin{document}
%--------------------------------------------
% Stuff that needs to be done before all else
%--------------------------------------------
%\pagestyle{plain}
%\nocite{*} % show all entries of bibliography, even if they are not cited.
%Titlepage
\maketitle
\input{bert}
\input{equations}
\input{computations}
%------------
%Bibliography
%------------
\bibliography{references}
\end{document}
\documentclass[12pt, a4paper, english, singlespacing, parskip]{scrartcl}