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 F which requires \nabla F=\rho \nabla \phi. Since \nabla^2 \phi = 4\pi G (\rho-\bar{\rho}) we can substitute for \rho and we have \nabla F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi which reduces to \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 F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \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?

\bfseries{\Large Source Terms }}

he original implementation of source terms:

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