Inconsistency in documentation of user-defined source terms?

[emprice]

I am implementing a custom source function in Athena++ and am confused by two wiki pages that seem to give conflicting information.

The example here and its surrounding text, "The conservative variables should be updated using these primitive variables and cell-centered magnetic fields.", suggest that prim should be used to update cons in the source function. However, the example here for passive scalars does the opposite, using cons to update cons.

I don't see how both can be correct if cons has not yet been updated to match prim. Is this a matter of taking an implicit vs. explicit step in the source function? If so, maybe that could be stated more explicitly in the text.

[felker]

@c-white looks like you made that edit to the Passive Scalars wiki page in September 2020, FYI

[c-white]

I added some clarifying text to the second page.

The code takes the following steps, starting with synced primitives and conserved quantities:

1. Use primitives to calculate fluxes.
2. Use fluxes to update conserved quantities.
3. Call user source function.
4. Update primitives to match conserved quantities.

For normal source terms (i.e. those that could sensibly be written on the RHS of the equations of motion) whose associated time scale is at least as large as the hydrodynamic step, the correct implementation is for the source function to use the primitives to calculate time derivatives of the conserved quantities, then add these changes into those quantities.

If the timescales are too short, then an implicit step should be taken. This could mean ignoring the primitives and simply updating the given conserved quantities using themselves. Alternatively, it could mean finding how the primitives would implicitly update using themselves, converting both the before and after primitives to conserved quantities, and then adding the difference between these conserved states to the conserved quantities given to the function.

If the source terms are weird (like in the passive scalars example, discontinuously altering values based on thresholds), I hesitate to give any hard rules for best practices. In this particular case the goal was to track failed cells, and the newly updated conserved quantities are the ones who may or may not be violating the condition of interest.

[shilpasarkar30]

This query is in continuation with the query by @emprice and answer by @c-white.
WIKI page: "Therefore, the updates are already reflected in the conservative variables but the primitive variables (prim and prim_scalar) and cell-centered magnetic fields are not updated yet. The conservative variables should be updated using these primitive variables and cell-centered magnetic fields".
I am not able to clearly understand but the second line looks just the reverse of the first. Also, the example given just below the wiki link, supports the second line-->prim used to update cons.

From the answer given by @c-white, I could understand that for hydrodynamical timescales the prim variables should be used to update cons variables. Thus, the cooling example given is for hydrodynamical timescales only?
For MHD cases, we cannot use that cooling function and we need to use conserved variables to update cons variables. Is this correct?

[c-white]

The goal of the function is not to update the primitives. In fact, any change to the primitives here will be discarded. The wiki means to say that the conserved variables have already had flux-divergences applied to them, but they still need any remaining components of the time derivative (i.e., source terms) to be applied. The fact that the flux-divergence update has been applied is important, since it means you have to add source terms into the conserved variables, rather than blindly overwrite them -- those flux divergences will not be added back later on.

---

Notation for below: the primitives are $\rho$, $p$, and $\vec{v}$. The conserved quantities are $d = \rho$, $e = \rho v^2 / 2 + p / (\Gamma - 1)$, and $\vec{m} = \rho \vec{v}$.

The example is indeed only for slow cooling. Ignoring the flux-divergences, the total energy equation it solves is
$$\frac{\mathrm{d}e}{\mathrm{d}t} = -\Lambda(\rho, p, \vec{v}) = \frac{1}{\tau} \max\Big(\frac{\rho}{\Gamma - 1} \Big(\frac{p}{\rho} - T_0\Big), 0\Big).$$
This is being differenced in forward-Euler style, where density and velocity/momentum are not going to change:
$$\frac{e^+ - e^-}{\Delta t} = -\Lambda(\rho, p^-, \vec{v}).$$
The change in conserved variables (LHS) is given by a function of primitives (RHS), and all is well.

Other cooling functions might be too fast for this approach; the safest approach is then to difference them in a backward-Euler way:
$$\frac{e^+ - e^-}{\Delta t} = -\Lambda(\rho, p^+, \vec{v}).$$
This leaves open the question of whether to implement this equation with conserved variables or with the given primitives. In the former case, we would simply write
$$\frac{e^+ - e^-}{\Delta t} = -\Lambda\Big(d, (\Gamma - 1) \Big(e^+ - \frac{m^2}{2 d}\Big), \frac{\vec{m}}{d}\Big),$$
which can be solved for $e^+$ once $\Lambda$ is specified. Alternatively, we could solve
$$\frac{p^+ - p^-}{(\Gamma - 1) \Delta t} = -\Lambda(\rho, p^+, \vec{v})$$
for $p^+$. Then we would update the conserved energy via
$$e^+ = e^- + \frac{p^+ - p^-}{\Gamma - 1}.$$

[shilpasarkar30]

Thank you so much for the clarifications!