diff --git a/DevicePersistent/HostMatrixAdapter/HostMatrixAdapter.H b/DevicePersistent/HostMatrixAdapter/HostMatrixAdapter.H
index a7c799a04..af2c44165 100644
--- a/DevicePersistent/HostMatrixAdapter/HostMatrixAdapter.H
+++ b/DevicePersistent/HostMatrixAdapter/HostMatrixAdapter.H
@@ -83,7 +83,8 @@ struct lduMatrixInitFunctor {
             interfaceBouCoeffs_, interfaceIntCoeffs_, interfaces_,
             solver_controls_, field_name_, verbose_);
 
-        persistent_device_matrix.swap(ret);
+        persistent_device_matrix->update_from(ret);
+
     }
 
     std::shared_ptr<HostMatrixWrapper> init() const
diff --git a/MatrixWrapper/Distributed/Distributed.H b/MatrixWrapper/Distributed/Distributed.H
index 5a164432a..d33c41d67 100644
--- a/MatrixWrapper/Distributed/Distributed.H
+++ b/MatrixWrapper/Distributed/Distributed.H
@@ -311,14 +311,10 @@ public:
         label nrows = host_A->get_local_nrows();
         label local_matrix_nnz = host_A->get_local_matrix_nnz();
         // size + padding has to be local_matrix_nnz
-        std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank
-                  << "update_impl  0\n";
         auto diag_comm_pattern = compute_send_recv_counts(
             exec_handler, ranks_per_gpu, nrows, local_matrix_nnz,
             local_matrix_nnz - nrows, 0);
 
-        std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank
-                  << "update_impl  1\n";
 
         label upper_nnz = host_A->get_upper_nnz();
         auto upper_comm_pattern = compute_send_recv_counts(
@@ -328,9 +324,6 @@ public:
             compute_send_recv_counts(exec_handler, ranks_per_gpu, upper_nnz,
                                      local_matrix_nnz, upper_nnz, nrows);
 
-        std::cout << __FILE__ << ":" << __LINE__ << " rank " << rank
-                  << "update_impl  2\n";
-
         scalar *local_ptr;
 
         // update main values
diff --git a/MatrixWrapper/LDUMatrix/HostMatrix.H b/MatrixWrapper/LDUMatrix/HostMatrix.H
index 96cf37ac5..667691c93 100644
--- a/MatrixWrapper/LDUMatrix/HostMatrix.H
+++ b/MatrixWrapper/LDUMatrix/HostMatrix.H
@@ -240,11 +240,11 @@ private:
 
     const lduAddressing &addr_;
 
-    const scalar *diag_;
+    mutable const scalar *diag_;
 
-    const scalar *upper_;
+    mutable const scalar *upper_;
 
-    const scalar *lower_;
+    mutable const scalar *lower_;
 
     // multiply the complete system by this factor, ie sAx=sb
     // NOTE this needed to avoid negative diagonal matrix entries, but
@@ -364,6 +364,23 @@ public:
             fieldName, verbose);
     };
 
+    void update_from(std::shared_ptr<HostMatrixWrapper> in) {
+        diag_ = in->get_diag();
+        upper_ = in->get_upper();
+        lower_ = in->get_lower();
+
+        interface_ptr_.clear();
+        auto interfaces = in->get_interfaces();
+        for (label i = 0; i < interfaces.size(); i++) {
+            if (interface_getter(interfaces, i) == nullptr) {
+                continue;
+            }
+            const auto iface{interface_getter(interfaces, i)};
+            interface_ptr_.push_back(in->get_interfaceBouCoeffs()[i].begin());
+        }
+
+    }
+
 
     /** Based on OpenFOAMs ldu matrix format this function computes two
      ** consecutive index arrays (local_sparsisty_.row_idxs and col_idxs) in row
@@ -401,6 +418,13 @@ public:
         return interface_ptr_[i];
     };
 
+    const FieldField<Field, scalar> & get_interfaceBouCoeffs()
+        const {return interfaceBouCoeffs_;};
+
+    const lduInterfaceFieldPtrsList &get_interfaces() const {
+        return interfaces_;
+    };
+
     bool get_verbose() const { return verbose_; }
 
     label get_local_nrows() const { return nrows_; }
diff --git a/lduLduBase/lduLduBase.H b/lduLduBase/lduLduBase.H
index 89837245c..bf7f204dd 100644
--- a/lduLduBase/lduLduBase.H
+++ b/lduLduBase/lduLduBase.H
@@ -163,8 +163,6 @@ public:
           solver_controls_(solverControls),
           field_name_(fieldName),
           verbose_(solverControls.lookupOrDefault<label>("verbose", 0)),
-          // FIXME NOTE the underlying exec_handler is invalid in the next call but the
-          // exec isnt
           exec_handler_{db_, solver_controls_, fieldName},
           host_A_{
               db_,       exec_handler_, solver_controls_,   fieldName,