Source terms in AstroBEAR. All source terms are currently handled in a strang split way. That is - each hydro advance step is bracketed by two source steps of a half time step. Out of curiosity I browsed the athena source code to see how their CTU implementation handles source terms (since I couldn't find it in their CTU paper)

1D CTU

• Reconstruct left and right states (wr,wl)
• Add gradphi to normal velocity components at interface states
  • wr(i)+=phi(i)-phi(i-1)
  • wl(i)+=phi(i)-phi(i-1)
• Add other source terms using the reconstructed left and right states
  • wl+=s(wl)
  • wr+=s(wr)
• get fluxes
• calculate cell centered time centered values of density, momentum, and energy using original values and fluxes
• 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

• 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
• For self gravity do something slightly more complicated……

\Delta p=-\frac{\Delta t}{\Delta x} \left(F_{i+1/2}-F_{i-1/2} \right ) which is just \dot{p} = f_{grav} where f_{grav}=-\nabla \cdot F which requires \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, \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

• But casting it this way gives you strict momentum conservation. If the source term were just \nabla \phi then momentum would be conserved but since in each cell there is a factor of \rho that currently enters, it modifies the differenced quantity from cell to cell. In theory these all cancel but only because \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 E_{Athena}=E_{Astrobear}+\rho\phi. This results in \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 \rho v \phi at cell edges, but \rho v is available from the mass flux and \phi can be averaged using adjacent cells.

• Add cooling terms using time centered density and pressure
• Add fluxes and store fluxes at edges as well as mass flux everywhere

• After calculating new \phi modify momentum and energy by effectively using time centered \phi. This allows for second order without storing \dot{\phi}.

I'm still puzzled about the time centered cell centered pressure being calculated without subtracting off the gravitational energy - and if \rho \phi is subtracted off earlier - then why no potential heating? Seems like this would underestimate the cooling strength… Possible bug?

2D CTU

• First 4 steps are identical as in 1D but now for both directinos
• Next calculate cell centered emf using initial values - and integrate to corner (upwind or otherwise) and update the interface magnetic field using CT
• Update interface states with transverse flux gradients
• Add source terms from transverse flux gradients
• Add source terms for self gravity arising from x2-flux gradient ie at left interface state…
  • v^y_{i+1/2,j} -= \frac{\Delta t}{2} \frac{\phi_{i,j+1}-\phi_{i,j-1}}{2 \Delta x}
  • E_{i+1/2,j} -= \frac{\Delta t}{2} \left [ ( \rho v )_{i,j+1/2} \frac{\phi_{i,j+1}-\phi_{i,j}}{2\Delta x} + ( \rho v )_{i,j-1/2} \frac{\phi_{i,j}-\phi_{i,j-1}}{2\Delta x} \right ]

• calculate cell centered time centered values of density, momentum, and energy using original values and fluxes
• 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

• Calculate new fluxes from corrected interface states

• 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

\bfseries{\Large Source Terms }}

he original implementation of source terms:

he combined 2nd--4th/5th order RK scheme: