limit artificial viscosity if it breaks species conservation

Source/driver/_cpp_parameters

@@ -39,6 +39,11 @@ allow_non_unit_aspect_zones  bool          0
 # the coefficient of the artificial viscosity
 difmag                       Real          0.1
 
+# when applying the artificial viscosity, do we limit the magnitude of
+# the diffusion coefficient to ensure that the species fluxes all keep
+# the same sign?
+limit_avisc_on_species       bool          1
+
 # the small density cutoff.  Densities below this value will be reset
 small_dens                   Real         -1.e200
 

Source/hydro/advection_util.cpp

@@ -320,6 +320,43 @@ Castro::apply_av(const Box& bx,
 
       div1 = diff_coeff * std::min(0.0_rt, div1);
 
+      if (div1 == 0.0_rt) {
+          return;
+      }
+
+      // we want to prevent the artificial viscosity from breaking
+      // species conservation (sum X_k = 1).  This can happen if the
+      // artificial viscosity flips the sign of some of the fluxes.
+      // Here we scale div1 to ensure that this is not an issue.
+
+      if (limit_avisc_on_species) {
+
+          amrex::Real limit_div1{div1};
+
+          for (int n = 0; n < NumSpec; ++n) {
+
+              Real du{};
+
+              if (idir == 0) {
+                  du = uin(i,j,k,UFS+n) - uin(i-1,j,k,UFS+n);
+              } else if (idir == 1) {
+                  du = uin(i,j,k,UFS+n) - uin(i,j-dg1,k,UFS+n);
+              } else {
+                  du = uin(i,j,k,UFS+n) - uin(i,j,k-dg2,UFS+n);
+              }
+
+              Real fnew = flux(i,j,k,UFS+n) + dx[idir] * div1 * du;
+              if (fnew * flux(i,j,k,UFS+n) < 0) {
+                  limit_div1 = std::min(limit_div1,
+                                        std::abs(flux(i,j,k,UFS+n) / (dx[idir] * du)));
+              }
+          }
+
+          div1 = limit_div1;
+
+      }
+
+
       for (int n = 0; n < NUM_STATE; ++n) {
 
           if (n == UTEMP) {
@@ -384,6 +421,10 @@ Castro::apply_av_rad(const Box& bx,
 
       div1 = diff_coeff * std::min(0.0_rt, div1);
 
+      if (div1 == 0.0_rt) {
+          return;
+      }
+
       for (int n = 0; n < ngroups; ++n) {
 
           Real div_var{};