diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp
index 709b7e6e8e4798d116aa78f9a9eb6979c6b898ae..5daa20d6aa3d340461a949620b57ef69e12173b7 100644
--- a/ProcessLib/ProcessOutput.cpp
+++ b/ProcessLib/ProcessOutput.cpp
@@ -118,75 +118,43 @@ void doProcessOutput(std::string const& file_name,
     }
 
 #ifndef USE_PETSC
-    // the following section is for the output of secondary variables
-
-    auto count_mesh_items = [](MeshLib::Mesh const& mesh,
-                               MeshLib::MeshItemType type) -> std::size_t {
-        switch (type)
-        {
-            case MeshLib::MeshItemType::Cell:
-                return mesh.getNumberOfElements();
-            case MeshLib::MeshItemType::Node:
-                return mesh.getNumberOfNodes();
-            default:
-                break;  // avoid compiler warning
-        }
-        return 0;
-    };
-
-    auto get_or_create_mesh_property = [&mesh, &count_mesh_items](
-        std::string const& property_name, MeshLib::MeshItemType type,
-        unsigned const num_components) {
-        // Get or create a property vector for results.
-        MeshLib::PropertyVector<double>* result = nullptr;
-
-        auto const N = count_mesh_items(mesh, type) * num_components;
-
-        if (mesh.getProperties().hasPropertyVector(property_name))
-        {
-            result = mesh.getProperties().template getPropertyVector<double>(
-                property_name);
-        }
-        else
-        {
-            result =
-                mesh.getProperties().template createNewPropertyVector<double>(
-                    property_name, type, num_components);
-            result->resize(N);
-        }
-        if (!result)
-            OGS_FATAL("Mesh property `%s' could not be created.",
-                      property_name.c_str());
-        if (result->size() != N)
-            OGS_FATAL(
-                "Mesh property `%s' has the wrong size. Actual: %i, expected: "
-                "%i",
-                property_name.c_str(), result->size(), N);
-
-        return result;
-    };
-
     auto add_secondary_var = [&](SecondaryVariable const& var,
                                  std::string const& output_name) {
         {
             DBUG("  secondary variable %s", output_name.c_str());
 
-            auto& result = *get_or_create_mesh_property(
-                output_name, MeshLib::MeshItemType::Node,
+            auto& nodal_values_mesh = *MeshLib::getOrCreateMeshProperty<double>(
+                mesh, output_name, MeshLib::MeshItemType::Node,
                 var.fcts.num_components);
-            assert(result.size() ==
-                   mesh.getNumberOfNodes() * var.fcts.num_components);
+            if (nodal_values_mesh.size() !=
+                mesh.getNumberOfNodes() * var.fcts.num_components)
+            {
+                OGS_FATAL(
+                    "Nodal property `%s' does not have the right number of "
+                    "components. Expected: %d, actual: %d",
+                    output_name.c_str(),
+                    mesh.getNumberOfNodes() * var.fcts.num_components,
+                    nodal_values_mesh.size());
+            }
 
             std::unique_ptr<GlobalVector> result_cache;
             auto const& nodal_values =
                 var.fcts.eval_field(t, x, dof_table, result_cache);
-            assert(nodal_values.size() == result.size());
+            if (nodal_values_mesh.size() !=
+                static_cast<std::size_t>(nodal_values.size()))
+            {
+                OGS_FATAL(
+                    "Secondary variable `%s' did not evaluate to the right "
+                    "number of components. Expected: %d, actual: %d.",
+                    var.name.c_str(), nodal_values_mesh.size(),
+                    nodal_values.size());
+            }
 
             // Copy result
             for (GlobalIndexType i = 0; i < nodal_values.size(); ++i)
             {
                 assert(!std::isnan(nodal_values[i]));
-                result[i] = nodal_values[i];
+                nodal_values_mesh[i] = nodal_values[i];
             }
         }
 
@@ -195,23 +163,38 @@ void doProcessOutput(std::string const& file_name,
             DBUG("  secondary variable %s residual", output_name.c_str());
             auto const& property_name_res = output_name + "_residual";
 
-            auto& result = *get_or_create_mesh_property(
-                property_name_res, MeshLib::MeshItemType::Cell,
+            auto& residuals_mesh = *MeshLib::getOrCreateMeshProperty<double>(
+                mesh, property_name_res, MeshLib::MeshItemType::Cell,
                 var.fcts.num_components);
-            assert(result.size() ==
-                   mesh.getNumberOfElements() * var.fcts.num_components);
+            if (residuals_mesh.size() !=
+                mesh.getNumberOfElements() * var.fcts.num_components)
+            {
+                OGS_FATAL(
+                    "Cell property `%s' does not have the right number of "
+                    "components. Expected: %d, actual: %d",
+                    property_name_res.c_str(),
+                    mesh.getNumberOfElements() * var.fcts.num_components,
+                    residuals_mesh.size());
+            }
 
             std::unique_ptr<GlobalVector> result_cache;
             auto const& residuals =
                 var.fcts.eval_residuals(t, x, dof_table, result_cache);
-
-            assert(residuals.size() == result.size());
+            if (residuals_mesh.size() !=
+                static_cast<std::size_t>(residuals.size()))
+            {
+                OGS_FATAL(
+                    "Thee residual of secondary variable `%s' did not evaluate "
+                    "to the right number of components. Expected: %d, actual: "
+                    "%d.",
+                    var.name.c_str(), residuals_mesh.size(), residuals.size());
+            }
 
             // Copy result
             for (GlobalIndexType i = 0; i < residuals.size(); ++i)
             {
                 assert(!std::isnan(residuals[i]));
-                result[i] = residuals[i];
+                residuals_mesh[i] = residuals[i];
             }
         }
     };
@@ -226,8 +209,6 @@ void doProcessOutput(std::string const& file_name,
         add_secondary_var(secondary_variables.get(external_variable_name),
                           external_variable_name);
     }
-
-    // secondary variables output end
 #else
     (void) secondary_variables;
 #endif // USE_PETSC