diff --git a/CMakeLists.txt b/CMakeLists.txt
index c8523e8f2078c43769bb82afc9e406639a4dec97..afc39bb68171c04e64c70010878655b13e48aca1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -182,7 +182,10 @@ if(OGS_USE_MPI)
 endif()
 
 # Eigen
-add_definitions(-DEIGEN_INITIALIZE_MATRICES_BY_ZERO) # TODO check if needed
+if(CMAKE_BUILD_TYPE STREQUAL "Debug")
+    add_definitions(-DEIGEN_INITIALIZE_MATRICES_BY_NAN)
+endif()
+
 if (EIGEN_NO_DEBUG)
     add_definitions(-DEIGEN_NO_DEBUG)
 endif()
diff --git a/MaterialLib/FractureModels/MohrCoulomb.cpp b/MaterialLib/FractureModels/MohrCoulomb.cpp
index 8a403be9230eb131fd7d39da4a79bc521299c8f0..dd68dfc9d618997450fe4d49362b3f58801e8756 100644
--- a/MaterialLib/FractureModels/MohrCoulomb.cpp
+++ b/MaterialLib/FractureModels/MohrCoulomb.cpp
@@ -66,8 +66,7 @@ void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
     MaterialPropertyValues const mat(_mp, t, x);
     Eigen::VectorXd const dw = w - w_prev;
 
