Changes between Version 11 and Version 12 of u/johannjc/scratchpad


Ignore:
Timestamp:
01/08/14 17:03:01 (11 years ago)
Author:
Jonathan
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • u/johannjc/scratchpad

    v11 v12  
    1 Momentum Conserving Self Gravity Source Terms in AstroBEAR
    2 * q - initial state
    3 * qLx - left edge of x - interface in predictor
    4 * qLy - lower edge of y - interface in predictor
    5 * qRy - upper edge of y - interface in predictor
    6 * q2 - intermediate cell centered quantities
    7 * q2Lx - updated left edge of x - interface used for final riemann solve
    8 * fx - predictor x flux
    9 * f2x - final x flux
    10 * q3 - new values
     1== Self-gravitating equilibrium profile using density profile and external pressure ==
    112
    12 Normal CTU:
    13 * qLx = Reconstruct(q), etc... (for all 6 qxx states)
    14 * fx=Riemann_Solve(qLx,qRx), etc... (for all 3 edges)
    15 * q2 = q+fx+fy+fz
    16 * q2Lx=qLx+fy+fz etc... (for all 6 q2xx states)
    17 * f2x=Riemann_Solve(q2Lx,q2Rx)  (for all 3 edges)
    18 * q3=q+f2x+f2y+f2z
    19 
    20 Strang split self-gravity:
    21 * '''q = q + SG(q, phi)*hdt'''
    22 * qLx = Reconstruct(q), etc... (for all 6 qxx states)
    23 * fx=Riemann_Solve(qLx,qRx), etc... (for all 3 edges)
    24 * q2 = q+fx+fy+fz
    25 * q2Lx=qLx+fy+fz etc... (for all 6 q2xx states)
    26 * f2x=Riemann_Solve(q2Lx,q2Rx)  (for all 3 edges)
    27 * q3=q+f2x+f2y+f2z
    28 * '''q3=q3+SG(q3, phi)*hdt'''
    29 
    30 Momentum conserving self-gravity
    31 * qLx = Reconstruct(q), etc... (for all 6 qxx states)
    32 * ''' qLx(px) =qLx(px)+SG_x(q,phi)*hdt'''
    33 * ''' qLy(py) = qLy(py)+SG_y(q,phi)*hdt'''
    34 * ''' qLy(pz) = qLy(pz)+SG_z(q,phi)*hdt'''
    35 * fx=Riemann_Solve(qLx,qRx), etc... (for all 3 edges)
    36 * q2 = q+fx+fy+fz '''+ SG(q,phi)'''
    37 * q2Lx=qLx+fy+fz etc... (for all 6 q2xx states)
    38 * '''q2Lx=q2Lx+SG_y(q,phi)+SG_z(q,phi)''' (for all 6 q2xx states)
    39 * f2x=Riemann_Solve(q2Lx,q2Rx)  (for all 3 edges)
    40 * '''f2x(p)=f2x(p)+SG_PFlux_x(phi) ''' (same for y and z)
    41 * q3=q+f2x+f2y+f2z'''+SG_EFlux_(f2x(rho),f2y(rho), f2z(rho), phi)'''
    42 
    43 SG_x has two non-zero terms...
    44 {{{
    45 #!latex
    46 \[
    47 \left [
    48 \begin{array}{c}
    49 0 \\
    50 \rho \nabla_x \phi \\
    51 0 \\
    52 0 \\
    53 \rho v_x \nabla_x \phi
    54 \end{array}
    55 \right ]
    56 \]
    57 }}}
    58 
    59 while SG_PFlux_x(phi) has only 1 non-zero term
    60 
    61 {{{
    62 #!latex
    63 \[
    64 \left [
    65 \begin{array}{c}
    66 0 \\
    67 \frac{1}{8G\pi} \left ( (\nabla_x \phi)^2 - (\nabla_y\phi)^2 - (\nabla_z\phi)^2 \right ) + \bar{\rho}\phi \ \\
    68 \frac{1}{4G\pi} \nabla_x \phi \nabla_y\phi \\
    69 \frac{1}{4G\pi} \nabla_x \phi \nabla_z\phi \\
    70 0
    71 \end{array}
    72 \right ]
    73 \]
    74 }}}
    75 
    76 and SG_EFlux(f2x, f2y, f2z, phi) is
    77 {{{
    78 #!latex
    79 \[
    80 \left [
    81 \begin{array}{c}
    82 0 \\
    83 0 \\
    84 0 \\
    85 0 \\
    86 \rho v_x \nabla_x \phi + \rho v_y \nabla_y \phi + \rho v_z \nabla_z \phi
    87 \end{array}
    88 \right ]
    89 \]
    90 
    91 }}}
     3Mathematically, we are solving the equation of hydrostatic equilibrium [[latex($\nabla P(r) = -\rho(r) \nabla \phi(r)$)]] subject to the constraint that  [[latex($P(R)=P_0$)]].  Since [[latex($\phi(r) = \frac{G M(r)}{r}$)]] we can rewrite the equation for HSE as [[latex($\frac{dP}{dr}=\frac{GM(r)\rho(r)}{r^2}$)]] and we can express the enclosed mass as [[latex($M(r)=4\pi\displaystyle\int_0^r{\rho(r') r'^2dr'}$)]]
    924
    935
    94 Making the momentum be conserved requires recasting the gravitational force as a total differential
    95 
    96 
    97 [[latex($\dot p=-\rho \nabla \phi = -\nabla \cdot F$)]] so  we need to find [[latex($F$)]] such that
    98 
    99 [[latex($\nabla \cdot F = \rho \nabla \phi$)]]
    100 
    101 Since [[latex($\nabla^2 \phi = 4\pi G (\rho-\bar{\rho})$)]]
    102 
    103 we can substitute for [[latex($\rho$)]] and we have
    104 
    105 [[latex($\nabla \cdot F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi$)]]
    106 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 ]$)]]
    107 
    108 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$)]]
    109 
    110 
    111 * In more than 1D we have
    112 
    113 {{{
    114 #!latex
    115 \begin{align*}
    116 \partial_t p_i&=\rho \partial_i \phi \\
    117 & = \left (\frac{\partial_j(\partial_j \phi)}{4 \pi G}+\bar{\rho} \right ) \partial_i \phi \\
    118 & = \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) \\
    119 & = \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) \\
    120 & = \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) \\
    121 & = \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) \\
    122 & = \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 )  \\
    123 & = \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) \\
    124 \end{align*}
    125 }}}