diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index c0387c5a90904d22d6943b44e42af3498581cc6a..3d9797f48008ce508b2b86a231ac812c9ee94f7a 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -215,6 +215,8 @@ public:
         _A.reset(_global_setup.createMatrix(_local_to_global_index_map->dofSize()));
         _x.reset(_global_setup.createVector(_local_to_global_index_map->dofSize()));
         _rhs.reset(_global_setup.createVector(_local_to_global_index_map->dofSize()));
+        _linearSolver.reset(new typename GlobalSetup::LinearSolver(*_A));
+
 
         setInitialConditions(*_hydraulic_head);
 
@@ -261,8 +263,7 @@ public:
         // Apply known values from the Dirichlet boundary conditions.
         MathLib::applyKnownSolution(*_A, *_rhs, _dirichlet_bc.global_ids, _dirichlet_bc.values);
 
-        typename GlobalSetup::LinearSolver linearSolver(*_A);
-        linearSolver.solve(*_rhs, *_x);
+        _linearSolver->solve(*_rhs, *_x);
 
         return true;
     }
@@ -325,6 +326,7 @@ private:
     std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
 
     GlobalSetup _global_setup;
+    std::unique_ptr<typename GlobalSetup::LinearSolver> _linearSolver;
     std::unique_ptr<typename GlobalSetup::MatrixType> _A;
     std::unique_ptr<typename GlobalSetup::VectorType> _rhs;
     std::unique_ptr<typename GlobalSetup::VectorType> _x;