Context Navigation

← Previous Post
Next Post →

Momentum Conserving Self Gravity

Making the momentum be conserved requires recasting the gravitational force as a total differential

\dot p=-\rho \nabla \phi = -\nabla \cdot F so we need to find F such that

\nabla \cdot F = \rho \nabla \phi

Since \nabla^2 \phi = 4\pi G (\rho-\bar{\rho})

we can substitute for \rho and we have

\nabla \cdot F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi which is equivalent (in 1D) to \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 F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi

In more than 1D we have

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 \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_j (\partial_i \phi )}{4 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_i (\partial_j \phi )}{4 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_i\frac{(\partial_j \phi \partial_j \phi)}{8 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_i\left ( \frac{ \partial_j \phi \partial_j \phi}{8 \pi G} +\bar{\rho}\phi \right) \\ & = \partial_j \frac{\left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_j \left ( \delta^i_j \left ( \frac{\partial_k \phi \partial_k \phi}{8 \pi G} + \bar{\rho}\phi \right) \right ) \\ & = \partial_j T_{ij} \mbox{ where } T_{ij} = \frac{\left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \delta^i_j\left(\frac{\partial_k \phi \partial_k \phi}{8 \pi G}+\bar{\rho}\phi \right) \\ \end{align*}

or we can write it as \mathbf{T}=\frac{1}{4\pi G}\left [ \nabla \phi \nabla \phi - \frac{1}{2} \left ( \nabla \phi \right ) \cdot \left ( \nabla \phi \right ) \mathbf{I} \right ]

Momentum Conserving Self Gravity Source Terms in AstroBEAR

q - initial state
qLx - left edge of x - interface in predictor
qLy - lower edge of y - interface in predictor
qRy - upper edge of y - interface in predictor
q2 - intermediate cell centered quantities
q2Lx - updated left edge of x - interface used for final riemann solve
fx - predictor x flux
f2x - final x flux
q3 - new values

Normal CTU:

qLx = Reconstruct(q), etc… (for all 6 qxx states)
fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
q2 = q+fx+fy+fz
q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
q3=q+f2x+f2y+f2z

Strang split self-gravity:

q = q + sg(q,phi)*hdt
qLx = Reconstruct(q), etc… (for all 6 qxx states)
fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
q2 = q+fx+fy+fz
q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
q3=q+f2x+f2y+f2z
q3=q3+sq(q3,phi)*hdt

where

sg=sg_x+sg_y+sg+z and sg_x has two non-zero terms…

[ \left [ \begin{array}{c} 0 \\ -\rho \nabla_x \phi \\ 0 \\ 0 \\ -\rho v_x \nabla_x \phi \end{array} \right ] \]

Momentum conserving self-gravity

qLx = Reconstruct(q), etc… (for all 6 qxx states)
qLx(px) =qLx(px)+sg_x(q,phi)*hdt
qLy(py) = qLy(py)+sg_y(q,phi)*hdt
qLy(pz) = qLy(pz)+sg_z(q,phi)*hdt
fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
q2 = q+fx+fy+fz + sg(q,phi)
q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
q2Lx=q2Lx+sg_y(q,phi)+sg_z(q,phi) (for all 6 q2xx states)
f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
f2x(p)=f2x(p)+SG_PFlux_x(phi) (same for y and z)
q3=q+f2x+f2y+f2z+SG_E(f2x(rho),f2y(rho), f2z(rho), phi)

The final updates for momentum include true 'fluxes'

SG_PFlux_x(phi) has only 1 non-zero term

[ \left [ \begin{array}{c} 0 \\ \frac{1}{8G\pi} \left ( (\nabla_x \phi)^2 - (\nabla_y\phi)^2 - (\nabla_z\phi)^2 \right ) + \bar{\rho}\phi \ \\ \frac{1}{4G\pi} \nabla_x \phi \nabla_y\phi \\ \frac{1}{4G\pi} \nabla_x \phi \nabla_z\phi \\ 0 \end{array} \right ] \]

and an energy update term that is derived from mass fluxes.

SG_E(f2x, f2y, f2z, phi)

[ \left [ \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ -\rho v_x \nabla_x \phi - \rho v_y \nabla_y \phi - \rho v_z \nabla_z \phi \end{array} \right ] \]

For a description of an algorithm that also conserves energy see http://adsabs.harvard.edu/abs/2013NewA...19...48J

Posted: 11 years ago (Updated: 11 years ago)
Author: Jonathan
Categories: (none)

Attachments (2)

Screen Shot 2013-10-22 at 4.35.47 PM.png (26.7 KB) - added by Jonathan 11 years ago.
Screen Shot 2013-10-22 at 4.35.24 PM.png (25.4 KB) - added by Jonathan 11 years ago.

astrobear

Context Navigation

Momentum Conserving Self Gravity

Momentum Conserving Self Gravity Source Terms in AstroBEAR

Attachments (2)

Comments

Download in other formats: