for simplified-SDC store only reactive part in plotfile (#1602)

we were storing the entire update, including the advection in the enuc update
AMReX-Astro · Mar 22, 2021 · 5059361 · 5059361
1 parent 994b57a
commit 5059361
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 3 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -1,3 +1,8 @@
+# 21.04
+
+   * For simplified-SDC, we now correctly store only the reactive
+     part of the update for rho_enuc, enuc, and rho_omegadot
+     plotfile variables (#1602)
 
 # 21.02
 

diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp
@@ -435,16 +435,26 @@ Castro::react_state(Real time, Real dt)
 #endif
 
                  if (react_src.contains(i,j,k)) {
+                     // store the reaction data for the plotfile.
+                     // Note, we want to just capture the reaction
+                     // portion here, so we subtract off the advective
+                     // part.
+
+                     // rho omegadot_k
                      for (int n = 0; n < NumSpec; ++n) {
-                         react_src(i,j,k,n) = (U_new(i,j,k,UFS+n) - U_old(i,j,k,UFS+n)) / dt;
+                         react_src(i,j,k,n) = (U_new(i,j,k,UFS+n) - U_old(i,j,k,UFS+n)) / dt - asrc(i,j,k,UFS+n);
                      }
 #if NAUX_NET > 0
+                     // rho auxdot_k
                      for (int n = 0; n < NumAux; ++n) {
-                         react_src(i,j,k,n+NumSpec) = (U_new(i,j,k,UFX+n) - U_old(i,j,k,UFX+n)) / dt;
+                         react_src(i,j,k,n+NumSpec) = (U_new(i,j,k,UFX+n) - U_old(i,j,k,UFX+n)) / dt - asrc(i,j,k,UFX+n);
                      }
 #endif
 
-                     react_src(i,j,k,NumSpec+NumAux) = (U_new(i,j,k,UEINT) - U_old(i,j,k,UEINT)) / dt;
+                     // rho enuc
+                     react_src(i,j,k,NumSpec+NumAux) = (U_new(i,j,k,UEINT) - U_old(i,j,k,UEINT)) / dt - asrc(i,j,k, UEINT);
+
+                     // burn weights
                      react_src(i,j,k,NumSpec+NumAux+1) = amrex::max(1.0_rt, static_cast<Real>(burn_state.n_rhs + 2 * burn_state.n_jac));
                  }