simplify computation of simplified-SDC reactive source (#1601)

Source/sources/Castro_sources.cpp

@@ -577,7 +577,8 @@ Castro::get_react_source_prim(MultiFab& react_src, Real time, Real dt)
 
     int ng = 0;
 
-    // Carries the contribution of all non-reacting source terms.
+    // Carries the contribution of all non-reacting source terms to
+    // the conserved state (including the hydro update)
 
     MultiFab A(grids, dmap, NUM_STATE, ng);
 
@@ -600,25 +601,6 @@ Castro::get_react_source_prim(MultiFab& react_src, Real time, Real dt)
 
     cons_to_prim(S_noreact, q_noreact, qaux_noreact, time);
 
-    // Compute the primitive version of the old state, q_old
-
-    MultiFab q_old(grids, dmap, NQ, ng);
-    MultiFab qaux_old(grids, dmap, NQAUX, ng);
-
-    cons_to_prim(S_old, q_old, qaux_old, time);
-
-    // Compute the effective advective update on the primitive state.
-    // A(q) = (q* - q_old)/dt
-
-    MultiFab A_prim(grids, dmap, NQ, ng);
-
-    A_prim.setVal(0.0);
-
-    if (dt > 0.0) {
-        MultiFab::Saxpy(A_prim,  1.0 / dt, q_noreact, 0, 0, NQ, ng);
-        MultiFab::Saxpy(A_prim, -1.0 / dt, q_old,     0, 0, NQ, ng);
-    }
-
     // Compute the primitive version of the new state.
 
     MultiFab q_new(grids, dmap, NQ, ng);
@@ -628,14 +610,19 @@ Castro::get_react_source_prim(MultiFab& react_src, Real time, Real dt)
 
     // Compute the reaction source term.
 
+    // I_q = (q^{n+1} - q^n) / dt - A(q)
+    //
+    // but A(q) = (q* - q^n) / dt -- that's the effect of advection w/o burning
+    //
+    // so I_q = (q^{n+1} - q*) / dt
+
     react_src.setVal(0.0, react_src.nGrow());
 
     if (dt > 0.0) {
         MultiFab::Saxpy(react_src,  1.0 / dt, q_new, 0, 0, NQ, ng);
-        MultiFab::Saxpy(react_src, -1.0 / dt, q_old, 0, 0, NQ, ng);
+        MultiFab::Saxpy(react_src, -1.0 / dt, q_noreact, 0, 0, NQ, ng);
     }
 
-    MultiFab::Saxpy(react_src, -1.0, A_prim, 0, 0, NQ, ng);
 
     // Now fill all of the ghost zones.
     Real cur_time = get_state_data(Simplified_SDC_React_Type).curTime();