-    Eigen::MatrixXd Ke(2,2);
-    Ke.setZero();
+    Eigen::MatrixXd Ke = Eigen::MatrixXd::Zero(2, 2);
     Ke(0,0) = mat.Ks;
     Ke(1,1) = mat.Kn;
 
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
index 8ece8e1dd01c3a3d0138e0f1b16fc3d5fe0474d3..11bb4d6f9a747028a8db10ac12a37bafa019ca52 100644
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ b/MaterialLib/SolidModels/Lubby2-impl.h
@@ -29,8 +29,9 @@ calculatedGdEBurgers()
 {
     Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
                   Lubby2<DisplacementDim>::KelvinVectorSize>
-        dGdE;
-    dGdE.setZero();
+        dGdE =
+            Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
+                          Lubby2<DisplacementDim>::KelvinVectorSize>::Zero();
     dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
                                 Lubby2<DisplacementDim>::KelvinVectorSize>()
         .diagonal()
diff --git a/MaterialLib/SolidModels/Lubby2.h b/MaterialLib/SolidModels/Lubby2.h
index 5c3516b7bfefc11294547b4b49228c91f312f83b..c58be6405a17bc4075787efe813fb1880b55e76f 100644
--- a/MaterialLib/SolidModels/Lubby2.h
+++ b/MaterialLib/SolidModels/Lubby2.h
@@ -109,10 +109,13 @@ public:
     {
         MaterialStateVariables()
         {
+            // Previous time step values are not initialized and are set later.
             eps_K_t.resize(KelvinVectorSize);
-            eps_K_j.resize(KelvinVectorSize);
             eps_M_t.resize(KelvinVectorSize);
-            eps_M_j.resize(KelvinVectorSize);
+
+            // Initialize current time step values
+            eps_K_j.setZero(KelvinVectorSize);
+            eps_M_j.setZero(KelvinVectorSize);
         }
 
         MaterialStateVariables& operator=(MaterialStateVariables const&) =
diff --git a/MeshLib/IO/VtkIO/VtuInterface-impl.h b/MeshLib/IO/VtkIO/VtuInterface-impl.h
index 6f8c904802dec858792066b320d6cb7e2da34f42..44da8a37f7c51de5f993c11e48c8cae958158a44 100644
--- a/MeshLib/IO/VtkIO/VtuInterface-impl.h
+++ b/MeshLib/IO/VtkIO/VtuInterface-impl.h
@@ -72,6 +72,10 @@ bool VtuInterface::writeVTU(std::string const& file_name,
     vtuWriter->SetNumberOfPieces(num_partitions);
     vtuWriter->SetStartPiece(rank);
     vtuWriter->SetEndPiece(rank);
+#else
+    // avoid unused parameter warnings.
+    (void)num_partitions;
+    (void)rank;
 #endif
 
     return (vtuWriter->Write() > 0);
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 7ec961e1075c62f4f4814782265aee571a736d25..298efe7c713322c328371ffe00f9c44ffb510e33 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -163,12 +163,12 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
             //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__specific_body_force}
             config.getConfigParameter<std::vector<double>>(
                 "specific_body_force");
-        if (specific_body_force.size() != DisplacementDim)
+        if (b.size() != DisplacementDim)
             OGS_FATAL(
                 "The size of the specific body force vector does not match the "
                 "displacement dimension. Vector size is %d, displacement "
                 "dimension is %d",
-                specific_body_force.size(), DisplacementDim);
+                b.size(), DisplacementDim);
 
         std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index 31797fd333dc8989a0fb9b436ae8cbf5cb4eacb8..e6a28f4f191ab9fa961f53601c8bc80a5e74bbf3 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -216,10 +216,13 @@ public:
                 _integration_method.getWeightedPoint(ip).getWeight() *
                 sm_u.integralMeasure * sm_u.detJ;
 
-            ip_data.sigma_eff.resize(kelvin_vector_size);
-            ip_data.sigma_eff_prev.resize(kelvin_vector_size);
-            ip_data.eps.resize(kelvin_vector_size);
+            // Initialize current time step values
+            ip_data.sigma_eff.setZero(kelvin_vector_size);
+            ip_data.eps.setZero(kelvin_vector_size);
+
+            // Previous time step values are not initialized and are set later.
             ip_data.eps_prev.resize(kelvin_vector_size);
+            ip_data.sigma_eff_prev.resize(kelvin_vector_size);
 
             ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
                 DisplacementDim, displacement_size>::Zero(DisplacementDim,
@@ -291,16 +294,19 @@ public:
                 displacement_size + pressure_size>>(
             local_rhs_data, displacement_size + pressure_size);
 
-        typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-        laplace_p.setZero(pressure_size, pressure_size);
+        typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+            ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                             pressure_size);
 
-        typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-        storage_p.setZero(pressure_size, pressure_size);
+        typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+            ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                             pressure_size);
 
         typename ShapeMatricesTypeDisplacement::template MatrixType<
             displacement_size, pressure_size>
-            Kup;
-        Kup.setZero(displacement_size, pressure_size);
+            Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+                displacement_size, pressure_size>::Zero(displacement_size,
+                                                        pressure_size);
 
         double const& dt = _process_data.dt;
 
diff --git a/ProcessLib/LIE/Common/Utils.cpp b/ProcessLib/LIE/Common/Utils.cpp
index 905c26cf8f472f6d31defa3a081b4a6b5b2b967f..b345f028d4b20370181381506360ed61387b6292 100644
--- a/ProcessLib/LIE/Common/Utils.cpp
+++ b/ProcessLib/LIE/Common/Utils.cpp
@@ -28,6 +28,7 @@ void computeNormalVector(MeshLib::Element const& e, unsigned const global_dim,
         auto v1 = (*e.getNode(1)) - (*e.getNode(0));
         element_normal[0] = -v1[1];
         element_normal[1] = v1[0];
+        element_normal[2] = 0;  // not used in 2d but needed for normalization
         element_normal.normalize();
     }
     else if (global_dim == 3)
diff --git a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 972eec043d4b8d9ca18e7db2e57dc483373a1064..618f338480c2fd1f70e963d048db691cd1b3823e 100644
--- a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -185,12 +185,12 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
             //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__specific_body_force}
             config.getConfigParameter<std::vector<double>>(
                 "specific_body_force");
-        if (specific_body_force.size() != GlobalDim)
+        if (b.size() != GlobalDim)
             OGS_FATAL(
                 "The size of the specific body force vector does not match the "
                 "displacement dimension. Vector size is %d, displacement "
                 "dimension is %d",
-                specific_body_force.size(), GlobalDim);
+                b.size(), GlobalDim);
 
         std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index b7a6c4c60ad02f774a957dd919ddcff425eea137..7001de80f5e4e4aa2b53cb278d1f1755a9b4d6b7 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -77,8 +77,8 @@ HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
             sm_u.detJ * sm_u.integralMeasure *
             integration_method.getWeightedPoint(ip).getWeight();
 
