strang splitting does not work?

[BenWibking]

It looks like the gravitational potential is hard-coded to be only a function of radius:

athenapk/src/hydro/srcterms/gravitational_field.hpp

Line 52 in ca22b73

const Real g_r = gravitationalField.g_from_r(r);

In general, it should be a function of all of the coordinates. (For my simulations, it should be a function of height.)

[pgrete]

Does it make sense for your usecase to work with GravitationalField gravitationalField or is your is it a (very simple) function?
In the former case, I could imagine to work with if constexpr inside the function to differentiate between different kinds of potentials.
In the latter case, it may even be easiest/cleanest to just add a problem specific, self-contained source function.

[BenWibking]

It's a simple function. Should it go in PostStepMeshUserWorkInLoop?

[pgrete]

I suggest to hook into the Hydro::ProblemSourceUnsplit or Hydro::ProblemSourceFirstOrder function (cf., the turbulence problem generator).
The "unplit" version is called during every stage of the integrator (i..e, the beta_dt should be used for $\Delta t$) whereas the "first order" one is only called once (so the full $\Delta t$ for the entire cycle should be used).

[BenWibking]

Is there a Strang-split version of this? Otherwise I'll just use Hydro::ProblemSourceFirstOrder.

[pgrete]

Is there a Strang-split version of this? Otherwise I'll just use Hydro::ProblemSourceFirstOrder.

Yes, it's called ProblemSourceStrangSplit 😄

Also this applies: https://github.com/parthenon-hpc-lab/athenapk/blob/main/src/hydro/hydro_driver.cpp#L160

[BenWibking]

Well, that's convenient. So I need to read from the cons variables and then update both cons and prim?

[BenWibking]

What is the reason for this restriction ? Could this be simplified by adding another conserved -> primitive conversion step in the task list?

[pgrete]

Well, that's convenient. So I need to read from the cons variables and then update both cons and prim?

yes (that's the "safe" option)

What is the reason for this restriction ? Could this be simplified by adding another conserved -> primitive conversion step in the task list?

The current implementation is mostly motivated by performance concerns, which doesn't mean it has to stay that way.
So far the source terms I've been using update the data they were reading, so I could easily update both prim and cons in those source terms without needing to read and write all data again.
For example, for cooling the update to the pressure is trivial based on required calculations in the kernel whereas a separate ConsToPrim would need to read all data (density, momentum, energy, and potentially magnetic fields) again, just to update the pressure variable.

[BenWibking]

Ok, that makes sense. I agree about there being a performance cost, but I'd expect it to be very small - have you measured this?

[pgrete]

Ok, that makes sense. I agree about there being a performance cost, but I'd expect it to be very small - have you measured this?

Yes, (but it also strongly depends on the total amount of physics included).
IIRC it's double digit percent for MHD because that kernel is extremely memory bound as there virtual no compute for cons to prim conversion.

[BenWibking]

I am getting some exceptionally weird (wrong) results when using Strang splitting for gravity:

Working version (first order splitting):

Bad version (Strang splitting):

Do you have any idea what could be going wrong here?

[pgrete]

That does indeed look strange.
If you could provide a reproducer (i.e., the code version/commit to use an a sample input file), I'll have a look.