User-defined auxiliary variable does not evolve as expected

[RemiLehe]

First of all, thanks again for your efforts in developing Castro and making it available on Github!

Overview

I am currently trying to use the auxiliary variables in Castro.
In order to make sure that I am using them properly, I am running a small test where I added an auxiliary variable that is supposed to mirror exactly eint_e, for the 1D cylindrical Sedov case. The full commit is here, but the key changes are:

• in species.net:

# Test auxiliary variable: create a variable that mirrors the internal energy
__aux_myenergy

• in problem_initialize_state_data.H:

state(i,j,k,UFX) = eint; // test variable that mirrors the internal energy

• in problem_source.H:

    Real rho_e = state(i,j,k,UFX);
    Real p = 2./3 * rho_e;
    src(i,j,k,UFX) = -p * divU;

(See the full commit to see the calculation of divU, which I copied from Castro/Source/sources/Castro_thermo.cpp)

However, when running the corresponding Sedov test, eint_e and this new auxiliary variable (my_energy) do not seem to evolve in the exact same way:

I assume my calculation in problem_source.H is too naive, but could you point to where it should be improved, e.g. I am using the wrong state? or should I explicitly call the eos? or something else?

How to reproduce this issue

Clone my fork from Castro and switch to the branch avoid_outofbounds_error (this branch is identical to the current main branch of Castro, but with the additional auxiliary variable mentioned above, as well as the change suggested by @WeiqunZhang to avoid out-of-bounds errors in auxiliary variables: #2982). Then run the 1d Sedov cylindrical test.

git clone https://github.com/RemiLehe/Castro
cd Castro
git checkout avoid_outofbounds_error
git submodule update --init
cd Exec/hydro_tests/Sedov
make DIM=1 -j 4
./Castro1d.gnu.MPI.ex inputs.1d.cyl

The above figure can then be obtained with the following Python code:

import yt
import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(6,6))

# Loop over different time of the simulation
for i, iteration in enumerate([0, 60, 120]):
    
    plt.subplot(311+i)
    ds = yt.load('./plt%05d/' %iteration)

    # Extract data
    level = 3
    ad0 = ds.covering_grid(level=level, left_edge=ds.domain_left_edge, dims=[ds.domain_dimensions[0]*2**level, 1, 1])
    r = np.linspace( ds.domain_left_edge[0], ds.domain_right_edge[0], ds.domain_dimensions[0]*2**level)

    # Plot the different curves
    X = ad0['eint_e'].to_ndarray()
    plt.plot(r, X.squeeze(), label='eint_e')
    X = ad0['myenerg'].to_ndarray()
    plt.plot(r, X.squeeze(), label='myenergy')

    plt.title(f't = %.1e' %ds.current_time.to_value())

    plt.legend(loc=0)
    plt.xlim(0, 0.3)
    plt.xlabel('Radius r')
    plt.tight_layout()

[RemiLehe]

Update: in the last commit the above-mentioned branch avoid_outofbounds_error, I now explicitly call the EOS with:

    eos_rep_t eos_state;
    eos_state.rho = state(i,j,k,URHO);
    eos_state.e = state(i,j,k,UFX) / state(i,j,k,URHO); // internal energy per unit volume
    for (int n = 0; n < NumSpec; n++) {
        eos_state.xn[n] = state(i,j,k,UFS+n)/state(i,j,k,URHO);
    }
#if NAUX_NET > 0
    for (int n = 0; n < NumAux; n++) {
        eos_state.aux[n] = state(i,j,k,UFX+n)/state(i,j,k,URHO);
    }
#endif
    eos(eos_input_re, eos_state);

But I still get the same result (i.e. eint_e and myenergy do not evolve in the same way.

[RemiLehe]

Hm, I think that this is related to @zingale's comment here:
#2679 (comment)

We want to do this by using the Godunov state to construct p h ∇ ⋅ U -- that way we will be consistent with ρ E and ρ e . To enable that we need to add a hook in consup_hydro that allows a problem to add additional term to the conserved update.