Issue with Derived Variable Outside of (0, 1) Bounds

[zanig]

I'm implementing a custom derived variable in pelelmex_derpseudo to compute a progress variable based on the H2 mass fraction $( Y_{H2} )$. The function uses a $\tanh$-based function $0.5 + \frac{\tanh(K (Y_{H2} - 0.5 Y_{H2,\text{max}}))}{2}$ to map $Y_{H2}$ / $Y_{H2,\text{max}}$ to approximately [0, 1] (with clipping to $$[10^{-8}, 1 - 10^{-8}] $$). However, I occasionally observe values slightly outside [0, 1] in the output when visualizing the solution with paraview or yt, which to my undertanding shouldn’t happen due to the clamping.

Code

Below is my implementation inside of PeleLMeX_DeriveFunc.cpp:

void pelelmex_derpseudo(
    PeleLM* /*a_pelelm*/,
    const Box& bx,
    FArrayBox& derfab,
    int dcomp,
    int ncomp,
    const FArrayBox& statefab,
    const FArrayBox& /*reactfab*/,
    const FArrayBox& /*pressfab*/,
    const Geometry& /*geom*/,
    Real /*time*/,
    const Vector<BCRec>& /*bcrec*/,
    int /*level*/
) {
    AMREX_ASSERT(ncomp == 1);
    auto const in_dat = statefab.array();
    auto const& deri = derfab.array(dcomp);

    //Constants
    const amrex::Real yh2_max = 0.0150375;
    const amrex::Real x0 = yh2_max / 2.0;
    const amrex::Real k = 10.0;
    const amrex::Real clip_min = 1e-8;
    const amrex::Real clip_max = 1.0 - 1e-8;

    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
        //get density
        amrex::Real rhoinv = 1.0 / in_dat(i, j, k, DENSITY);
        //get h2 massfrac
        amrex::Real yh2 = in_dat(i, j, k, FIRSTSPEC + H2_ID) * rhoinv;
        // calc yh2 / yh2_max
        amrex::Real p = yh2 / yh2_max;
        // Compute tanh function: 0.5 + 0.5 * tanh(k * (p - 0.5)) 
        amrex::Real value = 0.5 + 0.5 * std::tanh(k * (p - 0.5));
        // enforce that the value is inside bounds of 0,1
        deri(i, j, k) = amrex::Clamp(value, clip_min, clip_max);
    });
}

• Expected Behavior: The derived variable pseudo should be in $[10^{-8}, 1 - 10^{-8}]$.

• Observed Issue: Some cells at some timesteps have values slightly below 0 or above 1 (e.g., -0.86), see image.

Could this be due to GPU-specific issues such as issues with std::tanh(), or incorrect index access (e.g., DENSITY or H2_ID)?
Is the entry for that derived variable in derfab being modified elsewhere?

Thank you for your help!

PS: I am using MPI and CUDA for compiling.

[baperry2]

Nothing jumps out to me as wrong in your code. I don't see any issue with how you are accessing data. But also it shouldn't matter what data you access or what tanh does, the output of amrex::Clamp (which is device compatible) should always be clamped to the specified limits. See if you can simplify and isolate the issue, e.g., what happens if you just pass a constant value of say 1.5 to Clamp? It may be easier to debug on CPU before moving to GPU.

[zanig]

Thank you for your help!
Ive spotted a minor mistake in my code where i used the "k" both as an index and a constant inside the tanh but that shouldnt really affect the result.

Ive tried your suggestion by running two simulations again with setting value to 1.5 and -0.005 respectively and there the issue is not reproducable.

[esclapez]

I would have blamed Paraview interpolation since the extrema is located in the region of gradient, but you mention testing yt. Have you tried using the AMReX native visualization Amrvis ? It also give you easy access to the actual data table (locally on a user-selected region).
Do you have multiple AMR levels ? If this is not the finest level but maybe the coarsen version of a finer grid, something went wrong there ...