| 3 | | 1D CTU |
| 4 | | Reconstruct left and right states (wr,wl) |
| 5 | | Add gradphi to normal velocity components at interface states |
| 6 | | wr(i)+=phi(i)-phi(i-1) |
| 7 | | wl(i)+=phi(i)-phi(i-1) |
| 8 | | Add other source terms using the reconstructed left and right states wl+=s(wl) |
| 9 | | wr+=s(wr) |
| 10 | | get fluxes |
| 11 | | calculate cell centered time centered values of density and pressure using fluxes |
| 12 | | update cell centered time centered values with source terms just as for interface states (using original cell centered values for density etc) |
| | 3 | == 1D CTU == |
| | 4 | * Reconstruct left and right states (wr,wl) |
| | 5 | * Add gradphi to normal velocity components at interface states |
| | 6 | * wr(i)+=phi(i)-phi(i-1) |
| | 7 | * wl(i)+=phi(i)-phi(i-1) |
| | 8 | * Add other source terms using the reconstructed left and right states wl+=s(wl) |
| | 9 | * wr+=s(wr) |
| | 10 | * get fluxes |
| | 11 | * calculate cell centered time centered values of density, momentum, and energy using original values and fluxes |
| | 12 | * Update time centered momentum with gravo-source terms using original density and phi (no energy source term) pdV work done in fluxes??? - not sure what is going on here??? and calculate the time centered cell centered pressure |
| | 13 | |
| | 14 | * For a static potential update the momentum using the time centered density and the energy using the mass fluxes times grad phi at the cell edges |
| | 15 | * For self gravity do something slightly more complicated...... |
| | 16 | [[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 F$)]] which requires [[latex($\nabla 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 F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi$)]] which reduces to [[latex($\nabla F=\nabla \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 as [[latex($F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi$)]] |
| | 17 | * 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. |
| | 18 | |
| | 19 | * 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. |
| | 20 | |
| | 21 | |
| | 22 | * Add cooling terms using time centered density and pressure |
| | 23 | * Add fluxes and store fluxes at edges as well as mass flux everywhere |
| | 24 | |
| | 25 | * After calculating new [[latex($\phi$)]] modify momentum and energy by effectively using time centered [[latex($\phi$)]]. This allows for second order without storing [[latex($\dot{\phi}$)]]. |
| | 26 | |
| | 27 | I'm still puzzled about the time centered cell centered pressure being calculated without subtracting off the gravitational energy - and if [[latex($\rho \phi$)]] is subtracted off earlier - then why no potential heating? Seems like this would underestimate the cooling strength... Possible bug? |
| | 28 | |
| | 29 | |
| | 30 | |
