diff --git a/kernel/include/Operators/OperatorBase.tpp b/kernel/include/Operators/OperatorBase.tpp
index 45a4450d..ee7a391d 100644
--- a/kernel/include/Operators/OperatorBase.tpp
+++ b/kernel/include/Operators/OperatorBase.tpp
@@ -205,6 +205,7 @@ OperatorBase<T, DIM>::OperatorBase(const std::vector<std::string>& integrators,
     : mfem::Operator(this->compute_total_height(spatials), this->compute_total_width(spatials)),
       params_(default_params_),
       RHS(NULL),
+      rhs_solver_(nullptr),
       current_dt_(0.0),
       current_time_(0.0),
       height_(height),
@@ -241,6 +242,7 @@ OperatorBase<T, DIM>::OperatorBase(const std::vector<std::string>& integrators,
     : mfem::Operator(this->compute_total_height(spatials), this->compute_total_width(spatials)),
       params_(default_params_),
       RHS(NULL),
+      rhs_solver_(nullptr),
       current_dt_(0.0),
       current_time_(0.0),
       height_(height),
diff --git a/kernel/include/Operators/TransientOperator.hpp b/kernel/include/Operators/TransientOperator.hpp
index 9876a31c..48481277 100644
--- a/kernel/include/Operators/TransientOperator.hpp
+++ b/kernel/include/Operators/TransientOperator.hpp
@@ -84,9 +84,8 @@ class TransientOperator : public OperatorBase<T, DIM>, public mfem::TimeDependen
   /// Left-Hand-Side
   mfem::ParBlockNonlinearForm* LHS;
 
-  // CCI
-  mfem::HypreParMatrix* Mmat;
-  // CCI
+  std::vector<mfem::HypreParMatrix*> Mmat_;
+
   void build_mass_matrix(const std::vector<mfem::Vector>& u_vect);
   void build_lhs_nonlinear_form(const double dt, const std::vector<mfem::Vector>& u);
   bool constant_mass_matrix_{true};
diff --git a/kernel/include/Operators/TransientOperator.tpp b/kernel/include/Operators/TransientOperator.tpp
index 797bb3b5..d1754389 100644
--- a/kernel/include/Operators/TransientOperator.tpp
+++ b/kernel/include/Operators/TransientOperator.tpp
@@ -283,7 +283,10 @@ void TransientOperator<T, DIM>::build_mass_matrix(const std::vector<mfem::Vector
   }
 
   this->M_solver_.clear();
-
+  for (auto* mat : this->Mmat_) {
+    delete mat;
+  }
+  this->Mmat_.clear();
   for (unsigned int i = 0; i < u_vect.size(); i++) {
     if (M != nullptr) {
       delete M;
@@ -306,12 +309,13 @@ void TransientOperator<T, DIM>::build_mass_matrix(const std::vector<mfem::Vector
     M->Assemble(0);
     M->Finalize(0);
 
-    Mmat = M->ParallelAssemble();
-    std::unique_ptr<mfem::HypreParMatrix> Me(Mmat->EliminateRowsCols(this->ess_tdof_list_[i]));
+    this->Mmat_.emplace_back(M->ParallelAssemble());
+    std::unique_ptr<mfem::HypreParMatrix> Me(
+        this->Mmat_[i]->EliminateRowsCols(this->ess_tdof_list_[i]));
 
     this->mass_matrix_solver_ =
         std::make_shared<LSolver>(this->mass_solver_, this->mass_solver_params_,
-                                  this->mass_precond_, this->mass_precond_params_, *Mmat);
+                                  this->mass_precond_, this->mass_precond_params_, *this->Mmat_[i]);
     this->M_solver_.emplace_back(this->mass_matrix_solver_->get_solver());
   }
 }
@@ -700,6 +704,11 @@ void TransientOperator<T, DIM>::free_memory() {
   this->LHS = nullptr;
   delete this->reduced_oper;
   this->reduced_oper = nullptr;
+
+  for (auto* mat : this->Mmat_) {
+    delete mat;
+  }
+  this->Mmat_.clear();
 }
 
 /**
diff --git a/kernel/include/PostProcessing/postprocessing.tpp b/kernel/include/PostProcessing/postprocessing.tpp
index bef34eaf..424eb665 100644
--- a/kernel/include/PostProcessing/postprocessing.tpp
+++ b/kernel/include/PostProcessing/postprocessing.tpp
@@ -87,7 +87,7 @@ void PostProcessing<T, DC, DIM>::get_parameters() {
   this->enable_compute_energies_ =
       this->params_.template get_param_value_or_default<bool>("enable_compute_energies", true);
   this->enable_save_specialized_at_iter_ = this->params_.template get_param_value_or_default<bool>(
-      "enable_save_specialized_at_iter", false);
+      "enable_save_specialized_at_iter", true);
   this->force_clean_output_dir_ =
       this->params_.template get_param_value_or_default<bool>("force_clean_output_dir", false);
   if (this->params_.has_parameter("iso_val_to_compute")) {