Changes between Version 6 and Version 7 of CodeExplanation/SourceTerms

Timestamp:

01/12/12 13:23:03 (14 years ago)

Author:

Jonathan

Comment:

—

CodeExplanation/SourceTerms

@@ -17 +17 @@
 [[latex($\Delta p=-\frac{\Delta t}{\Delta x} \left(F_{i+1/2}-F_{i-1/2} \right )$)]]  which is just [[latex($\dot{p} = f_{grav}$)]] where [[latex($f_{grav}=-\nabla \cdot F$)]] which requires [[latex($\nabla \cdot F=\rho \nabla \phi$)]].  Since [[latex($\nabla^2 \phi = 4\pi G (\rho-\bar{\rho})$)]] we can substitute for [[latex($\rho$)]] and we have [[latex($\nabla \cdot F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi$)]] which is equivalent (in 1D) to [[latex($\nabla \cdot \left [\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G}  + \bar{\rho} \phi \right ]$)]] where we can identify the equivalent momentum flux tensor as [[latex($F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi$)]]
 
-* In more than 1D, [[latex($\nabla \cdot \left [\frac{1}{2} \frac{\left (\nabla \phi : \nabla \phi \right )}{4 \pi G}  + \bar{\rho} \phi I \right ]_i =\partial_j (\partial_j \phi \partial_i \phi) +\bar{\rho} \partial_i \phi$)]]
+* In more than 1D we have 
+
+
+{{{
+#!latex
+\begin{align*}
+\partial_t p_i&=\rho \partial_i \phi \\ 
+& = \left (\frac{\partial_j(\partial_j \phi)}{4 \pi G}+\bar{\rho} \right ) \partial_i \phi \\ 
+& = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_j (\partial_i \phi )+\partial_i (\bar{\rho}\phi) \\
+& = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_i (\partial_j \phi )+\partial_i (\bar{\rho}\phi) \\
+& = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - \frac{1}{2}\partial_i(\partial_j \phi \partial_j \phi)+\partial_i (\bar{\rho}\phi) \\
+& = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - \partial_i \left (\frac{\partial_j \phi \partial_j}{2} \phi+\bar{\rho}\phi \right) \\
+& = \partial_j \left ( \partial_j \phi \partial_i \phi \right) -  \partial_j \left ( \delta ^i_j\left(\frac{\partial_k \phi \partial_k \phi}{2}+\bar{\rho}\phi \right) \right) \\
+& = \partial_j T_{ij} \mbox{ where } T_{ij} = \left ( \partial_j \phi \partial_i \phi \right) -    \delta ^i_j\left(\frac{\partial_k \phi \partial_k \phi}{2}+\bar{\rho}\phi \right) \\
+\end{align*}
+}}}
 
 * But casting it this way gives you strict momentum conservation.  If the source term were just [[latex($\nabla \phi$)]] then momentum would be conserved but since in each cell there is a factor of [[latex($\rho$)]] that currently enters, it modifies the differenced quantity from cell to cell.  In theory these all cancel but only because [[latex($\nabla^2\phi=\rho$)]]...  In the above treatment - the density used in the source term is derived from the actual potential which allows the source term to be cast as a flux difference.  This allows for the strict momentum conservation.
 
-* Also athena stores combined kinetic/thermal/magnetic with gravitational...  So [[latex($E_{Athena}=E_{Astrobear}+\rho\phi$)]].  This results in [[latex($\dot{E}_{Athena}=\dot{E}_{Astrobear}+\dot{\rho}\phi = \rho \nabla \phi \cdot \mathbf{v}+\phi \nabla \cdot (\rho \mathbf{v}) = \nabla \cdot (\rho \mathbf{v} \phi)$)]] so the energy source term for self gravity can be written as a total divergence.  We need [[latex($\rho v \phi$)]] at cell edges, but [[latex($\rho v$)]] is available from the mass flux and [[latex($\phi$)]] can be averaged using adjacent cells.
+* As as the energy is concerned...  
+Currently the energy added is [[latex($\Delta E_i = \frac{\Delta t}{2\Delta x} \left[(\rho v)_{i-1/2}\Delta x \left(-\frac{\phi_{i}-\phi_{i-1}}{\Delta x} \right) - (\rho v)_{i+1/2}\Delta x \left(+\frac{\phi_{i+1}-\phi_{i}}{\Delta x}\right )\right ]$)]]
+
+Which is equivalent to
+
+[[latex($-\Delta t \left [ \frac{ (\rho v \nabla \phi)_{i+1/2} + (\rho v \nabla \phi)_{i-1/2}}{2} \right ] $)]]
+
+
+If we try to express the energy source term as a total flux we get
+
+[[latex($\dot{E_{kinetic}} = f_{grav} \cdot v = \rho \nabla \phi \cdot v = \rho v \cdot \nabla \phi = \nabla \cdot (\rho v \phi)-\phi \nabla \cdot {\rho v} = \nabla \cdot {\rho v \phi} - \phi \nabla \cdot {v \nabla^2\phi}$)]]
+
+The presence of the v will make it impossible to right this as a total flux.
+
+What if we store the potential in E?  Then [[latex($\dot{E}= \nabla \cdot (\rho v \phi) - \rho \dot{\phi}$)]]  
+
+We now have a flux as well as a source term.  However, the integral of the source term would need to be zero.  So our elliptic equation would have a constraint
+
+
+[[latex($\nabla^2\phi=4 \pi G \rho$)]]
+
+[[latex($\int{\rho \dot{\phi}} = 0$)]]
+
+
+
+
+
+We can conserve total energy before and after the hydro step if we add source terms to the Energy [[latex($\dot{E}=f_{grav} \cdot \mathbf{v} = \rho \nabla \phi \cdot \mathbf{v} = \rho \mathbf{v} \cdot \nabla \phi = \nabla \cdot \left (\rho \mathbf{v} \phi \right) - \phi \nabla \cdot \left ( \rho \mathbf{v} \right ) = \nabla \cdot \left (\rho \mathbf{v} \phi \right) - \phi \dot{\rho}$)]]
+
+
+Consider the difference in gravitational energy...
+
+[[latex($\dot{E_{grav}} = \dot{(\rho\phi)}=\dot{\rho}{\phi}+\rho\dot{\phi}= \nabla \cdot (\rho v) \phi + \rho \dot{\phi} = \nabla \cdot (\rho v \phi) - \rho v \cdot \nabla \phi + \rho \dot{\phi}$)]]
+
+As well as kinetic energy
+
+[[latex($\dot{E_{kinetic}} = f_{grav} \cdot v = \rho \nabla \phi \cdot v = \rho v \cdot \nabla \phi$)]]
+
+For the sum to be zero we need to have [[latex($\nabla \cdot (\rho v \phi) = \rho \dot{\phi}$)]]
+
+The first term [[latex($\dot{\rho}{\phi} = \nabla \cdot (\rho v) \phi = \nabla \cdot (\rho v \phi) - \rho v \cdot \nabla \phi$)]] so we could identify a gravitational energy flux term [[latex($=\rho v \phi$)]]
+
+
+