Philipe's Notes on SPH algorithm

Page last edited 3,347 days ago
From Wei-Liang Qian
Jump to navigation Jump to search

%\documentclass[a4paper, twocolumn]{article} \documentclass[a4paper]{article}

\usepackage[english]{babel} \usepackage[utf8x]{inputenc} %\usepackage{amsmath} \usepackage{graphicx} \usepackage[colorinlistoftodos]{todonotes} \usepackage{color} \usepackage{grffile}

\usepackage{nath} \delimgrowth=1 \mathindent=10pt


\title{Notes on SPH algorithm} \author{Ph. Mota}

\def\abs#1{\left|#1\right|} \def\ave#1{\langle#1\rangle} \def\aave#1{\lAngle#1\rAngle} \def\dTemplate:\rm d \def\nq{\\\quad} \def\nqq{\\\qquad}

\begin{document} \maketitle

\begin{abstract} We review Smoothed Paticle Dynamics (SPH) and discuss the precision, the stability as well as the efficiency of the algorithm. \end{abstract}



\section{The algorithm} In order to solve the hydrodynamic equations we apply the following steps.



\section{Time evolution} The hydrodynamic equations (the energy-momentum tensor conservation and the currents conservations) $$ T^{\mu\nu}{}_{;\mu} = 0 \\ N_A^\mu{}_{;\mu} = 0 $$ are reexpressed as follows. First, we write a reference density $\rho$ at point $a$ which is related to a conserved quantity $\nu$ of test particles $$ \rho_a = \sum_b \nu_b W_{ab} $$ which readily satisfies the following conservation equation $$ \frac{\partial \rho_a}{\partial t}+\nabla \cdot \left(\rho_a {\bf v_a}\right) = 0 $$ or $$ \partial_\mu(\rho_a \dot{r}_a^\mu) = 0 $$

The part of the energy-momentum conservation parallel to the 4-velocity is can be rewritten as the conservation of the entropy density current $$ \frac{1}{\tau}\partial_\mu(\tau \gamma s v^\mu) = 0 $$ Now if we identify $\dot{r}_a^\mu = v_a^\mu$ we have (for $\tau\neq0$), $$ \partial_\mu(\frac{\tau_a \gamma_a s_a \rho_a v^\mu_a}{\rho_a}) = 0 \\ \rho_a v^\mu_a \partial_\mu(\frac{\tau_a \gamma_a s_a }{\rho_a}) = 0 $$ Then, for non-vanishing reference density, we can write $$ \frac{\d}{\d\tau}(\frac{ \tau_a \gamma_a s_a }{\rho_a}) = 0 $$

In fact, the above equation can be used to define the weights $\nu_a$ of the reference density. In principle, since $\nu$ is not related to any particular physics, they can all be set to unity. However, this makes $\rho$ reflect the spatial density of test particles which is not guaranteed to be a smooth distribution. By using the entropy density as a proxy for the reference density, we ensure that the reference density is a smooth function regardless of the test particle distribution. This provides some stability to the method.

In fact, the particle positions $r_a$ can be set so that the reference density (even with unit weights) follows the entropy density distribution. However, in a multidimensional setting, this is a non-linear problem and cannot be straightforwardly solved. Instead, we propose that the test particle positions are randomly sampled against ({\color{red}what does this mean?}) the entropy density distribution and then the weights are corrected to produce a new reference density $\rho'$ using a Taylor expansion of the entropy density profile (given as a initial condition)

In principle the actual equation to solve for the evolution of the entropy density is $$ v^\mu\partial_\mu(\tau \gamma s) + \tau \gamma s\partial_\mu v^\mu = 0 $$ or $$ \frac{\d}{\d\tau}(\tau \gamma s) = - \tau \gamma s\partial_i v^i $$


For the energy momentm conservation we have $$ (\omega u^\mu u^\nu )_{;\nu} - g^{\mu\nu}\partial_\nu P = 0 \\ (\frac{s \omega u^\mu u^\nu}{s} )_{;\nu} - g^{\mu\nu}\partial_\nu P = 0 \\ s u^\nu (\frac{ \omega u^\mu }{s} )_{;\nu} - g^{\mu\nu}\partial_\nu P = 0 $$

For the space part we have $$ s u^\nu \partial_\nu(\frac{ \omega u^x }{s} ) = -\partial_x P \\ s u^\nu \partial_\nu(\frac{ \tau^2 \omega u^\eta }{s} ) = -\partial_\eta P $$ or $$ s u^\nu \partial_\nu(\frac{ -g^{ij} \omega u^j }{s} ) = -\partial_i P $$

The evolution equations are therefore $$ \frac{\d}{\d\tau}\bar{s} = - \bar{s}\partial_i v^i $$ and $$ \frac{\d}{\d\tau} Q_i = -\frac{\partial_i P}{\gamma s} $$ where $\bar{s}=\tau\gamma s$ and $Q_i = -\frac{g_{ij}\omega u^j}{s}=-\frac{g_{ij}\tau \omega \gamma u^j}{\bar{s}}$.

The first step is to calculate the velocity by inverting the non-linear equation $$ g_{ij} Q^i Q^j = -\frac{\omega^2 }{s^2} g_{ij} u^i u^j = f(\frac{\bar{s}}{\tau\gamma})^2 (\gamma^2 - 1) $$ From this we have $\gamma$ and we can calculate $$ s = \frac{\bar{s}}{\tau \gamma} $$ and all other thermodynamical quantities, plus $$ u^i = -\frac{s}{\omega}g^{ij}Q_j $$

