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 we have

begin{align*} \partial_t p_i&=\rho \partial_i \phi \\ & = \left (\frac{\partial_j(\partial_j \phi)}{4 \pi G}+\bar{\rho} \right ) \partial_i \phi \\ & = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_j (\partial_i \phi )+\partial_i (\bar{\rho}\phi) \\ & = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_i (\partial_j \phi )+\partial_i (\bar{\rho}\phi) \\ & = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - \frac{1}{2}\partial_i(\partial_j \phi \partial_j \phi)+\partial_i (\bar{\rho}\phi) \\ & = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - \partial_i \left (\frac{\partial_j \phi \partial_j}{2} \phi+\bar{\rho}\phi \right) \\ & = \partial_j \left ( \partial_j \phi \partial_i \phi \right) - \partial_j \left ( \delta ^i_j\left(\frac{\partial_k \phi \partial_k \phi}{2}+\bar{\rho}\phi \right) \right) \\ & = \partial_j T_{ij} \mbox{ where } T_{ij} = \left ( \partial_j \phi \partial_i \phi \right) - \delta ^i_j\left(\frac{\partial_k \phi \partial_k \phi}{2}+\bar{\rho}\phi \right) \\ \end{align*}

• 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.

• As as the energy is concerned…

Currently the energy added is \Delta E_i = \frac{\Delta t}{2\Delta x} \left[(\rho v)_{i-1/2}\Delta x \left(-\frac{\phi_{i}-\phi_{i-1}}{\Delta x} \right) - (\rho v)_{i+1/2}\Delta x \left(+\frac{\phi_{i+1}-\phi_{i}}{\Delta x}\right )\right ]

Which is equivalent to

-\Delta t \left [ \frac{ (\rho v \nabla \phi)_{i+1/2} + (\rho v \nabla \phi)_{i-1/2}}{2} \right ]

If we try to express the energy source term as a total flux we get

\dot{E_{kinetic}} = f_{grav} \cdot v = \rho \nabla \phi \cdot v = \rho v \cdot \nabla \phi = \nabla \cdot (\rho v \phi)-\phi \nabla \cdot {\rho v} = \nabla \cdot {\rho v \phi} - \phi \nabla \cdot {v \nabla^2\phi}

The presence of the v will make it impossible to right this as a total flux.

What if we store the potential in E? Then \dot{E}= \nabla \cdot (\rho v \phi) - \rho \dot{\phi}

We now have a flux as well as a source term. However, the integral of the source term would need to be zero. So our elliptic equation would have a constraint

\nabla^2\phi=4 \pi G \rho

\int{\rho \dot{\phi}} = 0

We can conserve total energy before and after the hydro step if we add source terms to the Energy \dot{E}=f_{grav} \cdot \mathbf{v} = \rho \nabla \phi \cdot \mathbf{v} = \rho \mathbf{v} \cdot \nabla \phi = \nabla \cdot \left (\rho \mathbf{v} \phi \right) - \phi \nabla \cdot \left ( \rho \mathbf{v} \right ) = \nabla \cdot \left (\rho \mathbf{v} \phi \right) - \phi \dot{\rho}

Consider the difference in gravitational energy…

\dot{E_{grav}} = \dot{(\rho\phi)}=\dot{\rho}{\phi}+\rho\dot{\phi}= \nabla \cdot (\rho v) \phi + \rho \dot{\phi} = \nabla \cdot (\rho v \phi) - \rho v \cdot \nabla \phi + \rho \dot{\phi}

As well as kinetic energy

\dot{E_{kinetic}} = f_{grav} \cdot v = \rho \nabla \phi \cdot v = \rho v \cdot \nabla \phi

For the sum to be zero we need to have \nabla \cdot (\rho v \phi) = \rho \dot{\phi}

The first term \dot{\rho}{\phi} = \nabla \cdot (\rho v) \phi = \nabla \cdot (\rho v \phi) - \rho v \cdot \nabla \phi so we could identify a gravitational energy flux term =\rho v \phi

• 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: