diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp
index 72be2fcc..9e37af02 100644
--- a/src/reactingFlow.cpp
+++ b/src/reactingFlow.cpp
@@ -1608,7 +1608,7 @@ void ReactingFlow::step() {
       // update wdot quantities at full substep in Yn/Tn state
       updateMixture();
       updateThermoP();
-      updateDensity(0.0);
+      updateDensity(0.0, false);
       speciesProduction();
       heatOfFormation();
 
@@ -2503,7 +2503,7 @@ void ReactingFlow::evaluatePlasmaConductivityGF() {
   sigma_gf_.SetFromTrueDofs(sigma_);
 }
 
-void ReactingFlow::updateDensity(double tStep) {
+void ReactingFlow::updateDensity(double tStep, bool update_mass_matrix) {
   Array<int> empty;
   Rmix_gf_.GetTrueDofs(tmpR0a_);
 
@@ -2541,9 +2541,11 @@ void ReactingFlow::updateDensity(double tStep) {
   }
   rn_gf_.SetFromTrueDofs(rn_);
 
-  MsRho_form_->Update();
-  MsRho_form_->Assemble();
-  MsRho_form_->FormSystemMatrix(empty, MsRho_);
+  if (update_mass_matrix) {
+    MsRho_form_->Update();
+    MsRho_form_->Assemble();
+    MsRho_form_->FormSystemMatrix(empty, MsRho_);
+  }
 
   // project to p-space in case not same as vel-temp
   R0PM0_gf_.SetFromTrueDofs(rn_);
diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp
index 7500c2f6..1387ffce 100644
--- a/src/reactingFlow.hpp
+++ b/src/reactingFlow.hpp
@@ -366,7 +366,7 @@ class ReactingFlow : public ThermoChemModelBase {
   void updateMixture();
   void updateThermoP();
   void extrapolateState();
-  void updateDensity(double tStep);
+  void updateDensity(double tStep, bool update_mass_matrix = true);
   void updateBC(int current_step);
   void updateDiffusivity();
   void computeSystemMass();