How to record the change due to UserExplicitSourceFunction properly?

[c-kj]

Background

Hi all. I'm dealing with a simulation with black hole accretion & cooling. I implemented the BH sink and cooling via the source term function and enroll them with UserExplicitSourceFunction .

Now I want to record the change due to my source functions, e.g. the total accreted mass, the work done by BH gravity and the total cooling loss, so that I can analyze the mass & energy conservation. To accomplish this, I record the value changes of cons onto a pmb->ruser_meshblock_data in the source term functions, and output them in hst output file.

The problem, theoretically

However, I realize that due to the multi-stage hydro integrator, my recorded value is always larger than the "actual" value. To figure this out, I read the Athena++ paper (Stone et.al. 2020) and I carefully check the source code of TimeIntegratorTaskList::IntegrateHydro. From eqn (11) in (Stone et.al. 2020), for each stage $s$, (here I add the change due to source term $S_s$ based on my understanding):

$$ \begin{align} u^{(1)} \leftarrow & u^{(1)} + \delta_s u^{(0)} \\ u^{(0)} \leftarrow & \gamma_{s0}u^{(0)} + \gamma_{s1}u^{(1)} + \beta_{s,s-1}\Delta t F(u^{(0)}) + \textcolor{blue}{S_s} \end{align} $$

the source term added at the previous stage will come into this stage with a coefficient $\gamma_{s0}$. For the integrator "rk2" I'm using, the actual change due to source term in each cycle is

$$ \Delta u = \gamma_{10} S_0+S_1, ~~~ (\gamma_{10}=1/2 \text{ for rk2}) $$

while the change I recorded is $S_0 + S_1$ , so recorded value is 1.5x of the actual value.

Tests I've done

I've tested this problem with periodic boundary conditions, so that the conservation is guaranteed, thus I could calculate the actual change of total mass/energy with the hst file. This difference appears in both the accreted mass and the total cooling loss, and the ratio between the "recorded values" and "actual values" is ~ 1.5, just as predicted above (small fluctuation due to the change between stages).

Current workarounds and their shortcoming

To avoid this, I have currently tried 2 methods:

1. When accumulating the values in the source term functions, multiply a extra factor $\gamma_{s0}$
2. Only let the BH sink & cooling source term work in the last stage of each cycle

I tested both methods, both giving the correct results. However, these methods suffer from some problems:

1. Method 1 need to get the ptlist->stage_wghts[stage-1].gamma1, however, stage_wghts is a private menber of TimeIntegratorTaskList. Or, I have to manually choose the weight $\gamma_{s0}$ for each integrator, which is cumbersome.
2. Method 2 skips all but the last stage in each cycle, the accuracy of the source term might be lower. This may be ok for BH sink, but it's unreasonable to let BH gravity and cooling only works in the last stage.

What's more, both methods need to get the information of the current stage, so I have to modify the main.cpp to get this value from the for loop.

The Question:

I'm wondering if there's any way to properly record the source terms' contribution, without affecting the accuracy and performance of the simulation?

I'd welcome all the discussions, and if anyone have some elegant trick, I'd be very appreciated.

[tomidakn]

I do understand what you want to do. We do not have an official way to do it, so I'd say anything is "proper" as long as it works. In this case, I'd recommend Method 1 because Method 2 is less accurate as you pointed out. The question is then how to do it with minimal modification to the main code. Honestly, I think it is OK for you to change stage_wghts to public or make Mesh a friend class of TimeIntegratorTaskList so that you can obtain gamma during the initialization. But if you want to do everything on the user side (i.e. within your pgen file), perhaps you need to read the integrator from ParameterInput and store gamma accordingly in Mesh::InitUserMeshData for example. I agree this is cumbersome, but you only need to write it once.

[c-kj]

Thanks Tomida san, it's nice to hear the developer's opinion, I'm more confident to modify the source code now : )
Following your suggestion, I add a getter method to TimeIntegratorTaskList, and use some trick to calculate the coefficient for each stage. Now the recorded mass and energy exactly follow conservation! Cheers.

[c-kj]

I'd like to share my solution to this problem, in case someone else may face the same issue.

The source term function in each stage has a different contribution to the final result in this cycle, but we couldn't get this "weight" easily, since it involves some iterative process of registers phydro->u and phydro->u1. Except the need to modify the source code, we also need to compute these weights, and I come up with the following trick:

  Real GetSourceTermWeight(int stage) {
    auto stage_wghts = ptlist->GetStageWeight();
    Real u = 0, u1 = 0;  
    for (int s = 1; s <= nstages; ++s) { 
      auto stage_wght = stage_wghts[s-1];
      if (stage_wght.main_stage) {  
        u1 = u1 + stage_wght.delta * u;
        u  = stage_wght.gamma_1 * u + stage_wght.gamma_2 * u1 + (s==stage ? 1 : 0);  
      }
    }
    return u;
  }

This piece of code imitate the procedure in TimeIntegratorTaskList::IntegrateHydro, but only inject a "source term" at the target stage, so we can extract the "weight" we want.

To use this, store these weights when enter the main loop, then in the source term functions, when accumulating the changes, multiply by this weight for the current stage, then the result will be correct.

Note: I tested this trick for vl2, rk2, rk3, rk4 integrator, but not sure if it's correct for "ssprk5_4" or orbital advection.