diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp
index 050b3151225..85583781104 100644
--- a/singularity-eos/eos/eos_gruneisen.hpp
+++ b/singularity-eos/eos/eos_gruneisen.hpp
@@ -447,10 +447,12 @@ Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const {
   if (slope_at_ref_density < 0) {
     // Something is very wrong in the construction of this EOS... just error out
     using PortsOfCall::printf;
+#ifndef NDEBUG
     printf("ERROR: The pressure is decreasing as density increases at the reference\n"
            "       density, %.15g, for temperature, %.15g. Check that the reference\n"
            "       temperature is set correctly. This is an unstable state.",
            _rho0, temperature);
+#endif // NDEBUG
     PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface");
   }
 
@@ -498,11 +500,13 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature(
     auto pres_max = PressureFromDensityTemperature(rho_upper, temp);
     if (press > pres_max) {
       // We're off the EOS surface
+#ifndef NDEBUG
       using PortsOfCall::printf;
       printf("ERROR: Requested pressure, %.15g, exceeds maximum, %.15g, for \n"
              "       temperature, %.15g.\n"
              "       setting pressure to maximum allowed value.\n",
              press, pres_max, temp);
+#endif // NDEBUG
       p_used = pres_max;
     }
     // Construct a reasonable guess for the density
diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp
index 69d472410f6..7a85b31f1a0 100644
--- a/singularity-eos/eos/eos_helmholtz.hpp
+++ b/singularity-eos/eos/eos_helmholtz.hpp
@@ -903,9 +903,11 @@ PORTABLE_INLINE_FUNCTION Real Helmholtz::lTFromRhoSie_(const Real rho, const Rea
           math_utils::pow10(electrons_.lTMax()), HELM_EOS_EPS, T, nullptr,
           options_.VERBOSE, false);
       if (status != RootFinding1D::Status::SUCCESS) {
+#ifndef NDEBUG
         if (options_.VERBOSE) {
           printf("Newton-Raphson failed to converge, falling back to regula falsi\n");
         }
+#endif // NDEBUG
         status = RootFinding1D::regula_falsi(
             [&](Real T) {
               Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
@@ -944,15 +946,19 @@ PORTABLE_INLINE_FUNCTION Real Helmholtz::lTFromRhoSie_(const Real rho, const Rea
   lT = std::log10(T);
   // Make sure the result is within the table bounds
   if (lT < electrons_.lTMin()) {
+#ifndef NDEBUG
     if (options_.VERBOSE) {
       printf("Temperature below table limit, setting lT = lTMin. (lT = %f)\n", lT);
     }
+#endif // NDEBUG
     lT = electrons_.lTMin();
   }
   if (lT > electrons_.lTMax()) {
+#ifndef NDEBUG
     if (options_.VERBOSE) {
       printf("Temperature above table limit, setting lT = lTMax. (lT = %f)\n", lT);
     }
+#endif // NDEBUG
     lT = electrons_.lTMax();
   }
   IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp
index ba6dc20d36e..b2c79b3cf83 100644
--- a/singularity-eos/eos/eos_mgusup.hpp
+++ b/singularity-eos/eos/eos_mgusup.hpp
@@ -209,8 +209,10 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const
   Real eta = 1.0 - robust::ratio(_rho0, rho);
   Real f1 = 1.0 - _s * eta;
   if (f1 <= 0.0) {
+#ifndef NDEBUG
     printf("f1, eta, rho, rho0, s = %.14e %.14e %.14e %.14e %.14e\n", f1, eta, rho, _rho0,
            _s);
+#endif // NDEBUG
     PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho "
                                    "give a negative argument for a logarithm.");
   }
@@ -241,7 +243,6 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const
       PORTABLE_WARN("Hugoniot Temperature not converged");
 #endif // NDEBUG
     }
-    // printf("sum=%e\n",sum);
     temp = _T0 - pf * ((G0os - 3.0) + exp(-G0os) * (G0os * G0os - 4.0 * G0os + 2.0) *
                                           (std::log(f1) + sum));
     temp = temp * exp(_G0 * eta);
diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp
index 084b6a840b9..0cfac330817 100644
--- a/singularity-eos/eos/get_sg_eos_functors.hpp
+++ b/singularity-eos/eos/get_sg_eos_functors.hpp
@@ -133,6 +133,7 @@ struct init_functor {
     any_bad_vals = any_bad_vals || error_utils::bad_value(spvol_v(i), "spvol");
 
     if (any_bad_vals) {
+#ifndef NDEBUG
       using PortsOfCall::printf;
       printf("### Bad Value Output state:\n");
       printf("  ~~Bulk state~~\n");
@@ -148,6 +149,7 @@ struct init_functor {
         printf(" %24.15g", frac_mass_v(i, m));
         printf("\n");
       }
+#endif // NDEBUG
       PORTABLE_ALWAYS_ABORT(
           "Bad values INPUT to singularity-eos interface. See output for details");
     }
diff --git a/singularity-utils/singularity-utils/root-finding-1d/root_finding.hpp b/singularity-utils/singularity-utils/root-finding-1d/root_finding.hpp
index 1e052e102fb..4163b9e83f0 100644
--- a/singularity-utils/singularity-utils/root-finding-1d/root_finding.hpp
+++ b/singularity-utils/singularity-utils/root-finding-1d/root_finding.hpp
@@ -143,10 +143,12 @@ PORTABLE_INLINE_FUNCTION bool set_bracket(const T &f, Real &a, const Real guess,
     }
   }
   // if we get here then we failed to bracket a root
+#ifndef NDEBUG
   if (verbose) {
     printf("set_bracket failed to bound a root! %.14e %.14e %.14e %.14e %.14e %.14e\n", a,
            guess, b, ya, yg, yb);
   }
+#endif // NDEBUG
   return false;
 }
 
@@ -175,9 +177,11 @@ PORTABLE_INLINE_FUNCTION Status regula_falsi(const T &f, const Real ytarget,
     } else {
       // ya, yg, and yb have the same sign
       if (!set_bracket(func, a, guess, b, ya, yg, yb, verbose)) {
+#ifndef NDEBUG
         if (verbose) {
           printf("regula_falsi failed! %.14e %.14e %.14e %.14e\n", ytarget, guess, a, b);
         }
+#endif // NDEBUG
         return Status::FAIL;
       }
     }
@@ -222,10 +226,12 @@ PORTABLE_INLINE_FUNCTION Status regula_falsi(const T &f, const Real ytarget,
   }
   auto status = Status::SUCCESS;
   if (iteration_count == max_iter) {
+#ifndef NDEBUG
     if (verbose) {
       printf("root finding reached the maximum number of iterations.  likely not "
              "converged\n");
     }
+#endif // NDEBUG
     status = Status::FAIL;
   }
   if (counts != nullptr) {
@@ -279,10 +285,12 @@ PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget,
     // Per default, we fail if the root is out of bounds controlled by
     // fail_on_bound_root.
     if ((_x <= a && _xold <= a) || (_x >= b && _xold >= b)) {
+#ifndef NDEBUG
       if (verbose) {
         printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess,
                a, b);
       }
+#endif // NDEBUG
       if (fail_on_bound_root) {
         status = Status::FAIL;
       }
@@ -291,10 +299,12 @@ PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget,
     _x = std::max(std::min(_x, b), a);
   }
   if (iter >= max_iter) {
+#ifndef NDEBUG
     if (verbose) {
       printf("root finding reached the maximum number of iterations.  likely not "
              "converged\n");
     }
+#endif // NDEBUG
     status = Status::FAIL;
   }