Custom source terms in non-barotropic flow

[hwhitehead]

I am attempting to include a custom equation of state that accounts for radiation pressures. As part of this I am making sure that my custom source terms are consistent with a non-barotropic equation of state (previously the flow was isothermal).

I want to check how source terms should update the cons(IEN) value. For instance I have a term that includes a point-source gravity term on the fluid. I have tried reading the src/hydro/srcterms files for reference. As a rough prediction, I expect that if

• $K = \frac{1}{2} \rho v^2$
• $\frac{dK}{dt} = \rho v \cdot \frac{dv}{dt}$

So if I have a source term that accelerates the gas in a cell from $v \rightarrow v' = v + dv$, as well as updating the gas momenta I should also update the total energy such that

• $E \rightarrow E' = E + dE$
• $dE = \rho v \cdot dv$

This seems to be consisted with the methodology in constant_acc.cpp, but the methods in other source files are more complicated, seeming to invoke an averaging of the flux over the cell. For instance in point_mass.cpp or shearing_box.cpp there appears to be an averaging of the density flux over the cell.

Could anyone give me some instruction as to how best to include custom source terms within a non-barotropic equation of state?

[c-white]

Gravity has the property that naive discretizations of the equations tend to not be well balanced, causing noticable errors to accumulate in even simple test problems. With a uniform gravitational field these issues don't appear, but for anything else it often pays to consider discretizations carefully.

As it turns out, a number of Athena developers have published results on this topic. See Jiang et al. 2013 and Mullen et al. 2021. Those papers are concerned with the more complicated case of self gravity, but still have lessons for fixed potentials. In particular, by moving some of the source terms from cell-centered values to face-centered, flux-dependent values, one can exactly conserve total (kinetic + internal + potential) energy.

I believe the code currently implements something like equations (10) and (11) from the earlier paper. If you happen to be in Cartesian coordinates, there is a pervasive issue of even well-balanced discretizations accumulating large errors in curl(g) (which should be identically 0). This is fixed by something like equation (57) in the latter paper.

In summary, your intuition is not wrong, and it might even work for your problem, but there are subtly better methods available.

[tomidakn]

For self-gravity, we (partly) switched to Mullen et al. since 640fdde .

[hwhitehead]

Thanks, I'll have a look at these papers and see if I can find an appropriate method.