have the HLL solver modify the flux passed in (#823)

this fixes a bug where we may have an uninitialized flux

Source/hydro/riemann.cpp

@@ -116,6 +116,13 @@ Castro::cmpflx_plus_godunov(const Box& bx,
           qr_zone[n] = qp(i,j,k,n);
         }
 
+        // pass in the current flux -- the
+        // HLL solver will overwrite this
+        // if necessary
+        for (int n = 0; n < NUM_STATE; n++) {
+          flx_zone[n] = flx(i,j,k,n);
+        }
+
         HLL(ql_zone, qr_zone, cl, cr,
             idir, coord,
             upass_map_p, qpass_map_p,

Source/hydro/riemann_solvers.cpp

@@ -1325,6 +1325,13 @@ Castro::HLL(const Real* ql, const Real* qr,
             const GpuArray<int, npassive>& qpass_map_p,
             Real* flux_hll) {
 
+  // This is the HLLE solver.  We should apply it to zone averages
+  // (not reconstructed states) at an interface in the presence of
+  // shocks to avoid the odd-even decoupling / carbuncle phenomenon.
+  //
+  // See: Einfeldt, B.  et al. 1991, JCP, 92, 273
+  //      Einfeldt, B. 1988, SIAM J NA, 25, 294
+
 
   constexpr Real small_hll = 1.e-10_rt;
 
@@ -1390,6 +1397,8 @@ Castro::HLL(const Real* ql, const Real* qr,
 
   if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) return;
 
+  // we'll overwrite the passed in flux with the HLL flux
+
   bd = 1.0_rt/bd;
 
   // compute the fluxes according to E91, eq. 4.4b -- note that the