NaN values encountered when trying to setup user-defined auxiliary variables

[RemiLehe]

First of all, thanks a lot for developing Castro ; the code is very useful to me!

I recently tried to use auxiliary variables (described in this page of the documentation, accessed with UFX, etc.). More specifically, I did a small test where I started from the 1D Sedov example and added an auxiliary variable that is supposed to mirror exactly eint_e(as a way to check that the variable behaves as expected): see changes in the commit below
RemiLehe@ccfd30f

However, when running this example, I get an error message:

amrex::Abort::0::State has NaNs in the rho_myenerg component::check_for_nan() !!!

Is this expected? Is there a way to work around this?

How to easily reproduce this issue

Clone my fork from Castro and switch to the branch test_aux_variables (this branch is identical to the current main branch of Castro, but with the additional commit listed above). Then run the 1d Sedov cylindrical test.

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

[yut23]

Thanks for the easy reproduction instructions, they help a lot!

My first step for debugging is always to build with DEBUG=TRUE and run with amrex.fpe_trap_{invalid,zero,overflow}=1. This adds runtime bounds checks, fills uninitialized arrays with NaNs, and makes any floating-point exceptions immediately trigger a crash. With those, the backtrace points to https://github.com/AMReX-Astro/Microphysics/blob/533f712d2d217ea4f4063f3a809eb16a839c9ef6/EOS/eos_composition.H#L73, which looks like division by zero. The problem is that your single species has Z=0, so sum is always zero.

If I set Z=1.0 in species.net, it makes it through the initialization. However, it still crashes at https://github.com/RemiLehe/Castro/blob/ccfd30fb9d6f236c493bc771305b8edc1fe7caea/Exec/hydro_tests/Sedov/problem_source.H#L36 with the error message amrex::Abort::0:: (-1,0,0,1) is out of bound (0:31,0:0,0:0,0:8) !!!, so it looks like state doesn't have any ghost cells. I'm not sure where this should be addressed myself, but one of the other developers should be able to help.

[RemiLehe]

Great, thanks a lot!
Let me know if there is any additional test I can to do to help.

[WeiqunZhang]

The issue is Castro::do_new_sources will fail if the state multifabs it gets does not have ghost cells and problem source function needs ghost cells. If you do the following changes, the test will run. But It's probably not the right way to fix it, because of the extra communications even when it's not needed. Maybe we only need to do that when UFX exists?

@zingale

diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
index a3ceec8f8..8df985d3f 100644
--- a/Source/driver/Castro.cpp
+++ b/Source/driver/Castro.cpp
@@ -3045,11 +3045,21 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep)
 
             if (getLevel(lev).apply_sources()) {
                 bool apply_sources_to_state = true;
+                MultiFab Sborder_old(S_old.boxArray(), S_old.DistributionMap(),
+                                     S_old.nComp(), 1);
+                MultiFab Sborder_new(S_new.boxArray(), S_new.DistributionMap(),
+                                     S_new.nComp(), 1);
+                getLevel(lev).expand_state(Sborder_old,
+                                           getLevel(lev).state[State_Type].prevTime(),
+                                           Sborder_old.nGrow());
+                getLevel(lev).expand_state(Sborder_new,
+                                           getLevel(lev).state[State_Type].curTime(),
+                                           Sborder_new.nGrow());
                 getLevel(lev).do_new_sources(
 #ifdef MHD
                                 Bx_new, By_new, Bz_new,
 #endif
-                                source, S_old, S_new, time, dt_advance_local, apply_sources_to_state);
+                                source, Sborder_old, Sborder_new, time, dt_advance_local, apply_sources_to_state);
             }
 
             if (use_retry && dt_advance_local < dt_amr && getLevel(lev).keep_prev_state) {
diff --git a/Source/sources/Castro_sources.cpp b/Source/sources/Castro_sources.cpp
index 1c4981a7d..68e0bbfea 100644
--- a/Source/sources/Castro_sources.cpp
+++ b/Source/sources/Castro_sources.cpp
@@ -266,7 +266,10 @@ Castro::do_new_sources (Real time, Real dt)
 {
     advance_status status {};
 
-    MultiFab& S_new = get_new_data(State_Type);
+    MultiFab Sborder_new(Sborder.boxArray(), Sborder.DistributionMap(), Sborder.nComp(), Sborder.nGrowVect());
+    expand_state(Sborder_new, time, Sborder_new.nGrow()); // time is the current time
+    // Sborder has the old state data with ghost cells.
+    // Sborder_new has the new state data with ghost cells.
 
     MultiFab& new_source = get_new_data(Source_Type);
 
@@ -280,7 +283,7 @@ Castro::do_new_sources (Real time, Real dt)
 #ifdef MHD
                    Bx_new, By_new, Bz_new,
 #endif
-                   new_source, Sborder, S_new, time, dt);
+                   new_source, Sborder, Sborder_new, time, dt);
 
     return status;
 }