-        ip_data.H_u.resize(GlobalDim,
-                           ShapeFunctionDisplacement::NPOINTS * GlobalDim);
+        ip_data.H_u.setZero(GlobalDim,
+                            ShapeFunctionDisplacement::NPOINTS * GlobalDim);
         computeHMatrix<
             GlobalDim, ShapeFunctionDisplacement::NPOINTS,
             typename ShapeMatricesTypeDisplacement::NodalRowVectorType, HMatrixType>(
@@ -86,11 +86,16 @@ HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
         ip_data.N_p = sm_p.N;
         ip_data.dNdx_p = sm_p.dNdx;
 
-        ip_data.w.resize(GlobalDim);
+        // Initialize current time step values
+        ip_data.w.setZero(GlobalDim);
+        ip_data.sigma_eff.setZero(GlobalDim);
+
+        // Previous time step values are not initialized and are set later.
         ip_data.w_prev.resize(GlobalDim);
-        ip_data.sigma_eff.resize(GlobalDim);
         ip_data.sigma_eff_prev.resize(GlobalDim);
+
         ip_data.C.resize(GlobalDim, GlobalDim);
+
         ip_data.aperture0 = (*frac_prop.aperture0)(0, x_position)[0];
         ip_data.aperture = ip_data.aperture0;
 
@@ -163,23 +168,24 @@ assembleBlockMatricesWithJacobian(
     // in a displacement vector
     auto const index_normal = GlobalDim - 1;
 
-    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-    laplace_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-    storage_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
     typename ShapeMatricesTypeDisplacement::template MatrixType<
         displacement_size, pressure_size>
-        Kgp;
-    Kgp.setZero(displacement_size, pressure_size);
+        Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
+            displacement_size, pressure_size>::Zero(displacement_size,
+                                                    pressure_size);
 
     using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
     using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
-    GlobalDimMatrix local_k_tensor;
-    local_k_tensor.setZero();
-    GlobalDimMatrix local_dk_db_tensor;
-    local_dk_db_tensor.setZero();
+    GlobalDimMatrix local_k_tensor = GlobalDimMatrix::Zero();
+    GlobalDimMatrix local_dk_db_tensor = GlobalDimMatrix::Zero();
 
     auto const& gravity_vec = _process_data.specific_body_force;
 
@@ -309,8 +315,7 @@ computeSecondaryVariableConcreteWithVector(
 
     using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
     using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
-    GlobalDimMatrix local_k_tensor;
-    local_k_tensor.setZero();
+    GlobalDimMatrix local_k_tensor = GlobalDimMatrix::Zero();
 
     auto const& gravity_vec = _process_data.specific_body_force;
 
@@ -364,10 +369,8 @@ computeSecondaryVariableConcreteWithVector(
 
     double ele_b = 0;
     double ele_k = 0;
-    Eigen::Vector2d ele_sigma_eff;
-    ele_sigma_eff.setZero();
-    Eigen::Vector2d ele_w;
-    ele_w.setZero();
+    Eigen::Vector2d ele_sigma_eff = Eigen::Vector2d::Zero();
+    Eigen::Vector2d ele_w = Eigen::Vector2d::Zero();
     double ele_Fs = - std::numeric_limits<double>::max();
     for (auto const& ip : _ip_data)
     {
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index e2cdc3eca41e577fa220cafc993eea637ccd55b7..fe74d623f9df55e6c78e8fa569049a928b86a727 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -83,9 +83,10 @@ HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
         ip_data.N_p = sm_p.N;
         ip_data.dNdx_p = sm_p.dNdx;
 
-        ip_data.sigma_eff.resize(kelvin_vector_size);
+        ip_data.sigma_eff.setZero(kelvin_vector_size);
+        ip_data.eps.setZero(kelvin_vector_size);
+
         ip_data.sigma_eff_prev.resize(kelvin_vector_size);
-        ip_data.eps.resize(kelvin_vector_size);
         ip_data.eps_prev.resize(kelvin_vector_size);
         ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
 
@@ -163,16 +164,19 @@ assembleBlockMatricesWithJacobian(
 {
     assert (this->_element.getDimension() == GlobalDim);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
-    laplace_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
-    typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
-    storage_p.setZero(pressure_size, pressure_size);
+    typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+                                                         pressure_size);
 
     typename ShapeMatricesTypeDisplacement::template MatrixType<
         displacement_size, pressure_size>
-        Kup;
-    Kup.setZero(displacement_size, pressure_size);
+        Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+            displacement_size, pressure_size>::Zero(displacement_size,
+                                                    pressure_size);
 
     double const& dt = _process_data.dt;
     auto const& gravity_vec = _process_data.specific_body_force;
@@ -363,12 +367,9 @@ computeSecondaryVariableConcreteWithBlockVectors(
         q.noalias() = - k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
     }
 
-    Eigen::Vector3d ele_stress;
-    ele_stress.setZero();
-    Eigen::Vector3d ele_strain;
-    ele_strain.setZero();
-    Eigen::Vector3d ele_velocity;
-    ele_velocity.setZero();
+    Eigen::Vector3d ele_stress = Eigen::Vector3d::Zero();
+    Eigen::Vector3d ele_strain = Eigen::Vector3d::Zero();
+    Eigen::Vector3d ele_velocity = Eigen::Vector3d::Zero();
 
     for (auto const& ip_data : _ip_data)
     {
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
index 50ce1cbfae3df8d5284c47bdca05219b4f03353e..f59c43b4774b3986999cf39a0c27941b2c530bf1 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
@@ -72,20 +72,24 @@ SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
         auto& ip_data = _ip_data[ip];
         ip_data._detJ = sm.detJ;
         ip_data._integralMeasure = sm.integralMeasure;
-        ip_data._h_matrices.resize(
-            DisplacementDim,
-            ShapeFunction::NPOINTS * DisplacementDim);
+        ip_data._h_matrices.setZero(DisplacementDim,
+                                    ShapeFunction::NPOINTS * DisplacementDim);
 
         computeHMatrix<
             DisplacementDim, ShapeFunction::NPOINTS,
             typename ShapeMatricesType::NodalRowVectorType, HMatrixType>(
             sm.N, ip_data._h_matrices);
 
-        ip_data._w.resize(DisplacementDim);
-        ip_data._w_prev.resize(DisplacementDim);
-        ip_data._sigma.resize(DisplacementDim);
+        // Initialize current time step values
+        ip_data._w.setZero(DisplacementDim);
+        ip_data._sigma.setZero(DisplacementDim);
+
+        // Previous time step values are not initialized and are set later.
         ip_data._sigma_prev.resize(DisplacementDim);
+        ip_data._w_prev.resize(DisplacementDim);
+
         ip_data._C.resize(DisplacementDim, DisplacementDim);
+
         ip_data._aperture0 = (*_fracture_property->aperture0)(0, x_position)[0];
         ip_data._aperture_prev = ip_data._aperture0;
 
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
index 3c9c657b3e20585a4ae094b00f53ee81d0e6b652..0b7e4a902055b407e9304bfeb84b5a7db6cdd4a3 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
@@ -70,12 +70,17 @@ SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
         ip_data.dNdx = sm.dNdx;
         ip_data._detJ = sm.detJ;
         ip_data._integralMeasure = sm.integralMeasure;
-        ip_data._sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
+
+        // Initialize current time step values
+        ip_data._sigma.setZero(KelvinVectorDimensions<DisplacementDim>::value);
+        ip_data._eps.setZero(KelvinVectorDimensions<DisplacementDim>::value);
+
+        // Previous time step values are not initialized and are set later.
         ip_data._sigma_prev.resize(
             KelvinVectorDimensions<DisplacementDim>::value);
-        ip_data._eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
         ip_data._eps_prev.resize(
             KelvinVectorDimensions<DisplacementDim>::value);
+
         ip_data._C.resize(KelvinVectorDimensions<DisplacementDim>::value,
                           KelvinVectorDimensions<DisplacementDim>::value);
 
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index 3ff08e674ecc3a95714cacebc57f9d35399ef5ec..978f8a1e5d73a3b515bfbad7a425e060114a7592 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -89,12 +89,16 @@ SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
         ip_data._detJ = sm.detJ;
         ip_data._integralMeasure = sm.integralMeasure;
 
-        ip_data._sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
+        // Initialize current time step values
+        ip_data._sigma.setZero(KelvinVectorDimensions<DisplacementDim>::value);
+        ip_data._eps.setZero(KelvinVectorDimensions<DisplacementDim>::value);
+
+        // Previous time step values are not initialized and are set later.
         ip_data._sigma_prev.resize(
             KelvinVectorDimensions<DisplacementDim>::value);
-        ip_data._eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
         ip_data._eps_prev.resize(
             KelvinVectorDimensions<DisplacementDim>::value);
+
         ip_data._C.resize(KelvinVectorDimensions<DisplacementDim>::value,
                           KelvinVectorDimensions<DisplacementDim>::value);
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 58e9e755aae2be0e2098af565aaa544f4e6c0c4e..139c1e29b99a244e0a5f5539aaf57a5bf2f0e266 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -184,11 +184,16 @@ public:
             ip_data.N = sm.N;
             ip_data.dNdx = sm.dNdx;
 
-            ip_data.sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
+            // Initialize current time step values
+            ip_data.sigma.setZero(
+                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data.eps.setZero(KelvinVectorDimensions<DisplacementDim>::value);
+
+            // Previous time step values are not initialized and are set later.
             ip_data.sigma_prev.resize(
                 KelvinVectorDimensions<DisplacementDim>::value);
-            ip_data.eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
-            ip_data.eps_prev.resize(KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data.eps_prev.resize(
+                KelvinVectorDimensions<DisplacementDim>::value);
 
             _secondary_data.N[ip] = shape_matrices[ip].N;
         }
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
index 00549b3d140bf6f38c2adf9b64571182c57ed227..8fbf7a92b69e2f9c6df5b6351120b16db8076184 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -100,8 +100,8 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
     auto Met = local_M.template block<temperature_size, temperature_size>(
         temperature_matrix_index, temperature_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index 63403dc4a6ed7bf6eb155214b941c432d81cb7f7..a35ab72aefc3f165ba576efe76ab51a61b786a32 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -155,12 +155,12 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
             //! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__specific_body_force}
             config.getConfigParameter<std::vector<double>>(
                 "specific_body_force");
-        if (specific_body_force.size() != DisplacementDim)
+        if (b.size() != DisplacementDim)
             OGS_FATAL(
                 "The size of the specific body force vector does not match the "
                 "displacement dimension. Vector size is %d, displacement "
                 "dimension is %d",
-                specific_body_force.size(), DisplacementDim);
+                b.size(), DisplacementDim);
 
         std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index 0b5d8737ba7e2490193249407fef2a0cdd9621b6..6433b367edd5afacc8bf637e064282b41ff88505 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -210,11 +210,11 @@ public:
                 _integration_method.getWeightedPoint(ip).getWeight() *
                 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
 
-            ip_data.sigma.resize(kelvin_vector_size);
-            ip_data.sigma_prev.resize(kelvin_vector_size);
-            ip_data.eps.resize(kelvin_vector_size);
-            ip_data.eps_m.resize(kelvin_vector_size);
-            ip_data.eps_m_prev.resize(kelvin_vector_size);
+            ip_data.sigma.setZero(kelvin_vector_size);
+            ip_data.sigma_prev.setZero(kelvin_vector_size);
+            ip_data.eps.setZero(kelvin_vector_size);
+            ip_data.eps_m.setZero(kelvin_vector_size);
+            ip_data.eps_m_prev.setZero(kelvin_vector_size);
 
             ip_data.N = shape_matrices[ip].N;
             ip_data.dNdx = shape_matrices[ip].dNdx;
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
index ee734e2f5e2f0a41acb0643c0e61f95d6c6b90c1..1ec88af26ab1f69a30efdd554754847b4ea41616 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -69,8 +69,8 @@ void TwoPhaseFlowWithPPLocalAssembler<
     auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
         cap_pressure_matrix_index, cap_pressure_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
index 7211da2309c26958eb80beae0690542040541750..14d0447ae3d31441891c4179af1af04afedaad44 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
@@ -51,8 +51,8 @@ void TwoPhaseFlowWithPrhoLocalAssembler<
     auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
         cap_pressure_matrix_index, cap_pressure_matrix_index);
 
-    NodalMatrixType laplace_operator;
-    laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+    NodalMatrixType laplace_operator =
+        NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
 
     auto Kgp =
         local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
diff --git a/Tests/MaterialLib/KelvinVector.cpp b/Tests/MaterialLib/KelvinVector.cpp
index 2adb278ef470f2788dc41d0e2c77eccc491521ad..6ff6c5c609639190f396b89e1dd3b5459b049d02 100644
--- a/Tests/MaterialLib/KelvinVector.cpp
+++ b/Tests/MaterialLib/KelvinVector.cpp
@@ -154,7 +154,7 @@ TEST_F(MaterialLibSolidsKelvinVector4, Inverse)
             (inverse(v) - tensorToKelvin<4>(kelvinToTensor(v).inverse()))
                 .norm();
         // The error is only weekly depending on the input vector norm.
-        return error < 1e-6 && 1e-8 * std::pow(v.norm(), 1.4);
+        return error < 1e-6 && error < 1e-8 * std::pow(v.norm(), 1.4);
     };
 
     ac::check<KelvinVector<4>>(
@@ -174,7 +174,7 @@ TEST_F(MaterialLibSolidsKelvinVector6, Inverse)
             (inverse(v) - tensorToKelvin<6>(kelvinToTensor(v).inverse()))
                 .norm();
         // The error is only weekly depending on the input vector norm.
-        return error < 1e-6 && 1e-8 * std::pow(v.norm(), 1.4);
+        return error < 1e-6 && error < 1e-8 * std::pow(v.norm(), 1.4);
     };
 
     ac::check<KelvinVector<6>>(
diff --git a/Tests/MaterialLib/TestFractureModels.cpp b/Tests/MaterialLib/TestFractureModels.cpp
index f7496e31c14ba3a88079f5a10c560b687c3ac4a8..f21b9bd8e6e6f88dd3f87b14a2d93f823149f308 100644
--- a/Tests/MaterialLib/TestFractureModels.cpp
+++ b/Tests/MaterialLib/TestFractureModels.cpp
@@ -35,11 +35,15 @@ TEST(MaterialLib_Fracture, LinearElasticIsotropic)
     std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
         fractureModel.createMaterialStateVariables());
 
-    Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+    Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
+    Eigen::Vector2d const sigma_prev = Eigen::Vector2d::Zero();
+    Eigen::Vector2d w(-1e-5, -1e-5);
+
+    // Result vectors, not initialized.
+    Eigen::Vector2d sigma;
     Eigen::Matrix2d C;
 
     ProcessLib::SpatialPosition x;
-    w << -1e-5, -1e-5;
     fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C, *state);
 
     ASSERT_NEAR(-1e4, sigma[0], eps_sigma);
@@ -74,12 +78,15 @@ TEST(MaterialLib_Fracture, MohrCoulomb)
     std::unique_ptr<FractureModelBase<2>::MaterialStateVariables> state(
         fractureModel.createMaterialStateVariables());
 
-    Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+    Eigen::Vector2d const w_prev = Eigen::Vector2d::Zero();
+    Eigen::Vector2d const sigma_prev(-3.46e6, -2e6);
+    Eigen::Vector2d const w(-1.08e-5, -0.25e-5);
+
+    // Result vectors, not initialized.
+    Eigen::Vector2d sigma;
     Eigen::Matrix2d C;
 
     ProcessLib::SpatialPosition x;
-    sigma_prev << -3.46e6, -2e6;
-    w << -1.08e-5, -0.25e-5;
     fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C, *state);
 
     ASSERT_NEAR(-3.50360e6, sigma[0], eps_sigma);
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 51582e59f1ad9f20857c41479910b892d751978f..17c699a838f467b118790fc76884c2d5d4937f7a 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -36,6 +36,7 @@ template <class T_VECTOR>
 void checkGlobalVectorInterface()
 {
     T_VECTOR x(10);
+    x.setZero();
 
     ASSERT_EQ(10u, x.size());
     ASSERT_EQ(0u, x.getRangeBegin());
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 387dcbbef857ae2979ea274e999a05f4ade18756..23f69595d56c0da5f0440b41c43aec9943e687a0 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -117,7 +117,9 @@ void checkLinearSolverInterface(T_MATRIX &A, BaseLib::ConfigTree const& ls_optio
 
     // set RHS and solution vectors
     T_VECTOR rhs(ex1.dim_eqs);
+    rhs.setZero();
     T_VECTOR x(ex1.dim_eqs);
+    x.setZero();
 
     // apply BC
     MathLib::applyKnownSolution(A, rhs, x, ex1.vec_dirichlet_bc_id, ex1.vec_dirichlet_bc_value);
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 13a71fefa27771f9c5b45044b2afdf73a8d0def2..e2b8cb390a38234bf32285da57d0ff1b2ec9eb47 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -191,8 +191,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
     FeType fe(*this->mesh_element);
 
     // evaluate a mass matrix M = int{ N^T D N }dA_e
-    NodalMatrix M(this->e_nnodes, this->e_nnodes);
-    M.setZero();
+    NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();
@@ -216,8 +215,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
     FeType fe(*this->mesh_element);
 
     // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
-    NodalMatrix K(this->e_nnodes, this->e_nnodes);
-    K.setZero();
+    NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();
@@ -241,10 +239,8 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
     FeType fe(*this->mesh_element);
 
     // evaluate both mass and laplace matrices at once
-    NodalMatrix M(this->e_nnodes, this->e_nnodes);
-    M.setZero();
-    NodalMatrix K(this->e_nnodes, this->e_nnodes);
-    K.setZero();
+    NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
+    NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
         shape.setZero();
diff --git a/Tests/NumLib/TestODEInt.cpp b/Tests/NumLib/TestODEInt.cpp
index 63245eb7d9b546de00a616aeaa3b65821b6af864..35d420776d4e0a769ceca16fd0ad6619ebcf8802 100644
--- a/Tests/NumLib/TestODEInt.cpp
+++ b/Tests/NumLib/TestODEInt.cpp
@@ -102,7 +102,9 @@ public:
         };
 
         if (num_timesteps > 0)
+        {
             EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb));
+        }
 
         for (auto& x :  sol.solutions)
             MathLib::LinAlg::setLocalAccessibleVector(x);
diff --git a/Tests/NumLib/TestSerialLinearSolver.cpp b/Tests/NumLib/TestSerialLinearSolver.cpp
index d929a3f0f8348eb04235736669beb7fcd878d46a..cb97f08ba8bc24e7abc6c5fd98a40d3109dd6852 100644
--- a/Tests/NumLib/TestSerialLinearSolver.cpp
+++ b/Tests/NumLib/TestSerialLinearSolver.cpp
@@ -69,8 +69,9 @@ TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
     auto A = MathLib::MatrixVectorTraits<GlobalMatrix>::newInstance(ms);
     A->setZero();
     auto rhs = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(ms);
+    rhs->setZero();
     auto x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(ms);
-    // TODO no setZero() for rhs, x?
+    x->setZero();
 
     using LocalAssembler = Example::LocalAssemblerData;
     // Initializer of the local assembler data.
diff --git a/scripts/cmake/CompilerSetup.cmake b/scripts/cmake/CompilerSetup.cmake
index 5ede98b33b3b9c6918383cd87b92a66fc3062935..e6ce8afa3a0f46507b97c676cdae0d34af2ff7b2 100644
--- a/scripts/cmake/CompilerSetup.cmake
+++ b/scripts/cmake/CompilerSetup.cmake
@@ -38,9 +38,6 @@ if(BUILD_SHARED_LIBS)
     set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS TRUE)
 endif()
 
-# Set additional user-given compiler flags
-set(CMAKE_CXX_FLAGS ${OGS_CXX_FLAGS})
-
 if(OGS_ENABLE_AVX2)
     set(CPU_FLAGS "-mavx2 -march=core-avx2")
 elseif(OGS_CPU_ARCHITECTURE STREQUAL "generic")
@@ -146,3 +143,8 @@ if(MSVC)
 else()
     add_definitions(-DOPENMP_LOOP_TYPE=unsigned)
 endif()
+
+# Set additional user-given compiler flags. The given flags must follow the
+# preceding cxx flags definition in order to override earlier flags, e.g. for
+# optimization.
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OGS_CXX_FLAGS}")