After this intermediary step we must compute the derivatives $\partial_iv^i$ and $\partial_i P$. For this we use the TaylorSPH algorithm. First we compute the following quantites $$ \rho_a = \sum_b W_{ab} \\ \partial_i \rho_a = \sum_b \partial_i W_{ab} \\ d_a^i = -\sum_b (r_a - r_b)^i W_{ab} \\ \partial_j d_a^i = -\sum_b (r_a - r_b)^i \partial_j W_{ab} \\ \ave{P}_a = \sum_b P_b W_{ab} \\ \ave{\partial_i P}_a = \sum_b P_b \partial_i W_{ab} \\ \ave{v^i}_a = \sum_b v^i_b W_{ab} \\ \ave{\partial_j v^i}_a = \sum_b v^i_b \partial_j W_{ab} $$ With these quantities calculated we can find the gradients by solving the linear system $$ M = [\begin{matrix} m_{11} & m_i \\ m'_j & m_{ij} \end{matrix}] = [\begin{matrix} \rho_a & \partial_i \rho_a \\ d_a^i & \partial_j d_a^i \end{matrix}]^{-1} $$

So that $$ \partial_i P = m'_i P + m_{ij} \ave{\partial_j P}_a \partial_i v^k = m'_i v^k + m_{ij} \ave{\partial_j v^k}_a $$

\section{Evolution} Given $r,s*,\frac{q}{s*}$ we calculate the derivatives as follows

We need to calculate



\section{Initial conditions} The initial conditions for the time integration need to be given in the form of the initial position $r_a$ of each test particle and its canonical momentum $p_a$. Apart from these one also needs a thermodynamic quantity, in our case we use the entropy density $s_a$.

\section{Setting up} Once these are determined, we compute the quantities needed to

\section{Fuse particles}

Here we introduce SPH particle fusing, which is a process to take place when there are too many SPH particles in a small and smooth region. Therefore it is advantageous to have these particles fused, in order to increase the efficincy in numerical calculation as well as to avoid unnecessary noise in SPH particle distribution ({\color{red} ?}).

We assume that most intensive quantities (such as the entropy density $s$, energy density $\omega$({\color{red}?}) and momentum density ${q^i}$) are taken to be the average over the particles to be fused. $$ s_{b'} = \frac{1}{N}\sum_b s_b \\ (q^i)_{b'} = (\omega u^i)_{b'} = \frac{1}{N}\sum_b (\omega u^i)_b \, $$ where the index $b$ only runs over those particles to be fused, $b'$ indicates the index of a new SPH particles resulting from the fusing. Another assumption is that that we want the total entropy and total momentum to be additive so that $$ S_{b'} = \sum_b \frac{\nu_b s_b}{\sigma_b} = \frac{\nu_{b'} s_{b'}}{\sigma_{b'}} \\ Q^i_{b'} = \sum_b \frac{\nu_b (\omega u^i)_{b}}{\sigma_b} = \frac{\nu_{b'} (\omega u^i)_{b'}}{\sigma_{b'}} $$ where $s_{b}$ and $\nu_{b}$ are the entropy density and conserved charge of the particles to be fused, the index $b$ only runs over those particles. Therefore one has $$ \frac{\nu_{b'}}{\sigma_{b'}} = \frac{S_{b'}}{s_{b'}} \\ \frac{\nu_{b'} (u^i)_{b'}}{\sigma_{b'}} = \frac{Q^i_{b'}}{\omega_{b'}} $$ where $\omega_{b'}$ is completely determined by $s_{b'}$ through the equation of state (EoS). The above expressions lead to $$ (u^i)_{b'} = \frac{s_{b'}}{S_{b'}} \frac{Q^i_{b'}}{\omega_{b'}} $$

One more step is needed to determine the corresponding kernel function for $b'$. First we express the reference density of the old particles $$ \rho_a = (\tau\gamma\sigma)_a = \sum_b \nu_b W_{ab} = \sum_b \nu_b W_{ab} + \sum_{c} \nu_c W_{ac} $$ where the index $b$ is for the particles which are fused and the index $c$ is for the particles which are not part of the fusing. We do not want to change the reference density of the particles that will not fuse, so the contribution from a particle which is not part of the fusing procedure remains unchanged $$ \rho_a = \sum_b \nu_b W_{ab} + \sum_{c} \nu_c W_{ac} = \nu_{b'} W_{ab'} + \sum_{c} \nu_c W_{ac} \\ \nu_{b'} = \sum_b \frac{ \nu_b W_{ab} }{ W_{ab'} } $$ the new kernel function is centered at $b'$ $$ W_{ab'} \equiv W(x_a - x_{b'}) $$ with the average $x_{b'} = \frac{(\sum_b \nu_b x_b)}{(\sum_b \nu_b)}$.

The properties of the new particle $b'$ are $$ \nu_{b'} = \sum_b \frac{ \nu_b W_{ab} }{ W_{ab'} }\\ s_{b'} = \frac{1}{N}\sum_b s_b \\ (q^i)_{b'} = \frac{1}{N}\sum_b (\omega u^i)_b \\ (u^i)_{b'} = \frac{s_{b'}}{S_{b'}} \frac{Q^i_{b'}}{\omega_{b'}} = \frac{s_{b'}}{\omega_{b'}}\frac{\sum_b \frac{\nu_b \omega_b u^i_b}{\sigma_b}}{\sum_b \frac{\nu_b s_b}{\sigma_b} } $$

\section{Smoothed Particle Hydrodynamics} These are some convenient quantites. We use the following notation we also use the notation $$ [f]_a = \sum_b f_b W_{ab} $$ and $$ [f]'_a = \sum_b f_b \partial W_{ab} $$ here we note that $b$ runs over all SPH particles.

The reference density $$ \rho_a = \sum_b \nu_b W_{ab} = [\nu]_a $$ its spatial derivative $$ \partial \rho_a = \sum_b \nu_b \partial W_{ab} = [\nu]'_a $$ The smoothed power $\lambda$ of the reference density $$ \Omega_a(\lambda) = \sum_b \nu_b \rho^\lambda_b W_{ab} = [\nu\rho^\lambda]_a = [\nu[\nu]^\lambda]_a $$ and its derivative $$ \partial \Omega_a(\lambda) = \sum_b \nu_b \rho^\lambda_b \partial W_{ab} = [\nu\rho^\lambda]'_a = [\nu[\nu]^\lambda]'_a $$

In what follows, we discuss two different types of implementation, each of which consists of different schemes. The first scenario, Type 1, is based on the original idea proposed by Monaghan $$ f_a = \rho_a^{-\lambda} \sum_b \nu_b f_b \rho_b^{\lambda-1} W_{ab} \\ = \Omega_a(0)^{-\lambda} \sum_b \nu_b f_b \rho_b^{\lambda-1} W_{ab}\\ = [\nu]_a^{-\lambda} [\nu f [\nu]^{\lambda-1}]_a $$ This means that $$ [\nu f [\nu]^{\lambda-1}]_a = [\nu]_a^{\lambda} f_a $$ where it reduces to the Monaghan's case for $\lambda=0$ $$ \lambda = 0\quad f_a = \sum_b \nu_b \frac{f_b}{\rho_b} W_{ab} = [\frac{\nu f}{\rho}]_a = [\frac{\nu f}{[\nu]}]_a $$ and also reduces to another proposed scheme for $\lambda=1$. This latter scheme ensures that linear functions are parametrized correctly by construction $$ \lambda = 1\quad f_a = \frac{1}{\rho_a} \sum_b \nu_b f_b W_{ab} = \frac{[\nu f]_a}{[\nu]_a} $$

Now we take the derivatives of these quantities to find that $$ \partial f_a = -\lambda \rho_a^{-\lambda-1}\partial \rho_a \sum_b \nu_b f_b \rho_b^{\lambda-1} W_{ab} + \rho_a^{-\lambda} \sum_b \nu_b f_b \rho_b^{\lambda-1} \partial W_{ab} \\ = -\lambda \rho_a^{-\lambda-1}\sum_c \nu_c \partial W_{ac} \sum_b \nu_b f_b \rho_b^{\lambda-1} W_{ab} + \rho_a^{-\lambda} \sum_c \nu_c f_c \rho_c^{\lambda-1} \partial W_{ac} \\ = \rho_a^{-\lambda-1} \sum_c \nu_c(\rho_a f_c \rho_c^{\lambda-1} -\lambda \sum_b \nu_b f_b \rho_b^{\lambda-1} W_{ab}) \partial W_{ac} $$ or in the notation $$ \partial f_a = [\nu]_a^{-\lambda} [\nu f [\nu]^{\lambda-1}]'_a - \lambda [\nu]_a^{-\lambda-1} [\nu]'_a [\nu f [\nu]^{\lambda-1}]_a $$ and using the relation above one gets $$ \partial f_a \equiv \sum_b F_{ab} = [\nu]_a^{-\lambda} [\nu f [\nu]^{\lambda-1}]'_a - \lambda [\nu]_a^{-1} [\nu]'_a f_a $$ $$ F_{ab} = \nu_b [\nu]_a^{-\lambda} ( f_b [\nu]^{\lambda-1}_b - \lambda [\nu]_a^{\lambda-1} f_a ) \partial W_{ab} $$ This reduces to the familiar cases for $\lambda=0$ $$ \lambda=0\quad \partial f_a = [\frac{\nu f}{[\nu]}]'_a $$ the sum is performed over $$ F_{ab} = \frac{\nu_b f_b}{[\nu]_b}\partial W_{ab} $$ which is the standard scheme for computing derivatives. $$ \lambda = 1\quad \partial f_a = \frac{1}{\rho_a^2} \sum_c \nu_c(\rho_a f_c - \sum_b \nu_b f_b W_{ab}) \partial W_{ac} \\ = [\nu]_a^{-1} [\nu f]'_a - [\nu]_a^{-2} [\nu]'_a [\nu f]_a \\ = [\nu]_a^{-2} ([\nu]_a [\nu f]'_a - [\nu]'_a [\nu f]_a) $$ or $$ \lambda = 1\quad \partial f_a = \frac{1}{[\nu]_a} ( [\nu f]'_a - [\nu]'_a f_a ) $$ where the sum is perfomed over $$ F_{ab} = \frac{\nu_b (f_b - f_a)}{[\nu]_a}\partial W_{ab} $$

For $\lambda=-1$ we have $$ \partial f_a = [\nu]_a [\nu f [\nu]^{-2}]'_a + [\nu]'_a [\nu f [\nu]^{-2}]_a \\ = [\nu]_a ( [ \frac{\nu f}{[\nu]^{2}}]'_a + \frac{[\nu]'_a f_a}{[\nu]_a^{2}} ) $$

Type 1a $$ \partial f_a = \rho_a^{-\lambda} \sum_c \nu_c(f_c \rho_c^{\lambda-1} -\lambda \rho_a^{\lambda-1} f_a) \partial W_{ac} $$ $$ \lambda = 1\quad \partial f_a = \frac{1}{\rho_a} \sum_c \nu_c( f_c - f_a ) \partial W_{ac} \\ \lambda = -1\quad \partial f_a = \rho_a \sum_c \nu_c( \frac{f_c}{\rho_c^2} + \frac{f_a}{\rho_a^2} ) \partial W_{ac} $$

One can also write $$ \partial f_a = \partial \frac{[\nu]^\alpha_a f_a}{[\nu]^\alpha_a} = [\nu]^\alpha_a \partial \frac{f_a}{[\nu]^\alpha_a} + \alpha \frac{f_a}{[\nu]_a} [\nu]'_a $$ By substituting $f_a \right f_a [\nu]_a^{-\lambda}$ into previous expression for $\partial f_a$ one has $$ \partial \frac{f_a}{[\nu]^\alpha_a} = [\nu]_a^{-\lambda} [\nu f [\nu]^{\lambda-\alpha-1}]'_a - \lambda [\nu]_a^{-\alpha-1} [\nu_a]' f_a $$ therefore $$ \partial f_a = [\nu]^\alpha_a ( [\nu]_a^{-\lambda} [\nu f [\nu]^{\lambda-\alpha-1}]'_a - \lambda [\nu]_a^{-\alpha-1} [\nu_a]' f_a ) + \alpha \frac{f_a}{[\nu]} [\nu]'_a \\ = [\nu]_a^{-\lambda+\alpha} [\nu f [\nu]^{\lambda-\alpha-1}]'_a - \lambda [\nu]_a^{-1} [\nu]'_a f_a + \alpha \frac{f_a}{[\nu]} [\nu]'_a \\ = [\nu]_a^{-\lambda+\alpha} [\nu f [\nu]^{\lambda-\alpha-1}]'_a - (\lambda-\alpha) [\nu]'_a \frac{f_a}{[\nu]_a} $$

so if we rewrite $\beta=\lambda-\alpha$ we have $$ \partial f_a = [\nu]_a^{-\beta} [\nu f [\nu]^{\beta-1}]'_a - \beta [\nu_a]' \frac{f_a}{[\nu]_a} $$ The resulting expression has the same form as the previous one, where $\lambda$ is replaced by $\beta$. ({\color{red} However, in derivation this expression, we have explicitly made use of the previous expression for $\partial f_a$, moreover, $\lambda$ is a fixed number for a given ``scheme, therefore the above derivation shows how two different schemes, namely $\lambda$ and $\beta \equiv \lambda-\alpha$ can be related. In other words, the current type is well defined up to a transformation connecting differnt values of $\lambda$.})

By summing SPH particles on both sides of the Taylor expansion of a function, one obtains the following expressions $$ f(x) = f_a + (x-x_a)\partial f_a \\ \int f(x) W(x-x_a) = f_a \int W(x-x_a) + \partial f_a \int (x-x_a) W(x-x_a) \\ \sum_b \nu_b f_b/\rho_b W_{ab} = f_a \sum_b \nu_b/\rho_b W_{ab} + \partial f_a \sum_b \nu_b/\rho_b (x_b-x_a)W_{ab} \\ f_a = \frac{\sum_b \nu_b f_b/\rho_b W_{ab}}{\sum_b \nu_b/\rho_b W_{ab}} \\ \sum_b \nu_b f_b/\rho_b \partial W_{ab} = f_a \sum_b \nu_b/\rho_b \partial W_{ab} + \partial f_a \sum_b \nu_b/\rho_b (x_b-x_a)\partial W_{ab} \\ \partial f_a = \frac{\sum_b \nu_b (f_b - f_a)/\rho_b \partial W_{ab}}{ \sum_b \nu_b/\rho_b (x_b-x_a)\partial W_{ab} } $$

Therefore, another family of interpolation schemes, Type 2, can be written as $$ f_a = \frac{\sum_b \nu_b f_b \rho_b^{-\lambda} W_{ab}}{\sum_b \nu_b \rho_b^{-\lambda} W_{ab}} = \frac{1}{\Omega_a(-\lambda)}\sum_b \nu_b f_b \rho_b^{-\lambda} W_{ab} = \frac{ [\nu f [\nu]^{-\lambda}]_a }{ [\nu [\nu]^{-\lambda}]_a } $$ and also $$ [\nu f [\nu]^{-\lambda}]_a = [\nu [\nu]^{-\lambda}]_a f_a $$ which reduces to $$ \lambda = 0\quad \frac{\sum_b \nu_b f_b W_{ab}}{\sum_b \nu_b W_{ab}} = \frac{[\nu f]_a}{[\nu]_a} $$ Now we take the derivatives $$ \partial f_a = \frac{\sum_b \nu_b f_b \rho_b^{-\lambda} \partial W_{ab}}{\sum_b \nu_b \rho_b^{-\lambda} W_{ab}} - \frac{\sum_b \nu_b f_b \rho_b^{-\lambda} W_{ab}}{( \sum_b \nu_b \rho_b^{-\lambda} W_{ab})^2} \sum_b \nu_b \rho_b^{-\lambda} \partial W_{ab} \\ = \frac{1}{( \sum_b \nu_b \rho_b^{-\lambda} W_{ab})^2}( \sum_b \nu_b \rho_b^{-\lambda} W_{ab} \sum_b \nu_b f_b \rho_b^{-\lambda} \partial W_{ab} - \sum_b \nu_b f_b \rho_b^{-\lambda} W_{ab} \sum_b \nu_b \rho_b^{-\lambda} \partial W_{ab} ) \\ = \frac{1}{\Omega_a(-\lambda)^2}( \Omega_a(-\lambda) \sum_b \nu_b f_b \rho_b^{-\lambda} \partial W_{ab} - \sum_b \nu_b f_b \rho_b^{-\lambda} W_{ab} \partial \Omega_a(-\lambda) ) \\ =\frac{1}{[\nu [\nu]^{-\lambda} ]_a^2}( [\nu[\nu]^{-\lambda}]_a[\nu f [\nu]^{-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a[\nu f[\nu]^{-\lambda}]_a ) \\ =\frac{1}{[\nu [\nu]^{-\lambda} ]_a}( [\nu f [\nu]^{-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a f_a ) $$

Now if we use $$ \partial f_a = [\nu]^\alpha_a \partial \frac{f_a}{[\nu]^\alpha_a} + \alpha \frac{f_a}{[\nu]_a} [\nu]'_a $$ we have $$ \partial \frac{f_a}{[\nu]^\alpha_a} = \frac{1}{[\nu [\nu]^{-\lambda} ]_a}( [\nu f [\nu]^{-\alpha-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a [\nu]_a^{-\alpha} f_a ) $$ so putting it all together $$ \partial f_a = [\nu]^\alpha_a (\frac{1}{[\nu [\nu]^{-\lambda} ]_a}( [\nu f [\nu]^{-\alpha-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a [\nu]_a^{-\alpha} f_a )) + \alpha \frac{f_a}{[\nu]_a} [\nu]'_a \\ = \frac{1}{[\nu [\nu]^{-\lambda} ]_a}( [\nu]^\alpha_a [\nu f [\nu]^{-\alpha-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a f_a ) + \alpha \frac{f_a}{[\nu]_a} [\nu]'_a $$

or one substitute $\nu \rightarrow \nu [\nu]^{-\beta}_\alpha$ in the above arguments $$ \partial f_a = [\nu[\nu]^{-\beta}]^\alpha_a \partial \frac{f_a}{[\nu[\nu]^{-\beta}]^\alpha_a} + \alpha \frac{f_a}{[\nu[\nu]^{-\beta}]_a} [\nu[\nu]^{-\beta}]'_a $$ $$ \partial \frac{f_a}{[\nu[\nu]^{-\beta}]^\alpha_a} = \frac{1}{[\nu [\nu]^{-\lambda} ]_a}( [\nu f [\nu[\nu]^{-\beta}]^{-\alpha}[\nu]^{-\lambda}]'_a - [\nu[\nu]^{-\lambda}]'_a [\nu[\nu]^{-\beta}]_a^{-\alpha} f_a ) $$ $$ \partial f_a = \frac{[\nu[\nu]^{-\beta}]^\alpha_a}{[\nu [\nu]^{-\lambda} ]_a}( 

 [\nu f [\nu[\nu]^{-\beta}]^{-\alpha}[\nu]^{-\lambda}]'_a  - [\nu[\nu]^{-\lambda}]'_a [\nu[\nu]^{-\beta}]_a^{-\alpha} f_a ) + \alpha \frac{f_a}{[\nu[\nu]^{-\beta}]_a} [\nu[\nu]^{-\beta}]'_a

\\ = \frac{[\nu[\nu]^{-\beta}]^\alpha_a}{[\nu [\nu]^{-\lambda} ]_a} [\nu f [\nu[\nu]^{-\beta}]^{-\alpha}[\nu]^{-\lambda}]'_a 

- [\nu[\nu]^{-\lambda}]'_a \frac{1}{[\nu [\nu]^{-\lambda} ]_a} f_a  + \alpha \frac{f_a}{[\nu[\nu]^{-\beta}]_a} [\nu[\nu]^{-\beta}]'_a

$$

if $\beta=\lambda$ $$ \partial f_a = [\nu[\nu]^{-\lambda}]^{\alpha-1}_a [\nu f [\nu[\nu]^{-\lambda}]^{-\alpha}[\nu]^{-\lambda}]'_a + (\alpha -1) \frac{f_a}{[\nu[\nu]^{-\lambda}]_a} [\nu[\nu]^{-\lambda}]'_a $$ ({\color{red} which ends up in the middle of nowhere. Many doubts at this point. For instance, the direvative $\partial f_a$ evaluated directly from $f_a$ is in fact, different from that at the beginning of the discussions of ``Type 2, where one obtains an expression by using $\sum_b\partial W_{ab} on both sides of the Taylor expansion$. So how to deal with so many different formulae.})


\section{Lagrangian coordinates} $$ f(x, t) = \int dx' f(x',t) \delta(x-x') $$ $$ f(r(t), t) = \int dr(t) \abs{\frac{dx}{dr(t)}} f(r(t),t) \delta(x-r(t)) $$ $$ d f(r(t), t) = d_t r^i(t)\partial_i f + \partial_t f $$



\section{Kernel}

({\color{red} I am kinda lost in this section})

Given the kernel $W$ and its property that $$ W(x-r(t);h) = \frac{1}{h^d} \phi(\frac{\abs{x-r(t)}}{h}) $$ where $\phi(\frac{\abs{x-r(t)}}{h}) = 0$ for $\abs{x-r(t)} > 2h$. We can expand the kernel around $q = \frac{\abs{x-r(t)}}{2h} < 1$. $$ \phi(q) = \phi(0) + q\phi'(0) + \dots $$ Now we multiply by $1/h^d$ $$ W(x-r(t);h) = \frac{1}{h^d}\phi(q) = \frac{1}{h^d}\phi(0) + \frac{q}{h^d}\phi'(0) + \dots $$ or inversely, one may expand a $\delta$ function ($h=0$) around a Gaussian kernel ($h=h$) $$ h^d \delta(x-r(t)) = h^d W(x-r(t);h=0) \nq = h^d W(x-r(t);h) 

- h \frac{\delta h^d W(x-r(t);h)}{\delta h}  
+ \frac{h^{2}}{2} \frac{\delta^2 h^d W(x-r(t);h)}{\delta h^2}

$$ $$ \int dx W(x-r(t);h) = 1 $$ Kernels of the Gaussian-type have the general form of $$ h^d W(x-r(t);h) = \phi(\frac{\abs{x-r(t)}}{h}) $$ therefore the derivative is $$ \frac{\delta \phi(\frac{\abs{x-r(t)}}{h}) }{\delta h} = - \frac{ \abs{x-r(t)} }{h^{2}} \phi'(\frac{\abs{x-r(t)}}{h}) $$

%On the other hand %$$ %\frac{\delta W(x-r(t);h)}{ \delta r(t) } = \frac{ x-r(t) }{ \abs{x-r(t)} h^{d+1} } \phi'(\frac{\abs{x-r(t)}}{h}) %$$

If we put the two first terms of the Taylor expansion together we get $$ h^d \delta(x-r(t)) = \phi(\frac{\abs{x-r(t)}}{h}) + \frac{ \abs{x-r(t)} }{h^{2}} \phi'(\frac{\abs{x-r(t)}}{h}) $$ or $$ W_2(x-r(t);h) = \frac{1}{h^d}\phi(\frac{\abs{x-r(t)}}{h}) + \frac{ \abs{x-r(t)} }{h^{d+2}} \phi'(\frac{\abs{x-r(t)}}{h}) $$ where the index 2 indicates that it is accurated up to second order

Another attempt is to reduce the width from $h$ to $h'$. $$ W(x-r(t);h) = W(x-r(t);h') + (h-h') \frac{\delta W(x-r(t);h')}{\delta h} + \frac{(h-h')^{2}}{2} \frac{\delta^2 W(x-r(t);h')}{\delta h^2} $$

If we put the first two terms together $$ W(x-r(t);h') 

= \frac{1}{h^{d}} \phi(\frac{\abs{x-r(t)}}{h})
+ (h-h') ( \frac{d}{h^{d+1}} \phi(\frac{\abs{x-r(t)}}{h}) 
+ \frac{ \abs{x-r(t)} }{h^{d+2}} \phi'(\frac{\abs{x-r(t)}}{h}) )

$$

Kernel derivatives $$ \partial_i W(x-r_a;h) = \frac{ \partial}{\partial x_i}( \frac{\abs{x-r_a}}{h} ) W'(x-r_a;h) 

= \frac{(x-r_a)_i}{h\abs{x-r_a}} W'(x-r_a;h)

$$ where we used that $$ \partial_i \abs{x-r_a} = \partial_i ((x-r_a)_k(x-r_a)_k )^{1/2} \nq = ((x-r_a)_k(x-r_a)_k )^{-1/2} (x-r_a)_k \delta_{ik} \nq = \frac{(x-r_a)_i}{\abs{x-r_a}} $$ $$ \partial_j \partial_i W(x-r_a;h) = \partial_j( \frac{ (x-r_a)_i }{ h\abs{x-r_a} } W'(x-r_a;h) ) \nq = \partial_j( \frac{ (x-r_a)_i }{ h\abs{x-r_a} } ) W'(x-r_a;h) 

+ \frac{(x-r_a)_i}{h\abs{x-r_a}} \partial_j W'(x-r_a;h)

\nq = ( \frac{ \delta_{ij} }{h\abs{x-r_a}} - \frac{(x-r_a)_i(x-r_a)_j}{h\abs{x-r_a}^3} ) W'(x-r_a;h) 

+ \frac{(x-r_a)_i(x-r_a)_j}{h^2\abs{x-r_a}^2} W(x-r_a;h)

$$ So $$ \partial^2 W(x-r_a;h) \nq = ( \frac{ 1 }{h\abs{x-r_a}} - \frac{(x-r_a)^2}{h\abs{x-r_a}^3} ) W'(x-r_a;h) 

+ \frac{ (x-r_a)^2 }{h^2\abs{x-r_a}^2} W(x-r_a;h)

\nq = \frac{ 1 }{h^2} W(x-r_a;h) $$ and for $i \ne j$ $$ \partial_j \partial_i W(x-r_a;h) \nq = - \frac{(x-r_a)_i(x-r_a)_j}{h\abs{x-r_a}^3} W'(x-r_a;h) 

+ \frac{(x-r_a)_i(x-r_a)_j}{h^2\abs{x-r_a}^2} W(x-r_a;h)

$$



\section{Taylor SPH}

From the Taylor series expansion we have $$ f(x) = \sum_n \frac{ (x-r_a)^n_i }{n!}(\partial^n_i f)_{r_a} $$ Now we multiply both sides by $W(x-r_a)$ and integrate over $x$ $$ \int dx f(x) W(x-r_a) = \sum_n \frac{ (\partial^n_i f)_{r_a} }{n!} \int dx (x-r_a)^n_i W(x-r_a) $$ If we use only the first order in the expansion we have $$ \int dx f(x) W(x-r_a) = f_a \int dx W(x-r_a) $$ In principle, the integral over the kernel is equal to one. However, once we introduce the particle approximation, namely, $$ \int dx f(x) W(x-r_a) \rightarrow \sum_b \frac{\nu_b f_b}{\rho_b} W(r_b-r_a) \\ \int dx W(x-r_a) \rightarrow \sum_b \frac{\nu_b }{\rho_b} W(r_b-r_a) $$ it is noting that the second expression does not exactly goes to ``1 in this approximation, therefore $$ \sum_b \frac{\nu_b f_b}{\rho_b} W(r_b-r_a) = f_a \sum_b \frac{\nu_b}{\rho_b} W(r_b-r_a) $$ Finally we have $$ f_a = \frac{\sum_b \frac{\nu_b f_b}{\rho_b} W(r_b-r_a)}{\sum_b \frac{\nu_b}{\rho_b} W(r_b-r_a)} $$ this reduces to the standard SPH form when the sum over the kernel is exactly one.

The same can be done for the gradient term as follows $$ \int dx f(x) \frac{x-r_a}{\abs{x-r_a}}W'(x-r_a) = \sum_n \frac{ (\partial^n_i f)_{r_a} }{n!} \int dx (x-r_a)^n_i \frac{x-r_a}{\abs{x-r_a}} W'(x-r_a) $$ Now we also keep up to the second order term $$ \int dx f(x) \frac{(x-r_a)_i}{\abs{x-r_a}}W'(x-r_a) = f_a \int dx \frac{(x-r_a)_i}{\abs{x-r_a}} W'(x-r_a) + \partial_j f_a \int dx (x-r_a)_j \frac{(x-r_a)_i}{\abs{x-r_a}} W'(x-r_a) $$ Using the particle approximation we get $$ \sum_b \frac{\nu_b f_b}{\rho_b} \frac{(r_b-r_a)_i}{\abs{r_b-r_a}}W'(r_b-r_a) \nq = f_a \sum_b \frac{\nu_b}{\rho_b} \frac{(r_b-r_a)_i}{\abs{r_b-r_a}} W'(r_b-r_a) 

+ \partial_j f_a \sum_b \frac{\nu_b}{\rho_b} (r_b-r_a)_j \frac{(r_b-r_a)_i}{\abs{r_b-r_a}} W'(r_b-r_a)

$$ $$ \partial_j f_a \sum_b \frac{\nu_b}{\rho_b} \frac{(r_b-r_a)_j (r_b-r_a)_i}{\abs{r_b-r_a}} W'(r_b-r_a) = \sum_b \frac{\nu_b (f_b - f_a)}{\rho_b} \frac{(r_b-r_a)_i}{\abs{r_b-r_a}}W'(r_b-r_a) $$ In one dimension this becomes $$ \partial_x f_a 

=

\frac{\sum_b \frac{\nu_b (f_b - f_a)}{\rho_b} \frac{(r_b-r_a)}{\abs{r_b-r_a}}W'(r_b-r_a) }{ \sum_b \frac{\nu_b}{\rho_b} \abs{r_b-r_a} W'(r_b-r_a)} 

= \frac{\sum_b \frac{\nu_b (f_b - f_a)}{\rho_b} F_{ab} }{ \sum_b \frac{\nu_b}{\rho_b} r_{ba} F_{ab} }
= \frac{ \aave{f} }{ \aave{r} }

$$ In two dimensions $$ [\begin{matrix} \partial_x f_a & \partial_y f_a \end{matrix}] \sum_b \frac{\nu_b}{\rho_b\abs{r_{ba}}} [\begin{matrix} (r_{ba})_x^2 & (r_{ba})_x (r_{ba})_y\\ (r_{ba})_y (r_{ba})_x & (r_{ba})_y^2 \end{matrix}] W'(r_b-r_a) \nq = \sum_b \frac{\nu_b (f_b - f_a)}{\rho_b\abs{r_{ba}}} [\begin{matrix} (r_b-r_a)_x & (r_b-r_a)_y \end{matrix}] W'(r_b-r_a) $$ Schematically we get $$ \partial_i f_a = M_{ij} \sum_b \frac{\nu_b (f_b - f_a)}{\rho_b} \frac{ (r_{ba})_j }{ \abs{r_{ba}} } W'(r_b-r_a) \nq = M_{ij} \sum_b \frac{\nu_b f_{ba}}{\rho_b} (F_{ba})_j \nq = M_{ij} \aave{ f }_j $$ where $$ (F_{ba})_j = \frac{ (r_{ba})_j }{ \abs{r_{ba}} } W'_{ba} $$

$$ M_{ij} = [\begin{matrix} \sum_b \frac{\nu_b}{\rho_b\abs{r_{ba}}} (r_{ba})_x^2 W'(r_b-r_a) & \sum_b \frac{\nu_b}{\rho_b\abs{r_{ba}}} (r_{ba})_x(r_{ba})_y W'(r_b-r_a)\\ \sum_b \frac{\nu_b}{\rho_b\abs{r_{ba}}} (r_{ba})_x(r_{ba})_y W'(r_b-r_a) & \sum_b \frac{\nu_b}{\rho_b\abs{r_{ba}}} (r_{ba})_y^2 W'(r_b-r_a) \end{matrix}]^{-1} \nq = [\begin{matrix} \sum_b \frac{\nu_b}{\rho_b} (r_{ba})_x (F_{ba})_x & \sum_b \frac{\nu_b}{\rho_b} (r_{ba})_x (F_{ba})_y\\ \sum_b \frac{\nu_b}{\rho_b} (r_{ba})_y (F_{ba})_x & \sum_b \frac{\nu_b}{\rho_b} (r_{ba})_y (F_{ba})_y \end{matrix}]^{-1} \nq = [\begin{matrix} \aave{r_x}_x & \aave{r_x}_y\\ \aave{r_y}_x & \aave{r_y}_y\\ \end{matrix}]^{-1} \nq = \frac{1}{\aave{r_x}_x \aave{r_y}_y - \aave{r_x}_y \aave{r_y}_x } [\begin{matrix} \aave{r_y}_y & -\aave{r_x}_y\\ -\aave{r_y}_x & \aave{r_x}_x\\ \end{matrix}] $$

Now, one may ask why the above approach did not consider the second order expansion in the case of kernel integral as in that of the gradient of kernel. This leads to the following approach, which gives rise to a system of coupled equations. One may invert simultaneously for $f_a$ and $f_{a,i}$ for the solution. Using the notation $$ \ave{v_i}_a = \sum_b \frac{\nu_b (v_i)_b}{\rho_b} W_{ab} $$ and $$ \ave{v_i}_{a,j} = \sum_b \frac{\nu_b (v_i)_b}{\rho_b} \frac{(r_{ab})_j}{\abs{r_{ab}}}W_{ab} \nq = \sum_b \frac{\nu_b (v_i)_b}{\rho_b} (F_j)_{ab} $$ we write $$ \ave{v_i}_a = (v_i)_a \ave{1}_a + (v_i)_{a,j} \ave{\Delta r_j}_a \\ \ave{v_i}_{a,j} = (v_i)_a \ave{1}_{a,j} + (v_i)_{a,k} \ave{\Delta r_k}_{a,j} $$ which in matrix form is $$ [\begin{matrix} \ave{v_i}_a \\ \ave{v_i}_{a,j} \end{matrix}] = [\begin{matrix} 

\ave{1}_a & \ave{\Delta r_k}_a \\
\ave{1}_{a,j} & \ave{\Delta r_k}_{a,j}

\end{matrix}] [\begin{matrix} (v_i)_a \\ (v_i)_{a,k} \end{matrix}] $$

We can see that in the case of standart SPH, the terms $\ave{\Delta r_k}_a$ and $\ave{1}_{a,j}$ are taken to be exactly zero and the terms $\ave{1}_a$ and $\ave{\Delta r_k}_{a,j}$ to be exactly one. This is true for the continuous case, i.e. ({\color{red} $-\delta_{jk}$?}) $$ \int dx W(x-r_a) = 1 \\ \int dx (x-r_a)_k \partial_j W(x-r_a) = \delta_{jk} $$ and $$ \int dx \partial_j W(x-r_a) = 0 \\ \int dx (x-r_a)_k W(x-r_a) = 0 $$ However, this is not exactly satisfied in the particle approximation.

The previous approach of taylor SPH also supposes that $\ave{\Delta r_k}_a = 0$. In this formulation all contributions are included. It is important to notice that the system above corresponds to $d+1$ variables and equations where $d$ is the spatial dimension of the system. The SPH quantities are then calculated to be $$ [\begin{matrix} (v_i)_a \\ (v_i)_{a,k} \end{matrix}] = [\begin{matrix} 

\ave{1}_a & \ave{\Delta r_k}_a 
\\
\ave{1}_{a,j} & \ave{\Delta r_k}_{a,j}

\end{matrix}]^{-1} [\begin{matrix} \ave{v_i}_a \\ \ave{v_i}_{a,j} \end{matrix}] $$

Extension to 3rd order we write $$ \ave{v_i}_a = (v_i)_a \ave{1}_a + (v_i)_{a,j} \ave{\Delta r_j}_a + \frac{(v_i)_{a,jk} }{2} \ave{\Delta r_j \Delta r_k }_a \\ \ave{v_i}_{a,j} = (v_i)_a \ave{1}_{a,j} + (v_i)_{a,k} \ave{\Delta r_k}_{a,j} + \frac{(v_i)_{a,kl} }{2} \ave{\Delta r_k \Delta r_l }_{a,j} \\ \ave{v_i}_{a,jk} = (v_i)_a \ave{1}_{a,jk} + (v_i)_{a,l} \ave{\Delta r_l}_{a,jk} + \frac{(v_i)_{a,lm} }{2} \ave{\Delta r_l \Delta r_m }_{a,jk} $$ which in matrix form is $$ [\begin{matrix} \ave{v_i}_a \\ \ave{v_i}_{a,j} \\ \ave{v_i}_{a,jk} \end{matrix}] = [\begin{matrix} \ave{1}_a & \ave{\Delta r_l}_a & \frac{1}{2} \ave{\Delta r_l \Delta r_m}_a\\ \ave{1}_{a,j} & \ave{\Delta r_l}_{a,j} & \frac{1}{2}\ave{\Delta r_l \Delta r_m}_{a,j}\\ \ave{1}_{a,jk} & \ave{\Delta r_l}_{a,jk} & \frac{1}{2} \ave{\Delta r_l \Delta r_m}_{a,jk} \end{matrix}] [\begin{matrix} (v_i)_a \\ (v_i)_{a,l} \\ (v_i)_{a,lm} \end{matrix}] $$

Since in most pratical applications, one only need the first order derivative, therefore the second order formalism is that useful. Nevertheless, one can makes use of the difference between first and second order formalism as an estimation of accuracy.



\section{Comparison}

In fig.\ref{fig:compare} we compare several schemes for both the function interpolation and its derivative. Below we describe the schemes plotted. We use the notation $$ \ave{v_i}_a = \sum_b \frac{\nu_b (v_i)_b}{\rho_b} W_{ab} \\ \ave{v_i}_{a,j} = \sum_b \frac{\nu_b (v_i)_b}{\rho_b} \frac{(r_{ab})_j}{\abs{r_{ab}}}W_{ab} \\ \ave{\Delta v_i}_a = \sum_b \frac{\nu_b ((v_i)_b-(v_i)_a) }{\rho_b} W_{ab} = \ave{v_i}_a - (v_i)_a \ave{1}_a $$

Some other convenient $$ \ave{v_i}_{a|b} = \frac{\nu_b (v_i)_b}{\rho_b} W_{ab} \\ \ave{v_i}_{a|b,j} = \frac{\nu_b (v_i)_b}{\rho_b} \frac{(r_{ab})_j}{\abs{r_{ab}}}W_{ab} \\ \ave{\Delta v_i}_{a|b} = \frac{\nu_b ((v_i)_b-(v_i)_a) }{\rho_b} W_{ab} = \ave{v_i}_{a|b} - (v_i)_a \ave{1}_{a|b} $$


\newpage \subsection{standardSPH} \begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{standardSPH_1f.png} \includegraphics[width=.4\columnwidth]{standardSPH_1df.png} \\ \includegraphics[width=.4\columnwidth]{standardSPH_0_125f.png} \includegraphics[width=.4\columnwidth]{standardSPH_0_125df.png} \caption{comparison of standardSPH interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare standardSPH} \end{figure}

The label standardSPH stands for $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 1 & 0\\ 0 & \mathbf{1} \end{matrix}]^{-1} [ \begin{matrix} \ave{f}_a\\ \ave{f}_{a,i} \end{matrix} ] = [ \begin{matrix} \ave{f}_a\\ \ave{f}_{a,i} \end{matrix} ] $$


\newpage \subsection{TaylorSPH} \begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{TaylorSPH_1f.png} \includegraphics[width=.4\columnwidth]{TaylorSPH_1df.png} \\ \includegraphics[width=.4\columnwidth]{TaylorSPH_0_125f.png} \includegraphics[width=.4\columnwidth]{TaylorSPH_0_125df.png} \caption{comparison of TaylorSPH interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare TaylorSPH} \end{figure}

\begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{TaylorSPH2nd_1f.png} \includegraphics[width=.4\columnwidth]{TaylorSPH2nd_1df.png} \\ \includegraphics[width=.4\columnwidth]{TaylorSPH2nd_0_125f.png} \includegraphics[width=.4\columnwidth]{TaylorSPH2nd_0_125df.png} \caption{comparison of TaylorSPH2nd interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare TaylorSPH2nd} \end{figure}


The label TaylorSPH uses $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 

\ave{1}_a & \ave{\Delta r_i}_a \\
\ave{1}_{a,j} & \ave{\Delta r_i}_{a,j}

\end{matrix}]^{-1} [ \begin{matrix} \ave{f}_a\\ \ave{f}_{a,j} \end{matrix} ] $$ which in 1d is expressed explicitly as $$ f_a = \frac{ \ave{\Delta r}_{a,x} \ave{f}_a - \ave{\Delta r}_{a} \ave{f}_{a,x} }{ \ave{1}_a\ave{\Delta r}_{a,x} - \ave{1}_{a,x}\ave{\Delta r}_{a} } \\ f_{a,x} = \frac{ \ave{1}_{a} \ave{f}_{a,x} - \ave{1}_{a,x} \ave{f}_{a} }{ \ave{1}_a\ave{\Delta r}_{a,x} - \ave{1}_{a,x}\ave{\Delta r}_{a} } $$ or $$ f_a = \frac{ \ave{\Delta r \ave{f}_a - \ave{\Delta r}_{a} f }_{a,x} }{ \ave{ \Delta r \ave{1}_a - \ave{\Delta r}_{a} }_{a,x} } \\ f_{a,x} = \frac{ \ave{ \ave{1}_{a} f - \ave{f}_{a} }_{a,x} }{ \ave{ \ave{1}_a \Delta r - \ave{\Delta r}_{a}}_{a,x} } $$


\newpage \subsection{truncTaylorSPH} \begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{truncTaylorSPH_1f.png} \includegraphics[width=.4\columnwidth]{truncTaylorSPH_1df.png} \\ \includegraphics[width=.4\columnwidth]{truncTaylorSPH_0_125f.png} \includegraphics[width=.4\columnwidth]{truncTaylorSPH_0_125df.png} \caption{comparison of truncTaylorSPH interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare truncTaylorSPH} \end{figure}

The label truncTaylorSPH $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 

\ave{1}_a & 0 \\
\ave{1}_{a,i} & \ave{\Delta r_j}_{a,i}

\end{matrix}]^{-1} [ \begin{matrix} \ave{f}_a\\ \ave{f}_{a,j} \end{matrix} ] $$ or, in 1d, $$ f_a = \frac{\ave{f}_a}{\ave{1}_a} \\ f_{a,i} = \frac{\ave{f}_{a,i} - f_a \ave{1}_{a,i}}{ \ave{\Delta r_i}_{a,i} } $$

The label invrho1 is expressed as $$ f_a = \frac{\ave{f}_a}{\ave{1}_a} \\ f_{a,i} = \frac{\ave{f}_{a,i} - f_a \ave{1}_{a,i} }{\ave{1}_a } $$ or $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 

\ave{1}_a & 0 \\
\ave{1}_{a,i} & \ave{1}_a

\end{matrix}]^{-1} [ \begin{matrix} \ave{f}_a\\ \ave{f}_{a,j} \end{matrix} ] $$


\subsection{invrho-1} The label invrho-1 is expressed as $$ f_a = \frac{\ave{\rho^2 f}_a}{\ave{\rho^2}_a} \\ f_{a,i} = \frac{\ave{\rho^2 f}_{a,i} - f_a \ave{\rho^2}_{a,i}}{ \ave{\rho^2}_a } $$ $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 

\ave{\rho^2}_a & 0 \\
\ave{\rho^2}_{a,i} & \ave{\rho^2}_a

\end{matrix}]^{-1} [ \begin{matrix} \ave{\rho^2 f}_a\\ \ave{\rho^2 f}_{a,j} \end{matrix} ] $$

\newpage \subsection{rho1} \begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{rho1_1f.png} \includegraphics[width=.4\columnwidth]{rho1_1df.png} \\ \includegraphics[width=.4\columnwidth]{rho1_0_125f.png} \includegraphics[width=.4\columnwidth]{rho1_0_125df.png} \caption{comparison of rho1 interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare rho1} \end{figure}

The label rho1 is expressed as $$ f_a = \frac{\ave{\rho f}_a}{\rho_a} \\ f_{a,i} = \frac{\ave{\rho f}_{a,i} - f_a \rho_{a,i}}{ \rho_a } $$ $$ [ \begin{matrix} f_a\\ f_{a,i} \end{matrix} ] = [\begin{matrix} 

\rho_a & 0 \\
\rho_{a,i} & \rho_a

\end{matrix}]^{-1} [ \begin{matrix} \ave{\rho f}_a\\ \ave{\rho f}_{a,j} \end{matrix} ] $$


\newpage \subsection{rho-1} \begin{figure}[!h] \centering \includegraphics[width=.4\columnwidth]{rho-1_1f.png} \includegraphics[width=.4\columnwidth]{rho-1_1df.png} \\ \includegraphics[width=.4\columnwidth]{rho-1_0_125f.png} \includegraphics[width=.4\columnwidth]{rho-1_0_125df.png} \caption{comparison of rho-1 interpolation performance for $h$-ratio 1(top) and 1/8(bottom). The function is shown on the left and the derivative on the right.} \label{fig:compare rho-1} \end{figure}

The label rho-1 is expressed as $$ f_a = \ave{\frac{f}{\rho}}_a \rho_a \\ f_{a,i} = \rho_a( \ave{\frac{f}{\rho}}_{a,i} + \frac{f_a}{\rho_a^2} \rho_{a,i} ) $$ $$ [ \begin{matrix} \frac{f_a}{\rho_a}\\ \frac{f_{a,i}}{\rho_a} \end{matrix} ] = [\begin{matrix} 

1 & 0 \\
-\frac{\rho_{a,i}}{\rho_a} & 1

\end{matrix}]^{-1} [ \begin{matrix} \ave{\frac{f}{\rho}}_a\\ \ave{\frac{f}{\rho}}_{a,j} \end{matrix} ] $$

This last one has the very important feature that $$ \sum_a \Omega_a f_{a,i} = \sum_a \Omega_a \rho_a( \ave{\frac{f}{\rho}}_{a,i} + \frac{f_a}{\rho_a^2} \rho_{a,i} ) \nq = \sum_{a,b} \nu_b \Omega_a \rho_a( \frac{f_b}{\rho_b^2} + \frac{f_a}{\rho_a^2} ) F_{ab,i} $$ and since $F_{ab,i} = - F_{ba,i}$ we conclude that if $\Omega_a = \frac{\nu_a}{\rho_a}$ $$ \sum_a \frac{\nu_a f_{a,i}}{\rho_a} = \sum_{a,b} \nu_a \nu_b( \frac{f_b}{\rho_b^2} + \frac{f_a}{\rho_a^2} ) F_{ab,i} = 0 $$ we therefore identify $$ \sum_a \frac{\nu_a f_{a,i}}{\rho_a} = \int dx \partial_i f = 0 $$ This allows us to express the integral in terms of the sum as $$ \int dx f = \sum_a \frac{\nu_a f_{a}}{\rho_a} $$ now using the first expression we have $$ \sum_a \frac{\nu_a f_{a}}{\rho_a} = \sum_{ab} \frac{\nu_a\nu_b f_{b}}{\rho_b^2} W_{ab} \nq = \sum_{b} \frac{\nu_b f_{b}}{\rho_b^2} \sum_a \nu_a W_{ab} \nq = \sum_{b} \frac{\nu_b f_{b}}{\rho_b} $$ which means that the integral is consistent with the first formula.

However, this is very inacurate, as seen in the figure. On of the reasons is that if we set $f=1$ we get, $$ f_a = \ave{\frac{1}{\rho}}_a \rho_a \neq 1 \\ f_{a,i} = \rho_a( \ave{\frac{1}{\rho}}_{a,i} + \frac{1}{\rho_a^2} \rho_{a,i} ) \neq 0 \\ \int dx f = \sum_a \frac{\nu_a}{\rho_a} $$ Furthermore, by setting $f=\rho$ we have $$ f_a = \ave{1}_a \rho_a \neq \rho_a \\ f_{a,i} = \rho_a\ave{1}_{a,i} + \rho_{a,i} \neq \rho_{a,i} \\ \int dx f = \sum_a \nu_a $$ In fact, for $f=\rho^2$ we have $$ f_a = \ave{\rho}_a \rho_a = \rho_a^2 \\ f_{a,i} = 2\rho_a \rho_{a,i} = {\rho_a^2}_{,i} \\ \int dx f = \sum_a \nu_a \rho_a $$

\begin{figure} \centering \includegraphics[width=\columnwidth]{compare.png} \caption{Comparison among several interpolation schemes described in the text. } \label{fig:compare} \end{figure}

The $y$ axis of the plot represents the relative displacement of the interpolated function and the actual function multiplied by the computanional time needed to perform the calculation $$ \Delta t \Delta f = \Delta t \sqrt{\frac{\sum_a (f^{`SPH}_a - f_a)^2}{\sum_a f_a^2}} $$ The solid lines are the interpolation of the function itself based on the weights and the values of the function in the particle points and the dashed line represents the estimate of the derivative also based only on the weights and the value of the function itself in the particle points. The analytical value of the derivative of the function is not used for the derivative estimation.

The TaylorSPH estimate of $f$ is the only that shows a rate of precision enhancing larger than the computational cost. It means that the time per precision is atually reduced for higher precision, $\frac{{\Delta f}_>}{{\Delta f}_<} < \frac{\Delta t_>}{\Delta t_<} $ The truncTaylorSPH, rho1 and invrho1 estimates for $f$ are very similar in precision. However, they reach a constant time per precision and do not improve for better precision, i.e., $\frac{{\Delta f}_>}{{\Delta f}_<} \sim \frac{\Delta t_>}{\Delta t_<} $.

The invrho-1 is a bit less precise but also resembles the previous cases. On the other hand the standarSPH and the rho-1 estimates of $f$ are inefficient in the sense that they take more time per precision the more precise it gets, the more expensive it becomes, namely, $\frac{{\Delta f}_>}{{\Delta f}_<} > \frac{\Delta t_>}{\Delta t_<} $.


\section{General expression}

We start with the Taylor expansion $$ f(x) = f_a + (x-r_a)_i f_{a,i} $$ Now we apply a general operation $$ \ave{f}_c \equiv \sum_b \phi_b f_b W_{cb} $$ where $\ave{f}$ is used to denote whatsoever operation in general, which essentially involves an approximation of the function $f$ at point $r_c$ by the summation of discrete SPH degrees of freedom at points $r_b$. ({\color{red} In the particular case discussed below}) $W_{ab}$ is the kernel fuction and $\phi_b$ is some multiple such as a weight in the average. Apply this to both sides of the Taylor expansion and we obtain $$ \ave{f}_a = \ave{1}_a f_a + \ave{(x-r_a)_i}_a f_{a,i} $$ Now we take the derivative with respect to $r_a$ $$ -\ave{f}_{a,j} = -\ave{1}_{a,j} f_a + \ave{1}_{a} f_{a,j} -\ave{(x-r_a)_i}_{a,j} f_{a,i} + \ave{(x-r_a)_i}_{a} f_{a,ij} $$ which natural leads to an infinite chain of equations. For the moment let us cut them by only considering terms up to the first-order derivative $$ [\begin{matrix} \ave{f}_a\\ \ave{f}_{a,i} \end{matrix}] = [\begin{matrix} \ave{1}_a & \ave{(x-r_a)_i}_a\\ \ave{1}_{a,j} & -\ave{1}_{a}\delta_{ij} + \ave{(x-r_a)_i}_{a,j} \end{matrix}] [\begin{matrix} f_a\\ f_{a,j} \end{matrix}] $$ At this point, one may compare with those of Taylor SPH obtained bofore and notice there is an extra term. The derivative in respect to $r_b$ is $$ \ave{f}_{a,j} = \ave{1}_{a,j} f_a + \ave{(x-r_a)_i}_{a,j} f_{a,i} $$ $$ [\begin{matrix} \ave{f}_a\\ \ave{f}_{a,j} \end{matrix}] = [\begin{matrix} \ave{1}_a & \ave{(x-r_a)_i}_a\\ \ave{1}_{a,j} & \ave{(x-r_a)_i}_{a,j} \end{matrix}] [\begin{matrix} f_a\\ f_{a,j} \end{matrix}] $$ In 1d it takes the following form $$ f_a = \frac{ \ave{x - r}_{a,x} \ave{f}_a - \ave{x - r}_{a} \ave{f}_{a,x} }{ \ave{1}_a\ave{x - r}_{a,x} - \ave{1}_{a,x}\ave{x - r}_{a} } \\ f_{a,x} = \frac{ \ave{1}_{a} \ave{f}_{a,x} - \ave{1}_{a,x} \ave{f}_{a} }{ \ave{1}_a\ave{x - r}_{a,x} - \ave{1}_{a,x}\ave{x - r}_{a} } $$ which can be condensed into $$ f_a = \frac{(\frac{\ave{f}}{\ave{x-r}})_{a,x}}{(\frac{\ave{1}}{\ave{x-r}})_{a,x}} 

= (\frac{\ave{f}_a}{\ave{1}_a})^2 \frac{(\frac{\ave{x-r}}{\ave{f}})_{a,x}}{(\frac{\ave{x-r}}{\ave{1}})_{a,x}}

\\ f_{a,x} = \frac{(\frac{\ave{f}}{\ave{1}})_{a,x}}{(\frac{\ave{x-r}}{\ave{1}})_{a,x}} 

= (\frac{\ave{f}_a}{\ave{x-r}_a})^2 \frac{(\frac{\ave{1}}{\ave{f}})_{a,x}}{(\frac{\ave{1}}{\ave{x-r}})_{a,x}}

$$

$$ f_a = \frac{ \ave{x}_{a,x}\ave{f}_a - r_a\ave{1}_{a,x}\ave{f}_a - \ave{x}_a \ave{f}_{a,x} + r_a \ave{f}_{a,x} } { \ave{1}_a\ave{x}_{a,x} - r_a\ave{1}_a\ave{1}_{a,x} - \ave{1}_{a,x}\ave{x}_{a} + r_a\ave{1}_{a}\ave{1}_{a,x} } \\ f_{a,x} = \frac{ \ave{1}_{a} \ave{f}_{a,x} - \ave{1}_{a,x} \ave{f}_{a} } { \ave{1}_a\ave{x}_{a,x} - r_a\ave{1}_a\ave{1}_{a,x} - \ave{1}_{a,x}\ave{x}_{a} + r_a\ave{1}_{a}\ave{1}_{a,x} } $$

$$ f_a = \frac{ \ave{x}_{a,x}\ave{f}_a - r_a\ave{1}_{a,x}\ave{f}_a - \ave{x}_a \ave{f}_{a,x} + r_a \ave{f}_{a,x} } { \ave{1}_a\ave{x}_{a,x} - \ave{1}_{a,x}\ave{x}_{a} } \\ f_{a,x} = \frac{ \ave{1}_{a} \ave{f}_{a,x} - \ave{1}_{a,x} \ave{f}_{a} } { \ave{1}_a\ave{x}_{a,x} - \ave{1}_{a,x}\ave{x}_{a} } $$

Now we apply this to the particular formula introduced above $$ f_a = \frac{ \sum_{bc} \phi_b\phi_c ( (r_c-r_a) f_b - (r_b-r_a) f_c )W_{ab}F_{ac,x}} { \sum_{de} \phi_d\phi_e (r_e - r_d) W_{ad} F_{ae,x} } \\ f_{a,x} = \frac{ \sum_{bc} \phi_b\phi_c ( f_c - f_b )W_{ab}F_{ac,x}} { \sum_{de} \phi_d\phi_e (r_e - r_d) W_{ad} F_{ae,x} } $$ ({\color{red} where}) $F_{ac,x}=W'_{ac}=W'(r_a-r_c)$

For any quantity $\psi$, one carries out the following average over the whole space ({\color{red} in order to obtain an quantity such as the ``action used in the Lagrangin formalism (?).}) $$ \sum_a \psi_a f_{a,x} = \sum_{bc} \phi_b\phi_c ( f_c - f_b ) \sum_a \psi_a \frac{ W_{ab}F_{ac,x}} { \sum_{de} \phi_d\phi_e (r_e - r_d) W_{ad} F_{ae,x} } $$ the second summation on the r.h.s. gives $$ \sum_a \psi_a \frac{ W_{ab}F_{ac,x}}{ \Omega_a } = \frac{\psi_b}{\phi_b\Omega_b} F_{bc,x} = (\frac{\psi_c}{\phi_c\Omega_c} W_{cb})_{,x} $$ where $$ \Omega_a = \sum_{de} \phi_d\phi_e (r_e - r_d) W_{ad} F_{ae,x} $$ $$ \sum_a \psi_a f_{a,x} = \sum_{bc} \frac{ \psi_b\phi_c }{\Omega_b} ( f_c - f_b ) F_{bc,x} $$ We set $\psi = \phi\Omega$ $$ \sum_a \phi_a\Omega_a f_{a,x} = \sum_{bc} \phi_b\phi_c ( f_c - f_b ) F_{bc,x} \\ =\sum_{bc} \phi_b\phi_c f_c F_{bc,x} - \sum_{bc} \phi_b\phi_c f_b F_{bc,x} \\ =\sum_{bc} \phi_b\phi_c f_c F_{bc,x} - \sum_{cb} \phi_c\phi_b f_c F_{cb,x} \\ =\sum_{bc} \phi_b\phi_c f_c F_{bc,x} + \sum_{cb} \phi_c\phi_b f_c F_{bc,x} \\ =2\sum_{bc} \phi_b\phi_c f_c F_{bc,x} = 2\sum_{c} \phi_c f_c \sum_b \phi_b F_{bc,x} $$ ({\color{red} a loose end?})

Now if we assume $$ \ave{ \frac{ f }{\rho} }_a = \sum_b \frac{f_b}{\rho^2_b} W_{ab} $$ $$ \ave{\frac{f}{\rho}}_a = \ave{1}_a (\frac{f}{\rho})_a + \ave{(x-r_a)_i}_a (\frac{f}{\rho})_{a,i} \\ \ave{\frac{f}{\rho}}_a = \ave{1}_a (\frac{f}{\rho})_a + \ave{(x-r_a)_i}_a (\frac{f_{a,i}}{\rho} - \frac{f_a \rho_{a,i}}{\rho^2_a}) $$ and $$ \ave{\frac{f}{\rho}}_{a,j} = \ave{1}_{a,j} (\frac{f}{\rho})_a + \ave{(x-r_a)_i}_{a,j} (\frac{f_{a,i}}{\rho} - \frac{f_a \rho_{a,i}}{\rho^2_a}) $$

In first order $$ \ave{\frac{f}{\rho}}_a = \ave{1}_a (\frac{f}{\rho})_a \\ f_a = \frac{\rho_a}{\ave{1}_a} \ave{\frac{f}{\rho}}_a $$ and

$$ [\begin{matrix} \ave{ \frac{f}{\rho} }_a\\ \ave{ \frac{f}{\rho} }_{a,j} \end{matrix}] = [\begin{matrix} \frac{\ave{1}_a}{\rho_a} - \frac{\ave{\Delta r_i}_a \rho_{a,i}}{\rho^2_a} & \frac{\ave{\Delta r_i}_a}{\rho_a} \\ \frac{\ave{1}_{a,j}}{\rho_a} - \frac{\ave{\Delta r_i}_{a,j} \rho_{a,i}}{\rho^2_a} & \frac{\ave{\Delta r_i}_{a,j}}{\rho_a} \end{matrix}] [\begin{matrix} f_a\\ f_{a,j} \end{matrix}] $$

$$ [\begin{matrix} f_a\\ f_{a,j} \end{matrix}] = \frac{1}{A} [\begin{matrix} \frac{\ave{\Delta r_i}_{a,j}}{\rho_a} & -\frac{\ave{\Delta r_i}_a}{\rho_a} \\ -\frac{\ave{1}_{a,j}}{\rho_a} + \frac{\ave{\Delta r_i}_{a,j} \rho_{a,i}}{\rho^2_a} & \frac{\ave{1}_a}{\rho_a} - \frac{\ave{\Delta r_i}_a \rho_{a,i}}{\rho^2_a} \end{matrix}] [\begin{matrix} \ave{ \frac{f}{\rho} }_a\\ \ave{ \frac{f}{\rho} }_{a,j} \end{matrix}] $$ where $$ A = (\frac{\ave{1}_a}{\rho_a} - \frac{\ave{\Delta r_i}_a \rho_{a,i}}{\rho^2_a}) \frac{\ave{\Delta r_i}_{a,j}}{\rho_a} \nqq - \frac{\ave{\Delta r_i}_a}{\rho_a}(\frac{\ave{1}_{a,j}}{\rho_a} - \frac{\ave{\Delta r_i}_{a,j} \rho_{a,i}}{\rho^2_a}) \nq = \frac{\ave{1}_a \ave{\Delta r_i}_{a,j} - \ave{1}_{a,j} \ave{\Delta r_i}_a}{\rho_a^2} $$ so that $$ f_a = \frac{1}{A}\left[\frac{\ave{\Delta r_i}_{a,j}}{\rho_a}\ave{ \frac{f}{\rho} }_a - \frac{\ave{\Delta r_i}_a}{\rho_a} \ave{ \frac{f}{\rho} }_{a,j}\right] \\ f_{a,j} =\frac{1}{A}\left[ -\frac{\ave{1}_{a,j}}{\rho_a}\ave{ \frac{f}{\rho} }_a + \frac{\ave{\Delta r_i}_{a,j} \rho_{a,i}}{\rho^2_a}\ave{ \frac{f}{\rho} }_a 

+ \frac{\ave{1}_a}{\rho_a}\ave{ \frac{f}{\rho} }_{a,j} - \frac{\ave{\Delta r_i}_a  \rho_{a,i}}{\rho^2_a}\ave{ \frac{f}{\rho} }_{a,j}\right]

$$



\section{Taylor revisited} In this section we explore another apporach to the Taylor expansion $$ \sum_b \phi_b f_b W_{ab} = \sum_b \phi_b f_a W_{ab} + \sum_b \phi_b ( r_b - r_a )_j f_{a,j} W_{ab} $$ one takes the derivative of the l.h.s. and expands afterwards $$ \sum_b \phi_b f_b F_{ab,i} = \sum_b \phi_b f_a F_{ab,i} + \sum_b \phi_b ( r_b - r_a )_jf_{a,j} F_{ab,i} $$ this leads to the following system $$ [\begin{matrix} \ave{\phi_b f_b}_a\\ \ave{\phi_b f_b}_{a,i} \end{matrix}] = [\begin{matrix} \ave{\phi_b}_a & \ave{\phi_b r_{bj}}_a - r_{aj} \ave{\phi_b}_a\\ \ave{\phi_b}_{a,i} & \ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i}\\ \end{matrix}] [\begin{matrix} f_a\\ f_{a,i} \end{matrix}] $$ ({\color{red} Though the notations are different, this is just TaylorSPH discussed previsously})

Now if we use $$ \phi_a \sum_b f_b W_{ab} = \phi_a \sum_b f_a W_{ab} + \phi_a \sum_b ( r_b - r_a )_j f_{a,j} W_{ab} \\ \phi_a \ave{f_b}_a = \phi_a f_a \ave{1}_a + \phi_a f_{a,j} (\ave{r_{bj}}_a - r_{aj}\ave{1}_a) $$ let us proceed by taking the derivative of the l.h.s. and expanding afterwards $$ \phi_{a,i} \sum_b f_b W_{ab} + \phi_a \sum_b f_b F_{ab,i} \nq = \phi_{a,i} \sum_b f_a W_{ab} + \phi_a \sum_b f_a F_{ab,i} \nqq + \phi_{a,i} \sum_b ( r_b - r_a )_jf_{a,j} W_{ab} + \phi_a \sum_b ( r_b - r_a )_jf_{a,j} F_{ab,i} $$ $$ \phi_{a,i} \ave{f_b} + \phi_a \ave{f_b}_{a,i} \nq = \phi_{a,i} f_a \ave{1}_a + \phi_a f_a \ave{1}_{a,i} \nqq + \phi_{a,i} f_{a,j} ( \ave{r_{bj}}_a - r_{aj}\ave{1}_a ) + \phi_a f_{a,j} ( \ave{r_{bj}}_{a,i} - r_{aj}\ave{1}_{a,i} ) $$ Using the original we have $$ \ave{f_b}_{a,i} = f_a \ave{1}_{a,i} + f_{a,j} ( \ave{r_{bj}}_{a,i} - r_{aj}\ave{1}_{a,i} ) $$ ({\color{red} This is just the results obtained in the above section of ``General Expression.})



\section{Taylor revisited II}

In this section we try to bring together the accuracy of the Taylor approach while keeping the conservative nature of the pressure gradient calculation. As we have seen in the previous sections, the Taylor approach provides remarkable results both for the function interpolation and the gradient estimation. However, in order to conserve momentum globally one needs to satisfy $$ \sum_a \frac{ \nu_a (\partial P)_a }{ \rho_a } = 0 $$ This property is satisfied by construction when using the variational principle which gives the expression for the gradient of pressure to be $$ (\partial P)_a = \rho_a \sum_b \nu_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2}) F_{ab} $$ On the other hand, this formula is found to perform very poorly compared to the Taylor corrected one.

Therefore, we try to include the conservative property into the Taylor correction method in order to achieve a high performance conservative expression for the gradient of pressure.

From the Taylor expansion we have $$ f_b = f_a + ( r_b - r_a ) \partial f_a $$ Now we want to introduce the conservative property, so we try $$ \sum_a \frac{ \nu_a (\partial f)_a }{ \rho_a } = \sum_a \frac{ \nu_a }{ \rho_a }\frac{ f_b - f_a }{ r_b - r_a } = 0 $$ This in turn let us write $$ f_b = \frac{ \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{ f_a }{ r_b - r_a } }{ \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } } $$

Now we apply the derivative to the whole formula to get $$ \partial_b f_b = -\frac{ \sum_{a\neq b} \frac{ \nu_a f_a }{ (\rho_a( r_b - r_a ) )^2 } \partial_b ( \rho_a( r_b - r_a ) ) }{ \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } } \nqq + \frac{ \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{ f_a }{ r_b - r_a } }{( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } )^2} \sum_{a\neq b} \frac{ \nu_a }{ (\rho_a(r_b - r_a))^2 } \partial_b ( \rho_a( r_b - r_a ) ) $$ To simplify the notation we write $$ \partial_b f_b = \frac{1}{( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } )^2} (-( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } ) \sum_{a\neq b} \frac{ \nu_a f_a }{ (\rho_a( r_b - r_a ) )^2 } \partial_b ( \rho_a( r_b - r_a ) ) \nqq + ( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{ f_a }{ r_b - r_a } ) \sum_{a\neq b} \frac{ \nu_a }{ (\rho_a( r_b - r_a ) )^2 }\partial_b ( \rho_a( r_b - r_a ) ) ) \nq = \frac{1}{( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } )^2} \sum_{a,c\neq b} \frac{ \nu_c \nu_a }{ \rho_c(r_b - r_c) (\rho_a( r_b - r_a ) )^2} ( f_c - f_a ) \partial_b ( \rho_a( r_b - r_a ) ) \nq = \frac{1}{( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } )^2} \sum_{a,c\neq b} \frac{ \nu_c \nu_a }{ \rho_c(r_b - r_c) (\rho_a( r_b - r_a ) )^2} ( f_c - f_a ) ( -( r_b - r_a )\nu_b F_{ab} + \rho_a ) $$

$$ \partial_b f_b = \frac{1}{( \sum_{a\neq b} \frac{ \nu_a }{ \rho_a }\frac{1}{ r_b - r_a } )^2} \sum_{a,c\neq b} \frac{ \nu_c \nu_a }{ \rho_c(r_b - r_c) (\rho_a( r_b - r_a ) )^2} ( f_c - f_a ) ( -( r_b - r_a )\nu_b F_{ab} + \rho_a ) $$

This calculation is not going anywhere... ({\color{red}Some note...}) It is understood that the ``conservation relation comes out naturally for any scalar function which vanishes (or have a constant value) at infinity. Since the derivative of a function in Taylor expansion method comes out naturally from the Taylor expansion, it is naturally to expect that the ``conservation (integral property of the derivative of a function) is accurate up to second order for a second order Taylor SPH method. In the context, maybe one can try to get the conservation relation rewritten according to the Taylor expansion, and compare the equality for each order.


\subsection{Symmetrization of TaylorSPH} We try a different approach to symmetrize the taylorSPH expression. Starting with the general taylor $$ [\begin{matrix} \ave{\phi_b f_b}_a\\ \ave{\phi_b f_b}_{a,i} \end{matrix}] = [\begin{matrix} \ave{\phi_b}_a & \ave{\phi_b r_{bj}}_a - r_{aj} \ave{\phi_b}_a\\ \ave{\phi_b}_{a,i} & \ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i}\\ \end{matrix}] [\begin{matrix} f_a\\ f_{a,i} \end{matrix}] $$ where $$ \ave{f}_a = \sum_b f_b W_{ab} \\ \ave{f}_{a,x} = \sum_b f_b F_{ab,x} $$ we have $$ [\begin{matrix} f_a\\ f_{a,i} \end{matrix}] = \Omega_a [\begin{matrix} \ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i} & - \ave{\phi_b r_{bj}}_a + r_{aj} \ave{\phi_b}_a\\ - \ave{\phi_b}_{a,i} & \ave{\phi_b}_a\\ \end{matrix}] [\begin{matrix} \ave{\phi_b f_b}_a\\ \ave{\phi_b f_b}_{a,i} \end{matrix}] $$ where $$ \Omega_a^{-1} = (\ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i})\ave{\phi_b}_a - (\ave{\phi_b r_{bj}}_a - r_{aj} \ave{\phi_b}_a)\ave{\phi_b}_{a,i} \nq = \ave{\phi_b r_{bj}}_{a,i}\ave{\phi_b}_a - \ave{\phi_b r_{bj}}_a\ave{\phi_b}_{a,i} \nq = \ave{\phi_b}_{a}^2 \partial_x( \frac{\ave{\phi_b r_{bj}}_{a}}{\ave{\phi_b}_a}) $$

which reads $$ f_a = \Omega_a ( (\ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i}) \ave{\phi f}_a - (\ave{\phi_b r_{bj}}_a - r_{aj} \ave{\phi_b}_a) \ave{\phi f}_{a,x} ) \\ f_{a,x} = \Omega_a ( \ave{\phi_b}_a \ave{\phi f}_{a,x} - \ave{\phi_b}_{a,x} \ave{\phi f}_a ) $$ Now we make the summation over the derivatives explict $$ f_a = \Omega_a \sum_b \phi_b ( (r_{bj} - r_{aj}) \ave{\phi f}_a - (\ave{\phi r_{j}}_a - r_{aj} \ave{\phi}_a ) f_b )F_{ab,x} \\ f_{a,x} = \Omega_a \sum_b\phi_b( \ave{\phi}_a f_b - \ave{\phi f}_a )F_{ab,x} $$ on the other hand $$ f_a = \Omega_a \sum_b \phi_b( (\ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i}) f_b - (r_{bj} - r_{aj}) \ave{\phi f}_{a,x} )W_{ab} \\ f_{a,x} = \Omega_a \sum_b( \ave{\phi}_a \phi_b f_b - \phi_b \ave{\phi f}_a )F_{ab,x} $$ Completely explicity $$ f_a = \Omega_a \sum_{bc} \phi_b\phi_c ( (r_{bj} - r_{aj}) \phi_c f_c - ( r_{cj} - r_{aj} ) f_b )W_{ac}F_{ab,x} \\ f_{a,x} = \Omega_a \sum_b\phi_b\phi_c( f_b - f_c )W_{ac} F_{ab,x} $$

Now let us say we want $$ \sum_a \phi_a f_{a,x} = 0 $$ $$ \sum_a \phi_a f_{a,x} = \sum_{ab} ( \phi_a\Omega_a\ave{\phi}_a \phi_b f_b - \phi_a\Omega_a\phi_b \ave{\phi f}_a )F_{ab,x} \nq = \sum_{abc} ( \phi_a\Omega_a\phi_c \phi_b f_b - \phi_a\Omega_a\phi_b \phi_c f_c )W_{ac}F_{ab,x} \nq = \sum_{abc} \phi_a\phi_b\phi_c ( \Omega_a f_b - \Omega_a f_c )W_{ac}F_{ab,x} \nq = \sum_{abc} \phi_a\phi_b\phi_c ( \Omega_a f_b W_{ac} + \Omega_b f_c W_{bc} )F_{ab,x} $$ We add and subtract the following term $$ \sum_a \phi_a f_{a,x} = \sum_{ab} ( \phi_a\Omega_a\ave{\phi}_a \phi_b f_b + \phi_b\Omega_b\ave{\phi}_b \phi_a f_a - \phi_b\Omega_b\ave{\phi}_b \phi_a f_a - \phi_a\Omega_a\phi_b \ave{\phi f}_a )F_{ab,x} \nq = \sum_{ab} ( \phi_a\Omega_a\ave{\phi}_a \phi_b f_b + \phi_b\Omega_b\ave{\phi}_b \phi_a f_a )F_{ab,x} - \sum_{ab} ( \phi_b\Omega_b\ave{\phi}_b \phi_a f_a + \phi_a\Omega_a\phi_b \ave{\phi f}_a )F_{ab,x} $$ the first term is symmetrical if we exchange $a$ and $b$ and therefore vanishes. We are left with $$ \sum_a \phi_a f_{a,x} = - \sum_{ab} ( \phi_b\Omega_b\ave{\phi}_b \phi_a f_a + \phi_a\Omega_a\phi_b \ave{\phi f}_a )F_{ab,x} \nq = - \sum_{abc} \phi_a\phi_b\phi_c\Omega_b f_a W_{cb} F_{ab,x} - \sum_{abc} \phi_a\phi_b\phi_c\Omega_a f_c W_{ca} F_{ab,x} \nq = - \sum_{abc} \phi_a\phi_b\phi_c\Omega_b f_a W_{cb} F_{ab,x} - \sum_{abc} \phi_a\phi_b\phi_c\Omega_b f_c W_{cb} F_{ba,x} \nq = \sum_{abc}\phi_a\phi_b\phi_c\Omega_b ( - f_a + f_c ) W_{cb} F_{ab,x} $$

Gave up at this point...



\subsection{} We try a slightly different system $$ f_b = f_a + f'_a (r_b-r_a) $$ the first equation is the same $$ \ave{f_b} = f_a\ave{1} + f'_a \ave{r_b-r_a} $$ but we modify the second equation so that $$ \ave{f_b - f_a}' = f'_a \ave{r_b - r_a}' $$

so the system looks like $$ [\begin{matrix} \ave{f_b}_a\\ \ave{f_b - f_a}' \end{matrix}] = [\begin{matrix} \ave{1}_a & \ave{\phi_b r_{bj}}_a - r_{aj} \ave{1}_a\\ 0 & \ave{\phi_b r_{bj}}'_{a} - r_{aj} \ave{1}'_{a}\\ \end{matrix}] [\begin{matrix} f_a\\ f'_{a} \end{matrix}] $$ and in 1d it reads as $$ f_a = \frac{ (\ave{\phi_b r_{bj}}_a - r_{aj} \ave{1}_a) \ave{a} }{ \ave{1}_a( \ave{\phi_b r_{bj}}'_{a} - r_{aj} \ave{1}'_{a} ) } $$

$$ f_a = \frac{ \ave{f}_a }{ \ave{1}_a } - \frac{\ave{ \Delta r}_{a} }{ \ave{\Delta r}'_{a} \ave{1}_a } \ave{f_b - f_a}'_{a} \\ f_{a,x} = \frac{ \ave{f_b - f_a }'_{a} }{ \ave{\Delta r}'_{a} } $$ this last expression can be rewritten as $$ \sum_a f_{a,x} = \sum_a\sum_b \frac{ (f_b - f_a) }{ \Omega_a } F_{ab} $$



\subsection{another try} $$ \ave{\phi_b f_b} = f_a\ave{\phi_b} + f'_a \ave{\phi_b (r_b-r_a) } $$ $$ f_b\ave{\phi_a}_b = \ave{\phi_a f_a}_b + \ave{ \phi_a f'a (r_b - r_a) } $$



\subsection{} The quantities calculated for the particle $A$ are $$ f_A = \Omega_A \sum_b \phi_b ( (r_{bj} - r_{Aj}) \ave{\phi f}_A - (\ave{\phi r_{j}}_A - r_{Aj} \ave{\phi}_A ) f_b )F_{Ab,x} \\ f_{A,x} = \Omega_A \sum_b\phi_b( \ave{\phi}_A f_b - \ave{\phi f}_A )F_{Ab,x} $$ now we separate the term corresponding to particle $B$ from the sum $$ f_A = \Omega_A ( \phi_B ( (r_{Bj} - r_{Aj}) \ave{\phi f}_A - (\ave{\phi r_{j}}_A - r_{Aj} \ave{\phi}_A ) f_B )F_{AB,x} \nqq + \sum_{b\neq B} \phi_b ( (r_{bj} - r_{Aj}) \ave{\phi f}_A - (\ave{\phi r_{j}}_A - r_{Aj} \ave{\phi}_A ) f_b )F_{Ab,x} ) \\ f_{A,x} = \Omega_A (\phi_B( \ave{\phi}_A f_B - \ave{\phi f}_A )F_{AB,x} + \sum_{b \neq B}\phi_b( \ave{\phi}_A f_b - \ave{\phi f}_A )F_{Ab,x}) $$ now we compute the same for particle $B$ $$ f_B = \Omega_B ( \phi_A ( (r_{Aj} - r_{Bj}) \ave{\phi f}_B - (\ave{\phi r_{j}}_B - r_{Bj} \ave{\phi}_B ) f_A )F_{BA,x} \nqq + \sum_{b\neq A} \phi_b ( (r_{bj} - r_{Bj}) \ave{\phi f}_B - (\ave{\phi r_{j}}_B - r_{Bj} \ave{\phi}_B ) f_b )F_{Bb,x} ) \\ f_{B,x} = \Omega_B (\phi_A( \ave{\phi}_B f_A - \ave{\phi f}_B )F_{BA,x} + \sum_{b \neq A}\phi_b( \ave{\phi}_B f_b - \ave{\phi f}_B )F_{Bb,x}) $$

An ad hoc symmetrization of the derivative takes the form of $$ f_{a,x} = \frac{1}{2} \sum_b\phi_b( \Omega_a\ave{\phi}_a f_b + \Omega_b\ave{\phi}_b f_a - \Omega_a\ave{\phi f}_a - \Omega_b\ave{\phi f}_b )F_{ab,x} $$

$$ f_{a,x} = \frac{1}{2} ( \Omega_a\ave{\phi}_a \ave{\phi_b f_b}_{a,x} + \ave{\Omega_b\phi_b \ave{\phi}_b}_{a,x} f_a - \Omega_a\ave{\phi f}_a\ave{\phi}_{a,x} - \ave{\phi_b\Omega_b\ave{\phi f}_b}_{a,x} ) $$




\subsection{From the integral formulation}

In the standard SPH we have that the integral of a quantity is given by $$ \sum_a \nu_a \frac{f_a}{\rho_a} $$ In this case, in order to ensure that a quantity is conservative, one needs that $$ \sum_a \nu_a \frac{ f_{a,i} }{\rho_a} = 0 $$ Now if we manipulate this expression we arrive at $$ \sum_a \nu_a ( (\frac{ f_{a} }{\rho_a})_{,i} + \frac{f_a}{\rho_a^2}\rho_{a,i} ) = 0 $$ from where we can identify $$ \frac{ f_{a,i} }{\rho_a} = (\frac{ f_{a} }{\rho_a})_{,i} + \frac{f_a}{\rho_a^2}\rho_{a,i} $$ or $$ f_{a,i} = \rho_a \sum_b \nu_b ( \frac{ f_{b} }{\rho_b^2} + \frac{f_a}{\rho_a^2}) F_{ab,i} $$

However, if we take the standard SPH interpolation $$ f_a = \sum_b \frac{\nu_b f_b}{\rho_b}W_{ab} $$ and we use the obove definition of the integral we get $$ \sum_a \nu_a \frac{f_a}{\rho_a} = \sum_b \frac{\nu_b f_b}{\rho_b} \sum_a \frac{\nu_a }{\rho_a} W_{ab} \nq = \sum_b \frac{\nu_b f_b I_b}{\rho_b} $$ which shows the inconsistency since $$ \sum_a \frac{\nu_a }{\rho_a} W_{ab} \neq 1 $$

In order to get a consistent definition we use a general SPH formulation $$ f_a = \psi_a \sum_b \phi_b f_b W_{ab} $$ so that $$ \sum_a \Omega_a f_a = \sum_b \phi_b f_b \sum_a \Omega_a \psi_a W_{ab} $$ in order to make both sides equal we have to satisfy $$ \frac{\Omega_a}{\phi_a} = \sum_b \Omega_b \psi_b W_{ab} $$ if we take $$ \Omega_a = \frac{\nu_a}{\psi_a} $$ we have $$ \Omega_a = \phi_a \rho_a \\ \frac{ \nu_a }{ \phi_a\psi_a } = \sum_b \nu_b W_{ab} = \rho_a \\ \psi_a = \frac{ \nu_a }{\phi_a\rho_a} $$ Now if we compute $$ \rho_a = \psi_a \sum_b \phi_b \rho_b W_{ab} \nq = \frac{ \nu_a }{\phi_a\rho_a} \sum_b \phi_b \rho_b W_{ab} $$ we get that $$ \phi_a = \frac{\nu_a}{\rho_a} $$

Summarizing $$ \psi_a = 1 \\ \Omega_a = \nu_a $$

$$ f_a = \sum_b \frac{\nu_b f_b}{\rho_b} W_{ab} $$ the integral form is then $$ \sum_a \nu_a f_a $$ which does not correspond to the original formula for the integral

We list here the constraints for a general SPH formalution including the unity description $$ f_a = \psi_a \sum_b \phi_b f_b W_{ab} \\ \rho_a = \psi_a \sum_b \phi_b \rho_b W_{ab} = \sum_b \nu_b W_{ab} \\ \psi_a \sum_b \phi_b W_{ab} = 1 $$ from the last one we conclude that $$ \psi_a = \frac{1}{\sum_b \phi_b W_{ab}} $$ so that $$ \sum_b \phi_b \rho_b W_{ab} = \sum_b \phi_b W_{ab} \sum_b \nu_b W_{ab} $$ this means that with the current approach we cannot satisfy unity and the self-consistecy of $\rho$ at the same time.

Let us try yet another approach $$ f_a = \psi_a \sum_b \phi_b f_b W_{ab} + \psi'_a \sum_b \phi'_b f_b W_{ab} $$ in this case we have $$ \rho_a = \psi_a \sum_b \phi_b \rho_b W_{ab} + \psi'_a \sum_b \phi'_b \rho_b W_{ab} \\ \psi_a \sum_b \phi_b W_{ab} + \psi'_a \sum_b \phi'_b W_{ab} = 1 $$

Let us go back to general formula for the integral $$ \sum_a \Omega_a^\lambda f_a $$ now let us assume we want to have a conservative expression $$ \sum_a \Omega_a^\lambda f'_a = 0 \nq \sum_a (\Omega_a^\lambda f_a)' - \lambda \Omega_a^{\lambda-1} f_a \Omega'_a = 0 $$ so that $$ f'_a = \frac{1}{ \Omega_a^\lambda } ( (\Omega_a^\lambda f_a)' - \lambda \Omega_a^{\lambda-1} f_a \Omega'_a ) $$ Now we can use TaylorSPH expressions for these quantites $$ f'_a = \frac{ \ave{1}\ave{f}' - \ave{1}'\ave{f}}{\ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r}} $$ $$ (\Omega^\lambda f)'_a = \frac{ \ave{1}\ave{\Omega^\lambda f}' - \ave{1}'\ave{ \Omega^\lambda f }}{\ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r}} $$

$$ f'_a = \frac{1}{\Omega_a^\lambda }( \frac{ \ave{1}\ave{\Omega^\lambda f}' - \ave{1}'\ave{ \Omega^\lambda f }}{\ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r}} - \lambda \Omega^{\lambda-1}_a f_a \frac{ \ave{1}\ave{\Omega}' - \ave{1}'\ave{ \Omega }}{\ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r}} ) $$

In case we use $$ \Omega = \frac{1}{\rho} $$ we get $$ f'_a = \frac{\rho_a}{ \ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r} }( \ave{1}\ave{\frac{f}{\rho} }' - \ave{1}'\ave{ \frac{f}{\rho} } + \frac{f_a}{\rho^2}( \ave{1}\ave{\Omega}' - \ave{1}'\ave{ \rho } ) ) \nq = \frac{\rho_a}{ \ave{1}\ave{\Delta r}' - \ave{1}'\ave{\Delta r} }( \ave{1}(\ave{ \frac{f}{\rho} }' + \frac{f_a}{\rho^2}\ave{ \rho }' ) - \ave{1}'\ave{ \frac{f}{\rho} } + \frac{f_a}{\rho^2}( \ave{1}\ave{ \rho }' - \ave{1}'\ave{ \rho } ) ) $$



\subsection{Taylor revisited again}

$$ f_b = f_a + (r_b - r_a) f'_a $$ so we rewrite $$ \rho^\lambda_b f_b = \rho^\lambda_a f_a + (r_b - r_a) (\rho^\lambda_a f_a)' \nq = \rho^\lambda_a f_a + \rho^\lambda_a (r_b - r_a) ( f'_a + \lambda\frac{f_a}{\rho_a}\rho_a') \nq = \rho^\lambda_a( 1 + \lambda\frac{(r_b - r_a)}{\rho_a}\rho_a' )f_a + \rho^\lambda_a (r_b - r_a) f'_a $$

$$ \frac{\rho^\lambda_b f_b - \rho^\lambda_a f_a}{(r_b - r_a)} - \lambda\rho^\lambda \frac{f_a}{\rho_a}\rho_a' = \rho^\lambda_a f'_a $$

$$ \sum_b \rho^{\lambda-1}_b f_b W_{ab} = \rho^\lambda_a f_a( \sum_b \frac{W_{ab}}{\rho_b} + \lambda\rho_a' \sum_b\frac{(r_b - r_a)}{\rho_a\rho_b}W_{ab} ) + \rho^\lambda_a f'_a \sum_b\frac{(r_b - r_a)}{\rho_b}W_{ab} $$

$$ [ \begin{matrix} \sum_b \rho^{\lambda-1}_b f_b W_{ab} \\ \sum_b \rho^{\lambda-1}_b f_b F_{abi} \end{matrix} ] = [ \begin{matrix} \rho^\lambda_a( \sum_b \frac{W_{ab}}{\rho_b} + \lambda\rho_a' \sum_b\frac{(r_b - r_a)}{\rho_a\rho_b}W_{ab} ) & \rho^\lambda_a\sum_b\frac{(r_b - r_a)}{\rho_b}W_{ab} \\ \rho^\lambda_a( \sum_b \frac{F_{ab}}{\rho_b} + \lambda\rho_a' \sum_b\frac{(r_b - r_a)}{\rho_a\rho_b}F_{ab} ) & \rho^\lambda_a\sum_b\frac{(r_b - r_a)}{\rho_b}F_{ab} \end{matrix} ] [ \begin{matrix} f_a \\ f'_a \end{matrix} ] $$

in 1d the determinant is $$ \det(M) = $$


\section{The Lagrangian Formalism} In this section we explore different Lagrangian formulations to extract a consistent and efficient equation of motion from the SPH.




\section{explicit errors} with the explicit errors introduced by the truncation we have $$ [\begin{matrix} \ave{\phi_b f_b}_a\\ \ave{\phi_b f_b}_{a,i} \end{matrix}] = [\begin{matrix} \ave{\phi_b}_a & \ave{\phi_b r_{bj}}_a - r_{aj} \ave{\phi_b}_a\\ \ave{\phi_b}_{a,i} & \ave{\phi_b r_{bj}}_{a,i} - r_{aj} \ave{\phi_b}_{a,i}\\ \end{matrix}] [\begin{matrix} f_a\\ f_{a,i} \end{matrix}] + [\begin{matrix} e_1\\ e_2 \end{matrix}] $$


\subsection{The standard SPH} First we apply the variational method to the standard SPH formulation which assumes that $$ \int d^3 f(x) = \sum_b \frac{\nu_b f_b}{ \rho_b } \\ \int d^3 f(x) W(x-r_a) = \sum_b \frac{\nu_b f_b}{ \rho_b } W( r_b - r_a ) \\ \varepsilon_a = \sum_b \frac{\nu_b\varepsilon_b}{ \rho_b } W_{ab} \\ \varepsilon_{a,i} = \sum_b \frac{\nu_b\varepsilon_b}{ \rho_b } F_{ab,i} \\ \rho_a = \sum_b \nu_b W_{ab} \\ \rho_{a,i} = \sum_b \nu_b F_{ab,i} $$

The action is $$ L = \int dx^3 \sqrt{-g} \varepsilon(x) \nq = \sum_a \frac{\sqrt{-g}_a\nu_a\varepsilon_a}{ \rho_a } \int dx^3 W(x-r_a) \nq = \sum_a \frac{\sqrt{-g}_a\nu_a\varepsilon_a}{ \rho_a } $$ in the last line we used the assumption that the kernel is normalized.

Now we compute ({\color{red} No $\sqrt{-g}$ in the following expressions?}) $$ \partial_{r_b} L = \sum_a \frac{\nu_a}{\rho_a^2}(\rho_a\partial_{r_b}\varepsilon_a - \varepsilon_a \partial_{r_b}\rho_a) \nq = \sum_a \frac{\nu_a}{\rho_a^2}(\rho_a\frac{\partial \varepsilon_a }{\partial \rho_a } - \varepsilon_a ) \partial_{r_b}\rho_a \nq = \sum_a A_a \partial_{r_b}\rho_a $$

$$ \partial_{r_b}\rho_a = \sum_c \nu_c \partial_{r_b} W_{ac} \nq = \sum_c \nu_c \frac{\partial W_{ac}}{\partial (r_a - r_c)} \partial_{r_b} (r_a - r_c) \nq = \sum_c \nu_c \frac{\partial W_{ac}}{\partial (r_a - r_c)} (\delta_{ab} - \delta_{bc}) $$

$$ \frac{\partial W_{ac}}{\partial (r_a - r_c)} = \frac{W'_{ac}}{h} \frac{\partial \abs{r_a - r_c}}{\partial (r_a - r_c)} = \frac{W'_{ac}}{h} \frac{r_a - r_c}{\abs{r_a - r_c}} = F_{ac} = - F_{ca} $$

$$ \partial_{r_b} L = \sum_a A_a \sum_c \nu_c F_{ac} (\delta_{ab} - \delta_{bc}) \nq = \sum_{ac}\nu_c A_a F_{ac} (\delta_{ab} - \delta_{bc}) \nq = \sum_{c}\nu_c A_b F_{bc} - \sum_{a}\nu_b A_a F_{ab} \nq = \sum_{a}(\nu_a A_b F_{ba} - \nu_b A_a F_{ab} ) \nq = \sum_{a}(\nu_a A_b + \nu_b A_a ) F_{ba} $$ which is by definition symmetric.

On the other hand $$ \partial_{\dot{r}_b} L = \sum_a \frac{\nu_a}{\rho_a}\partial_{\dot{r}_b}\varepsilon_a 

= \sum_a \frac{\nu_a}{\rho_a}\frac{\partial \varepsilon_a}{\partial \dot{r}_a} \delta_{ab}
= \frac{\nu_b}{\rho_b}\frac{\partial \varepsilon_b}{\partial \dot{r}_b}

$$

Therefore one arrives at the equation of motion (EoM): $$ \frac{d}{dt} ( \frac{\nu_b}{\rho_b}\frac{\partial \varepsilon_b}{\partial \dot{r}_b} ) 

= \sum_{a}(\nu_a A_b  + \nu_b A_a ) \frac{\partial W_{ba}}{\partial (r_b - r_a)}

$$

\subsection{Working with the densities} We can work straight with the densitites $$ \mathcal{L}(r_a) = \sqrt{-g} \varepsilon(r_a) = \sum_b \frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b } W_{ab} $$

$$ \partial_{\dot{r}_c} \mathcal{L}_a = \frac{\sqrt{-g}_c\nu_c}{ \rho_c } \frac{\partial \varepsilon_c}{\partial \dot{r}_c} W_{ac} $$

$$ \partial_{r_c} \mathcal{L}_a = \sum_b \partial_{r_c}(\frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b }) W_{ab} 

+ \sum_b \frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b } \partial_{r_c} W_{ab}

\nq = \sum_b A_b W_{ab} \partial_{r_c}\rho_b 

+ \sum_b \frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b }  \partial_{r_c} W_{ab}

\nq = \sum_b A_b W_{ab} \sum_d \nu_d F_{bd} (\delta_{bc} - \delta_{dc}) 

+ \sum_b \frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b } F_{ab} (\delta_{ac} - \delta_{bc})

\nq = \sum_{b}( A_c W_{ac} \nu_b F_{cb} - A_b W_{ab} \nu_d F_{bc}) 

+ \delta_{ac} \sum_b \frac{\sqrt{-g}_b\nu_b\varepsilon_b}{ \rho_b } F_{cb} - \frac{\sqrt{-g}_c\nu_c\varepsilon_c}{ \rho_c } F_{ac}

$$ Mathematically, this does not work because of the extra derivatives of the kernel. This, in turn is due to the fact that we applied the variational principle to the densities. We must therefore, apply the principle to the integrated quantities, as in the previous case to avoid this problem. ({\color{red} It does not work since we tried to apply the variational principle to a certain point in space, in stead of the action, the latter is not a function of coordinate but an integral of the whole system. The action is an integral over space (in the case of a field) and time, and summation over all degrees of freedom of the system.})


\subsection{Self-consistency}

SPH formalism as an approximation, one may simply apply it, namely, to obtain an EoM of the continue in terms of SPH degree of freedom. However, it was shown by Kodama and Hama that the EoM obtained above is in fact consistent with another approach, namely, apply the variational principle directly to the Lagrangian expressed in terms of SPH particles in first place. This is shown below. This is not too surprising due to the fact that hydrodynamical equation itself can be derived using variational principle. Unfortunately, it seems that the abve conclusion is only valid for the case of ideal hydrodynamics.

In the discussions of previous sections, we have introduced some more general forms for the SPH formalism, it is therefore interesting to see how the variational principle is affected. The study can be carried out in two approaches. Firstly, we start from hydrodynamical equation, which is known to be exact, and substitute physical quantites by using the above expressions. On the other hand, one may start from the Lagragian, which its form yet to be determined (to the least, it should be consistent with the corresponding SPH formalism), calculate its variation and take care of extra non-vanishing derivatives, and hopefully, one arrives at the same destiny.

$$ \int d^3x f(x) W(x-r_a) = \frac{1}{\rho_a} \sum_b \nu_b f_b W( r_b - r_a ) \\ \int d^3x W(x-r_a) = \frac{1}{\rho_a} \sum_b \nu_b W( r_b - r_a ) = 1 $$

by integrating $x'$ first we get $$ \int d^3x f(x) = \int d^3x d^3x' f(x) W(x-x') \nq = \int d^3x \frac{f(x)}{ \rho(x) } \sum_a \nu_a W( x - r_a ) $$ on the other hand if we integrate $x$ first $$ \int d^3x f(x) = \int d^3x d^3x' f(x) W(x-x') \nq = \int d^3x' \frac{1}{\rho(x')}\sum_b \nu_b f_b W( r_b - x' ) \nq = \sum_b \nu_b f_b \frac{1}{\rho_b} \sum_a \frac{\nu_a}{\rho_a} W( r_b - r_a ) \nq = \sum_b \frac{\nu_b I_b f_b}{\rho_b} $$

$$ \partial_{\dot{r}_a} L = \sqrt{-g}_a \frac{\nu_a I_a }{\rho_a} \partial_{\dot{r}_a} \varepsilon_a $$

({\color{red} Maybe the expression of the Lagragian should contain already the first-order derivative.})

$$ \partial_{r_a} L = \sum_b \frac{\nu_b}{\rho_b^2} (\rho_b I_b \partial_{r_a} f_b + \rho_b f_b \partial_{r_a} I_b \nqq - I_b f_b \partial_{r_a} \rho_b ) \nq = \sum_b \frac{\nu_b}{\rho_b^2} (I_b (\rho_b \frac{\partial f_b}{\partial \rho_b} - f_b ) \partial_{r_a} \rho_b + \rho_b f_b \partial_{r_a} I_b ) $$

$$ \sum_b A_b \partial_{r_a} \rho_b = \sum_b A_b \sum_c \nu_c F_{bc} (\delta_{ab} - \delta_{ac}) \nq = \sum_b A_a \nu_b F_{ab} - A_b \nu_a F_{ba} \nq = \sum_b (A_a \nu_b + A_b \nu_a) F_{ab} $$

$$ \sum_b A_b \partial_{r_a} I_b = \sum_b A_b \partial_{r_a} \sum_c \frac{\nu_c}{\rho_c} W_{bc} \nq = \sum_b A_b \sum_c \frac{\nu_c}{\rho_c} \partial_{r_a} W_{bc} \nqq -\sum_b A_b \sum_c \frac{\nu_c}{\rho_c^2} W_{bc} \partial_{r_a} \rho_c \nq = \sum_b A_b \sum_c \frac{\nu_c}{\rho_c} F_{bc} (\delta_{ab} - \delta_{ac}) \nqq -\sum_b A_b \sum_c \frac{\nu_c}{\rho_c^2} W_{bc} \sum_d \nu_d F_{cd}(\delta_{ac} - \delta_{ad}) \nq = \sum_b (A_a\frac{\nu_b}{\rho_b} F_{ab} - A_b\frac{\nu_a}{\rho_a} F_{ba}) \nqq -\sum_b A_b \sum_c (\frac{\nu_a}{\rho_a^2} W_{ba}\nu_d F_{ac} - \frac{\nu_c}{\rho_c^2} W_{bc}\nu_a F_{ca}) \nq = \sum_b (A_a\frac{\nu_b}{\rho_b} F_{ab} - A_b\frac{\nu_a}{\rho_a} F_{ba}) \nqq -\sum_b A_b \sum_c (\frac{\nu_a}{\rho_a^2} W_{ba}\nu_c F_{ac} - \frac{\nu_c}{\rho_c^2} W_{bc}\nu_a F_{ca}) \nq = \sum_b (A_a\frac{\nu_b}{\rho_b} + A_b\frac{\nu_a}{\rho_a} )F_{ab} \nqq -\sum_b A_b \sum_c (\frac{\nu_a\nu_c}{\rho_a^2} W_{ba} + \frac{\nu_a\nu_c}{\rho_c^2} W_{bc})F_{ac} $$

$$ \partial_{r_a} L \nq = \sum_b \frac{\nu_b}{\rho_b^2} (I_b (\rho_b \frac{\partial f_b}{\partial \rho_b} - f_b ) \partial_{r_a} \rho_b \nqq + \rho_b f_b \partial_{r_a} I_b ) \nq = \sum_b \nu_b\nu_a ( \frac{I_a}{\rho_a^2} (\rho_a \frac{\partial f_a}{\partial \rho_a} - f_a ) + \frac{I_b}{\rho_b^2} (\rho_b \frac{\partial f_b}{\partial \rho_b} - f_b)) F_{ab} \nqq + \sum_b (\frac{\nu_a\nu_b f_a}{\rho_a\rho_b} + \frac{\nu_b\nu_a f_b}{\rho_b\rho_a} )F_{ab} \nqq -\sum_b \frac{\nu_b f_b}{\rho_b} \sum_c (\frac{\nu_a\nu_c}{\rho_a^2} W_{ba} + \frac{\nu_a\nu_c}{\rho_c^2} W_{bc})F_{ac} $$

we note that, $$ \rho_a \frac{\partial \varepsilon_a}{\partial \rho_a} - \varepsilon_a = \sigma_a \frac{\partial \varepsilon_a}{\partial \sigma_a} - \varepsilon_a = \omega_a - \varepsilon_a = P_a $$ so that

$$ \partial_{r_a} L \nq = \sum_b \nu_b\nu_a ( \frac{I_a P_a}{\rho_a^2} + \frac{I_bP_b}{\rho_b^2}) F_{ab} \nqq + \sum_b \frac{\nu_b\nu_a}{\rho_a\rho_b}( \varepsilon_a + \varepsilon_b )F_{ab} \nqq -\sum_c \nu_a\nu_c (\frac{1}{\rho_a^2}\sum_b \frac{\nu_b f_b}{\rho_b} W_{ba} \nqq + \frac{1}{\rho_c^2}\sum_b \frac{\nu_b f_b}{\rho_b} W_{bc})F_{ac} $$

where $$ \ave{\frac{f}{\rho}}_a = \frac{1}{\rho_a} \sum_b \frac{\nu_b f_b}{\rho_b} W_{ab} $$

so $$ \partial_{r_a} L \nq = \sum_b \nu_b\nu_a ( \frac{I_a P_a}{\rho_a^2} + \frac{I_bP_b}{\rho_b^2}) F_{ab} \nqq + \sum_b \frac{\nu_b\nu_a}{\rho_a\rho_b}( \varepsilon_a + \varepsilon_b )F_{ab} \nqq -\sum_b \nu_a\nu_b (\frac{1}{\rho_a}\ave{\frac{\varepsilon}{\rho}}_a + \frac{1}{\rho_b}\ave{\frac{\varepsilon}{\rho}}_b)F_{ab} $$

\subsection{Self-consistency} $$ \int d^3x f(x) W(x-r_a) = \sum_b \frac{\nu_b f_b}{ \rho_b } W( r_b - r_a ) \\ \int d^3x W(x-r_a) = \sum_b \frac{ \nu_b }{ \rho_b } W( r_b - r_a ) $$

therefore if we integrate $x'$ first $$ \int d^3x f(x) = \int d^3x d^3x' f(x) W(x-x') \nq = \int d^3x f(x) \sum_a \frac{ \nu_a }{ \rho_a } W( x - r_a ) \nq = \sum_a \frac{ \nu_a }{ \rho_a } \int d^3x f(x) W( x - r_a ) \nq = \sum_a \frac{ \nu_a }{ \rho_a } \sum_b \frac{\nu_b f_b}{ \rho_b } W( r_b - r_a ) \nq = \sum_{ab} \frac{ \nu_a\nu_b }{ \rho_a\rho_b } f_b W( r_b - r_a ) \nq = \sum_{b} \frac{ \nu_b I_b }{ \rho_b } f_b $$ if we integrate $x$ first we have $$ \int d^3x f(x) = \int d^3x d^3x' f(x) W(x-x') \nq = \sum_b \frac{\nu_b f_b}{ \rho_b } \int d^3x' W( r_b - x' ) \nq = \sum_b \frac{\nu_b f_b}{ \rho_b } \sum_a \frac{ \nu_a }{ \rho_a } W( r_b - r_a ) \nq = \sum_{ab} \frac{\nu_a\nu_b }{ \rho_a\rho_b } f_b W( r_b - r_a ) \nq = \sum_{b} \frac{\nu_b I_b }{ \rho_b } f_b $$

now we apply this to the action $$ L = \int dx^3 \sqrt{-g} \varepsilon(x) \nq = \sum_{ab} \frac{\nu_a\nu_b }{ \rho_a\rho_b } \sqrt{-g}_a \varepsilon_a W( r_a - r_b ) $$ we have a kernel in the lagrangean, but now there is a double sum, so the deltas are all going to be computed.

$$ \partial_{\dot{r}_c} L = \sum_{ab} \frac{\nu_a\nu_b }{ \rho_a\rho_b } \sqrt{-g}_a \partial_{\dot{r}_c} \varepsilon_a W( r_a - r_b ) \nq = \sum_{b} \frac{\nu_c\nu_b }{ \rho_c\rho_b } \sqrt{-g}_c \partial_{\dot{r}_c} \varepsilon_c W( r_c - r_b ) \nq = \frac{\nu_c}{\rho_c } \sqrt{-g}_c \partial_{\dot{r}_c} \varepsilon_c \sum_{b} \frac{\nu_b }{\rho_b } W( r_c - r_b ) \nq = \frac{\nu_c I_c}{\rho_c } \sqrt{-g}_c \partial_{\dot{r}_c} \varepsilon_c $$

$$ \partial_{r_c} L = $$


\subsection{}

we aproximate the lagrangean as $$ I = \int d^4 \varepsilon(x) = \int dt\int d^3 \varepsilon(x) = \int dt\int d^3 \frac{\sum_a \nu_a \varepsilon_a \rho_a^{\lambda} W(x - x_a)}{\sum_a \nu_a \rho_a^{\lambda} W(x - x_a)} \\ = \int dt \sum_a \nu_a \varepsilon_a \rho_a^{\lambda} \int d^3x \frac{W(x - x_a)}{\sum_b \nu_b \rho_b^{\lambda} W(x - x_b)} $$

\subsection{Intrinsic relativistic} In this section we shall explore the relativistic aspect of the distance measured between particles. In a non-relativistic case, the kernel depends on the relative distance of the test particles as $$ q_{ab} = \frac{\abs{x_a - x_b}}{h} = \sqrt{ \frac{(x_a - x_b)_i}{h} \frac{(x_a - x_b)_i}{h} } $$

In the relativistic case, however, this distance corresponds to a very specific frame of reference, the laboratory frame. In the local rest frame of each particle, effects of Lorentz contraction are important to the direction of relative motion. We can decompose the distance between two particles in $$ \abs{\Delta x} \nq = g_{\mu\nu}\Delta x^\mu \Delta x^\nu = u_\mu u_\nu \Delta x^\mu \Delta x^\nu \nqq + (g_{\mu\nu} - u_\mu u_\nu) \Delta x^\mu \Delta x^\nu \nq = (u_\mu \Delta x^\mu)^2 + \Delta_{\mu\nu} \Delta x^\mu \Delta x^\nu $$ The first term in the rhs is the part which is parallel to the relative proper velocity $u^\mu$. The second term represents the perpendicular part. Therefore, the first part must be Lorentz contracted.

Let us use the example of a particle at the origin moving in the $x$-direction with $u_x=\sqrt{0.5}$. Compared to a particle at rest at $x=1$ its distance is without relativity corrections $$ \abs{\Delta x} = 1 $$ Taking into account the lorentz factor, its distance should be contracted with $\gamma = \sqrt{1 + u^2}$. $$ \abs{\Delta x} = \frac{1}{\sqrt{1 + 0.5}} \approx 0.8 $$ However, another particle at position $(1,1)$ give the following values without relativity $$ \abs{\Delta x} = \sqrt{2} \approx 1.41 $$ With it $$ \abs{\Delta x} = \sqrt{\frac{\Delta x^2}{\gamma^2} + \Delta y^2} \nqq = \sqrt{ \frac{1}{\sqrt{1 + 0.5}} + 1 } \approx 1.29 $$


\end{document}

Retrieved from "https://wlqian.shoutwiki.com/w/index.php?title=Philipe%27s_Notes_on_SPH_algorithm&oldid=598"