Changes between Version 8 and Version 9 of u/johannjc/scratchpad


Ignore:
Timestamp:
10/22/13 15:53:32 (11 years ago)
Author:
Jonathan
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • u/johannjc/scratchpad

    v8 v9  
    1 So Marshak boundary conditions are mixed
     1Momentum 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
    211
    3 [[latex($\left . \left ( E-\frac{2}{3\kappa}\nabla{E} \right ) \right |_0 = \frac{4}{c} F$)]]
     12Normal 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
     20Strang 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
     30Momentum 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
     43SG_x has two non-zero terms...
     44{{{
     45#!latex
     46\[
     47\left [
     48\begin{array}{c}
     490 \\
     50\rho \nabla_x \phi \\
     510 \\
     520 \\
     53\rho v_x \nabla_x \phi
     54\end{array}
     55\right ]
     56\]
     57}}}
     58
     59while SG_PFlux_x(phi) has only 1 non-zero term
     60
     61{{{
     62#!latex
     63\[
     64\left [
     65\begin{array}{c}
     660 \\
     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 \\
     700
     71\end{array}
     72\right ]
     73\]
     74}}}
     75
     76and SG_EFlux(f2x, f2y, f2z, phi) is
     77{{{
     78#!latex
     79\[
     80\left [
     81\begin{array}{c}
     820 \\
     830 \\
     840 \\
     850 \\
     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}}}
    492
    593
    6 One could instead do
    7 
    8 [[latex($E=\frac{4}{c}F$)]]
    9 
    10 or
    11 
    12 [[latex($ \frac{c}{3\kappa} \nabla E = F$)]]
    13 
    14 so that
    15 
    16 [[latex($\frac{\partial E}{\partial t} = \nabla \cdot \left ( \frac{c}{3\kappa} \nabla E \right ) = \nabla \cdot F$)]]
    17 
    18 In any even, if we discretize this mixed equation we get (E_g is ghost zone, E_i is internal zone at left edge)
    19 
    20 [[latex($E_g + E_i - \frac{4}{3\kappa \Delta x} \left ( E_i - E_g \right ) = \frac{8}{c} F$)]]
    21 
    22 which we can solve for [[latex($E_g$)]] as
    23 
    24 [[latex($E_g = \frac{\left ( 1 - \frac{3\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}$)]]
    25 
    26 and we have
    27 [[latex($E_g - E_i = \frac{\left ( -\frac{6\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}$)]]
     94Making the momentum be conserved requires recasting the gravitational force as a total differential
    2895
    2996
    30 If we then plug this into the time evolution equation ignoring portions that are unchanged we get:
     97[[latex($\dot p=-\rho \nabla \phi = -\nabla \cdot F$)]] so  we need to find [[latex($F$)]] such that
    3198
    32 [[latex($E^{n+1}_i-E^{n}_i = \alpha \left ( E_g^*-E_i^* \right )$)]]
     99[[latex($\nabla \cdot F = \rho \nabla \phi$)]]
    33100
    34 where [[latex($\alpha = \frac{c \Delta t}{3 \kappa \Delta x^2}$)]]
    35 we get
     101Since [[latex($\nabla^2 \phi = 4\pi G (\rho-\bar{\rho})$)]]
    36102
    37 [[latex($E^{n+1}_i-E^{n}_i = \alpha \left ( \frac{-\frac{6\kappa \Delta x}{4} E_i^* + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}} \right )$)]]
     103we 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$)]]
     106which 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
     108where 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$)]]
    38109
    39110
    40 where [[latex($E^* = \psi E^{n+1} + \bar{\psi} E^n$)]] is the time averaged energy and depends on the time stepping (crank-nicholson or forward euler)
     111* In more than 1D we have
    41112
    42 
    43 and finally we arrive at
    44 
    45 [[latex($\left ( 1+\frac{2 \psi \alpha}{1+\frac{4}{3\kappa\Delta x}} \right ) E^{n+1}_i = \left ( 1-\frac{2 \bar{\psi} \alpha}{1+\frac{4}{3\kappa\Delta x}} \right )  E^{n}_i+ \frac{\frac{8\alpha}{c} F}{1+\frac{4}{3\kappa \Delta x}}$)]]
    46  
    47 This is very similar to what we normally get   ie
    48 
    49 [[latex($\left ( 1+\psi \alpha \right ) E^{n+1}_i - \psi \alpha E^{n+1}_g = \left ( 1-\bar{\psi} \alpha \right )  E^{n}_i + \bar{\psi} \alpha E^{n}_g$)]]
    50 
    51 
    52 and noting that [[latex($\alpha = \frac{c \Delta t}{3 \kappa \Delta x^2}$)]] for the limiter [[latex($\lambda=1/3$)]]
    53 
    54 or flux is just
    55 
    56 [[latex($F \frac{\Delta t}{\Delta x} \rightarrow F \frac{\Delta t}{\Delta x} \frac{2}{1+\frac{3\kappa\Delta x}{4}}$)]]
    57 
    58 
    59 
    60 Note that [[latex($E^{n+1}_i = E^{n}_i$)]] for
    61 
    62 [[latex($E = \frac{4F}{c}$)]]
    63 
    64 and that if [[latex($3 \kappa \Delta x = 4$)]] we revert to what we would expect if there was no Radiation energy in the ghost zone.
    65 
     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\frac{ \left (\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 \frac{\left ( \delta ^i_j\left(\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}}}