From 43c0c3989bec91d1c185d64802c0d63a4f0b66f5 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 16 Feb 2024 17:35:37 +0100
Subject: [PATCH] ENH: handle matrix demand-driven internals with unique_ptr

- align member order between template and non-template versions
- move construct, move assignment
---
 .../matrices/LduMatrix/LduMatrix/LduMatrix.C  | 298 +++++++---------
 .../matrices/LduMatrix/LduMatrix/LduMatrix.H  | 192 +++++-----
 .../LduMatrix/LduMatrix/LduMatrixOperations.C |  60 ++--
 .../lduMatrix/lduAddressing/lduAddressing.C   |  38 +-
 .../lduMatrix/lduAddressing/lduAddressing.H   |  42 +--
 .../matrices/lduMatrix/lduMatrix/lduMatrix.C  | 335 +++++++++---------
 .../matrices/lduMatrix/lduMatrix/lduMatrix.H  | 150 ++++----
 .../lduMatrix/lduMatrix/lduMatrixOperations.C |  77 ++--
 8 files changed, 600 insertions(+), 592 deletions(-)

diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
index dbdaa0fbc3f..c76026bd535 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,109 +34,86 @@ License
 template<class Type, class DType, class LUType>
 Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const lduMesh& mesh)
 :
-    lduMesh_(mesh),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr),
-    lowerPtr_(nullptr),
-    sourcePtr_(nullptr),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
+    lduMesh_(mesh)
 {}
 
 
 template<class Type, class DType, class LUType>
 Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const LduMatrix& A)
 :
-    lduMesh_(A.lduMesh_),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr),
-    lowerPtr_(nullptr),
-    sourcePtr_(nullptr),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
+    lduMesh_(A.lduMesh_)
 {
     if (A.diagPtr_)
     {
-        diagPtr_ = new Field<DType>(*(A.diagPtr_));
+        diagPtr_ = std::make_unique<Field<DType>>(*(A.diagPtr_));
     }
 
     if (A.upperPtr_)
     {
-        upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+        upperPtr_ = std::make_unique<Field<LUType>>(*(A.upperPtr_));
     }
 
     if (A.lowerPtr_)
     {
-        lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+        lowerPtr_ = std::make_unique<Field<LUType>>(*(A.lowerPtr_));
     }
 
     if (A.sourcePtr_)
     {
-        sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+        sourcePtr_ = std::make_unique<Field<Type>>(*(A.sourcePtr_));
     }
 }
 
 
 template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reuse)
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix&& A)
 :
     lduMesh_(A.lduMesh_),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr),
-    lowerPtr_(nullptr),
-    sourcePtr_(nullptr),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
+    diagPtr_(std::move(A.diagPtr_)),
+    lowerPtr_(std::move(A.lowerPtr_)),
+    upperPtr_(std::move(A.upperPtr_)),
+    sourcePtr_(std::move(A.sourcePtr_))
 {
-    if (reuse)
-    {
-        if (A.diagPtr_)
-        {
-            diagPtr_ = A.diagPtr_;
-            A.diagPtr_ = nullptr;
-        }
+    // Clear the old interfaces?
+}
 
-        if (A.upperPtr_)
-        {
-            upperPtr_ = A.upperPtr_;
-            A.upperPtr_ = nullptr;
-        }
 
-        if (A.lowerPtr_)
-        {
-            lowerPtr_ = A.lowerPtr_;
-            A.lowerPtr_ = nullptr;
-        }
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reuse)
+:
+    lduMesh_(A.lduMesh_)
+{
+    if (reuse)
+    {
+        // Move assignment
+        diagPtr_ = std::move(A.diagPtr_);
+        upperPtr_ = std::move(A.upperPtr_);
+        lowerPtr_ = std::move(A.lowerPtr_);
+        sourcePtr_ = std::move(A.sourcePtr_);
 
-        if (A.sourcePtr_)
-        {
-            sourcePtr_ = A.sourcePtr_;
-            A.sourcePtr_ = nullptr;
-        }
+        // Clear the old interfaces?
     }
     else
     {
+        // Copy assignment
         if (A.diagPtr_)
         {
-            diagPtr_ = new Field<DType>(*(A.diagPtr_));
+            diagPtr_ = std::make_unique<Field<DType>>(*(A.diagPtr_));
         }
 
         if (A.upperPtr_)
         {
-            upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+            upperPtr_ = std::make_unique<Field<LUType>>(*(A.upperPtr_));
         }
 
         if (A.lowerPtr_)
         {
-            lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+            lowerPtr_ = std::make_unique<Field<LUType>>(*(A.lowerPtr_));
         }
 
         if (A.sourcePtr_)
         {
-            sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+            sourcePtr_ = std::make_unique<Field<Type>>(*(A.sourcePtr_));
         }
     }
 }
@@ -152,48 +130,51 @@ Foam::LduMatrix<Type, DType, LUType>::LduMatrix
     diagPtr_(new Field<DType>(is)),
     upperPtr_(new Field<LUType>(is)),
     lowerPtr_(new Field<LUType>(is)),
-    sourcePtr_(new Field<Type>(is)),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
+    sourcePtr_(new Field<Type>(is))
 {}
 
 
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
+Foam::word Foam::LduMatrix<Type, DType, LUType>::matrixTypeName() const
 {
     if (diagPtr_)
     {
-        delete diagPtr_;
+        return
+        (
+            (!upperPtr_)
+          ? (!lowerPtr_ ? "diagonal" : "diagonal-lower")
+          : (!lowerPtr_ ? "symmetric" : "asymmetric")
+        );
     }
 
-    if (upperPtr_)
-    {
-        delete upperPtr_;
-    }
+    // is empty (or just wrong)
+    return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined");
+}
 
-    if (lowerPtr_)
-    {
-        delete lowerPtr_;
-    }
 
-    if (sourcePtr_)
+template<class Type, class DType, class LUType>
+const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
+{
+    if (!diagPtr_)
     {
-        delete sourcePtr_;
+        FatalErrorInFunction
+            << "diagPtr_ unallocated"
+            << abort(FatalError);
     }
-}
 
+    return *diagPtr_;
+}
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
 Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
 {
     if (!diagPtr_)
     {
-        diagPtr_ = new Field<DType>(lduAddr().size(), Zero);
+        diagPtr_ =
+            std::make_unique<Field<DType>>(lduAddr().size(), Foam::zero{});
     }
 
     return *diagPtr_;
@@ -201,127 +182,115 @@ Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
 
 
 template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
 {
-    if (!upperPtr_)
+    if (upperPtr_)
     {
-        if (lowerPtr_)
-        {
-            upperPtr_ = new Field<LUType>(*lowerPtr_);
-        }
-        else
+        return *upperPtr_;
+    }
+    else
+    {
+        if (!lowerPtr_)
         {
-            upperPtr_ = new Field<LUType>
-            (
-                lduAddr().lowerAddr().size(),
-                Zero
-            );
+            FatalErrorInFunction
+                << "lowerPtr_ and upperPtr_ unallocated"
+                << abort(FatalError);
         }
-    }
 
-    return *upperPtr_;
+        return *lowerPtr_;
+    }
 }
 
 
 template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
 {
-    if (!lowerPtr_)
+    if (!upperPtr_)
     {
-        if (upperPtr_)
+        if (lowerPtr_)
         {
-            lowerPtr_ = new Field<LUType>(*upperPtr_);
+            upperPtr_ = std::make_unique<Field<LUType>>(*lowerPtr_);
         }
         else
         {
-            lowerPtr_ = new Field<LUType>
-            (
-                lduAddr().lowerAddr().size(),
-                Zero
-            );
+            upperPtr_ =
+                std::make_unique<Field<LUType>>
+                (
+                    lduAddr().lowerAddr().size(),
+                    Foam::zero{}
+                );
         }
     }
 
-    return *lowerPtr_;
+    return *upperPtr_;
 }
 
 
 template<class Type, class DType, class LUType>
-Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
 {
-    if (!sourcePtr_)
+    if (lowerPtr_)
     {
-        sourcePtr_ = new Field<Type>(lduAddr().size(), Zero);
+        return *lowerPtr_;
     }
-
-    return *sourcePtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
-{
-    if (!diagPtr_)
+    else
     {
-        FatalErrorInFunction
-            << "diagPtr_ unallocated"
-            << abort(FatalError);
-    }
+        if (!upperPtr_)
+        {
+            FatalErrorInFunction
+                << "lowerPtr_ and upperPtr_ unallocated"
+                << abort(FatalError);
+        }
 
-    return *diagPtr_;
+        return *upperPtr_;
+    }
 }
 
 
 template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
 {
-    if (!lowerPtr_ && !upperPtr_)
+    if (!lowerPtr_)
     {
-        FatalErrorInFunction
-            << "lowerPtr_ or upperPtr_ unallocated"
-            << abort(FatalError);
+        if (upperPtr_)
+        {
+            lowerPtr_ = std::make_unique<Field<LUType>>(*upperPtr_);
+        }
+        else
+        {
+            lowerPtr_ = std::make_unique<Field<LUType>>
+            (
+                lduAddr().lowerAddr().size(),
+                Foam::zero{}
+            );
+        }
     }
 
-    if (upperPtr_)
-    {
-        return *upperPtr_;
-    }
-    else
-    {
-        return *lowerPtr_;
-    }
+    return *lowerPtr_;
 }
 
 
 template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
+const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
 {
-    if (!lowerPtr_ && !upperPtr_)
+    if (!sourcePtr_)
     {
         FatalErrorInFunction
-            << "lowerPtr_ or upperPtr_ unallocated"
+            << "sourcePtr_ unallocated"
             << abort(FatalError);
     }
 
-    if (lowerPtr_)
-    {
-        return *lowerPtr_;
-    }
-    else
-    {
-        return *upperPtr_;
-    }
+    return *sourcePtr_;
 }
 
 
 template<class Type, class DType, class LUType>
-const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
+Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
 {
     if (!sourcePtr_)
     {
-        FatalErrorInFunction
-            << "sourcePtr_ unallocated"
-            << abort(FatalError);
+        sourcePtr_ =
+            std::make_unique<Field<Type>>(lduAddr().size(), Foam::zero{});
     }
 
     return *sourcePtr_;
@@ -330,43 +299,50 @@ const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
+// template<class Type, class DType, class LUType>
+// Foam::Ostream& Foam::operator<<
+// (
+//     Ostream& os,
+//     const InfoProxy<Type, DType, LUType>& iproxy
+// )
+// {
+//     const auto& mat = *iproxy;
+//
+//     ...
+//
+//     os.check(FUNCTION_NAME);
+//     return os;
+// }
+
+
 template<class Type, class DType, class LUType>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const LduMatrix<Type, DType, LUType>& ldum
+    const LduMatrix<Type, DType, LUType>& mat
 )
 {
-    if (ldum.diagPtr_)
+    if (mat.hasDiag())
     {
-        os  << "Diagonal = "
-            << *ldum.diagPtr_
-            << endl << endl;
+        os  << "Diagonal = " << mat.diag() << nl << nl;
     }
 
-    if (ldum.upperPtr_)
+    if (mat.hasUpper())
     {
-        os  << "Upper triangle = "
-            << *ldum.upperPtr_
-            << endl << endl;
+        os  << "Upper triangle = " << mat.upper() << nl << nl;
     }
 
-    if (ldum.lowerPtr_)
+    if (mat.hasLower())
     {
-        os  << "Lower triangle = "
-            << *ldum.lowerPtr_
-            << endl << endl;
+        os  << "Lower triangle = " << mat.lower() << nl << nl;
     }
 
-    if (ldum.sourcePtr_)
+    if (mat.hasSource())
     {
-        os  << "Source = "
-            << *ldum.sourcePtr_
-            << endl << endl;
+        os  << "Source = " << mat.source() << nl << nl;
     }
 
     os.check(FUNCTION_NAME);
-
     return os;
 }
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
index 96e9fd14f4e..407f7ee3226 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -93,25 +93,32 @@ class LduMatrix
         const lduMesh& lduMesh_;
 
         //- Diagonal coefficients
-        Field<DType> *diagPtr_;
+        std::unique_ptr<Field<DType>> diagPtr_;
 
         //- Off-diagonal coefficients
-        Field<LUType> *upperPtr_, *lowerPtr_;
+        std::unique_ptr<Field<LUType>> upperPtr_;
+
+        //- Off-diagonal coefficients
+        std::unique_ptr<Field<LUType>> lowerPtr_;
 
         //- Source
-        Field<Type> *sourcePtr_;
+        std::unique_ptr<Field<Type>> sourcePtr_;
 
         //- Field interfaces (processor patches etc.)
         LduInterfaceFieldPtrsList<Type> interfaces_;
 
         //- Off-diagonal coefficients for interfaces
-        FieldField<Field, LUType> interfacesUpper_, interfacesLower_;
+        FieldField<Field, LUType> interfacesUpper_;
+
+        //- Off-diagonal coefficients for interfaces
+        FieldField<Field, LUType> interfacesLower_;
 
 
 public:
 
     friend class SolverPerformance<Type>;
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for LduMatrix solvers
     class solver
     {
@@ -274,6 +281,7 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for LduMatrix smoothers
     class smoother
     {
@@ -371,6 +379,7 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for LduMatrix preconditioners
     class preconditioner
     {
@@ -466,6 +475,8 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
+
     // Static Data
 
         // Declare name of the class and its debug switch
@@ -476,126 +487,125 @@ public:
 
         //- Construct given an LDU addressed mesh.
         //  The coefficients are initially empty for subsequent setting.
-        LduMatrix(const lduMesh&);
+        //  Not yet 'explicit' (legacy code may rely on implicit construct)
+        LduMatrix(const lduMesh& mesh);
 
-        //- Construct as copy
+        //- Copy construct
         LduMatrix(const LduMatrix<Type, DType, LUType>&);
 
+        //- Move construct
+        LduMatrix(LduMatrix<Type, DType, LUType>&&);
+
         //- Construct as copy or re-use as specified.
         LduMatrix(LduMatrix<Type, DType, LUType>&, bool reuse);
 
         //- Construct given an LDU addressed mesh and an Istream
-        //  from which the coefficients are read
-        LduMatrix(const lduMesh&, Istream&);
+        //- from which the coefficients are read
+        LduMatrix(const lduMesh& mesh, Istream& is);
 
 
     //- Destructor
-    ~LduMatrix();
+    ~LduMatrix() = default;
 
 
     // Member Functions
 
-        // Access to addressing
+    // Addressing
 
-            //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const noexcept
-            {
-                return lduMesh_;
-            }
+        //- Return the LDU mesh from which the addressing is obtained
+        const lduMesh& mesh() const noexcept
+        {
+            return lduMesh_;
+        }
 
-            //- Return the LDU addressing
-            const lduAddressing& lduAddr() const
-            {
-                return lduMesh_.lduAddr();
-            }
+        //- Return the LDU addressing
+        const lduAddressing& lduAddr() const
+        {
+            return lduMesh_.lduAddr();
+        }
 
-            //- Return the patch evaluation schedule
-            const lduSchedule& patchSchedule() const
-            {
-                return lduAddr().patchSchedule();
-            }
+        //- Return the patch evaluation schedule
+        const lduSchedule& patchSchedule() const
+        {
+            return lduMesh_.lduAddr().patchSchedule();
+        }
 
-            //- Return interfaces
-            const LduInterfaceFieldPtrsList<Type>& interfaces() const
-            {
-                return interfaces_;
-            }
+        //- Const access to the interfaces
+        const LduInterfaceFieldPtrsList<Type>& interfaces() const noexcept
+        {
+            return interfaces_;
+        }
 
-            //- Return interfaces
-            LduInterfaceFieldPtrsList<Type>& interfaces()
-            {
-                return interfaces_;
-            }
+        //- Non-const access to the interfaces
+        LduInterfaceFieldPtrsList<Type>& interfaces() noexcept
+        {
+            return interfaces_;
+        }
 
 
-        // Access to coefficients
+    // Coefficients
 
-            Field<DType>& diag();
-            Field<LUType>& upper();
-            Field<LUType>& lower();
-            Field<Type>& source();
+        const Field<DType>& diag() const;
+        const Field<LUType>& upper() const;
+        const Field<LUType>& lower() const;
+        const Field<Type>& source() const;
 
-            FieldField<Field, LUType>& interfacesUpper()
-            {
-                return interfacesUpper_;
-            }
+        Field<DType>& diag();
+        Field<LUType>& upper();
+        Field<LUType>& lower();
+        Field<Type>& source();
 
-            FieldField<Field, LUType>& interfacesLower()
-            {
-                return interfacesLower_;
-            }
 
+    // Interfaces
 
-            const Field<DType>& diag() const;
-            const Field<LUType>& upper() const;
-            const Field<LUType>& lower() const;
-            const Field<Type>& source() const;
+        const FieldField<Field, LUType>& interfacesUpper() const noexcept
+        {
+            return interfacesUpper_;
+        }
 
-            const FieldField<Field, LUType>& interfacesUpper() const
-            {
-                return interfacesUpper_;
-            }
+        const FieldField<Field, LUType>& interfacesLower() const noexcept
+        {
+            return interfacesLower_;
+        }
 
-            const FieldField<Field, LUType>& interfacesLower() const
-            {
-                return interfacesLower_;
-            }
+        FieldField<Field, LUType>& interfacesUpper() noexcept
+        {
+            return interfacesUpper_;
+        }
 
+        FieldField<Field, LUType>& interfacesLower() noexcept
+        {
+            return interfacesLower_;
+        }
 
-            bool hasDiag() const noexcept
-            {
-                return (diagPtr_);
-            }
 
-            bool hasUpper() const noexcept
-            {
-                return (upperPtr_);
-            }
+    // Characteristics
 
-            bool hasLower() const noexcept
-            {
-                return (lowerPtr_);
-            }
+        //- The matrix type (empty, diagonal, symmetric, ...)
+        word matrixTypeName() const;
 
-            bool hasSource() const noexcept
-            {
-                return (sourcePtr_);
-            }
+        bool hasDiag() const noexcept { return bool(diagPtr_); }
+        bool hasUpper() const noexcept { return bool(upperPtr_); }
+        bool hasLower() const noexcept { return bool(lowerPtr_); }
+        bool hasSource() const noexcept { return bool(sourcePtr_); }
 
-            bool diagonal() const noexcept
-            {
-                return (diagPtr_ && !lowerPtr_ && !upperPtr_);
-            }
+        //- Matrix has diagonal only
+        bool diagonal() const noexcept
+        {
+            return (diagPtr_ && !lowerPtr_ && !upperPtr_);
+        }
 
-            bool symmetric() const noexcept
-            {
-                return (diagPtr_ && (!lowerPtr_ && upperPtr_));
-            }
+        //- Matrix is symmetric
+        bool symmetric() const noexcept
+        {
+            return (diagPtr_ && !lowerPtr_ && upperPtr_);
+        }
 
-            bool asymmetric() const noexcept
-            {
-                return (diagPtr_ && lowerPtr_ && upperPtr_);
-            }
+        //- Matrix is asymmetric (ie, full)
+        bool asymmetric() const noexcept
+        {
+            return (diagPtr_ && lowerPtr_ && upperPtr_);
+        }
 
 
         // Operations
@@ -651,8 +661,12 @@ public:
 
     // Member Operators
 
+        //- Copy assignment
         void operator=(const LduMatrix<Type, DType, LUType>&);
 
+        //- Move assignment
+        void operator=(LduMatrix<Type, DType, LUType>&&);
+
         void negate();
 
         void operator+=(const LduMatrix<Type, DType, LUType>&);
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
index 7945a8cd2da..46e7ad7b605 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,7 +92,7 @@ Foam::LduMatrix<Type, DType, LUType>::H(const Field<Type>& psi) const
 {
     auto tHpsi = tmp<Field<Type>>::New(lduAddr().size(), Foam::zero{});
 
-    if (lowerPtr_ || upperPtr_)
+    if (hasLower() || hasUpper())
     {
         Type* __restrict__ HpsiPtr = tHpsi.ref().begin();
 
@@ -169,32 +169,30 @@ void Foam::LduMatrix<Type, DType, LUType>::operator=(const LduMatrix& A)
         return;  // Self-assignment is a no-op
     }
 
-    if (A.diagPtr_)
+    if (A.hasDiag())
     {
         diag() = A.diag();
     }
 
-    if (A.upperPtr_)
+    if (A.hasUpper())
     {
         upper() = A.upper();
     }
-    else if (upperPtr_)
+    else
     {
-        delete upperPtr_;
-        upperPtr_ = nullptr;
+        upperPtr_.reset(nullptr);
     }
 
-    if (A.lowerPtr_)
+    if (A.hasLower())
     {
         lower() = A.lower();
     }
-    else if (lowerPtr_)
+    else
     {
-        delete lowerPtr_;
-        lowerPtr_ = nullptr;
+        lowerPtr_.reset(nullptr);
     }
 
-    if (A.sourcePtr_)
+    if (A.hasSource())
     {
         source() = A.source();
     }
@@ -204,6 +202,24 @@ void Foam::LduMatrix<Type, DType, LUType>::operator=(const LduMatrix& A)
 }
 
 
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator=(LduMatrix&& A)
+{
+    if (this == &A)
+    {
+        return;  // Self-assignment is a no-op
+    }
+
+    diagPtr_ = std::move(A.diagPtr_);
+    upperPtr_ = std::move(A.upperPtr_);
+    lowerPtr_ = std::move(A.lowerPtr_);
+    sourcePtr_ = std::move(A.sourcePtr_);
+
+    interfacesUpper_ = std::move(A.interfacesUpper_);
+    interfacesLower_ = std::move(A.interfacesLower_);
+}
+
+
 template<class Type, class DType, class LUType>
 void Foam::LduMatrix<Type, DType, LUType>::negate()
 {
@@ -235,12 +251,12 @@ void Foam::LduMatrix<Type, DType, LUType>::negate()
 template<class Type, class DType, class LUType>
 void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
 {
-    if (A.diagPtr_)
+    if (A.hasDiag())
     {
         diag() += A.diag();
     }
 
-    if (A.sourcePtr_)
+    if (A.hasSource())
     {
         source() += A.source();
     }
@@ -265,7 +281,7 @@ void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
     }
     else if (asymmetric() && A.symmetric())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             lower() += A.upper();
             upper() += A.upper();
@@ -284,12 +300,12 @@ void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
     }
     else if (diagonal())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             upper() = A.upper();
         }
 
-        if (A.lowerPtr_)
+        if (A.hasLower())
         {
             lower() = A.lower();
         }
@@ -312,12 +328,12 @@ void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
 template<class Type, class DType, class LUType>
 void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
 {
-    if (A.diagPtr_)
+    if (A.hasDiag())
     {
         diag() -= A.diag();
     }
 
-    if (A.sourcePtr_)
+    if (A.hasSource())
     {
         source() -= A.source();
     }
@@ -342,7 +358,7 @@ void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
     }
     else if (asymmetric() && A.symmetric())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             lower() -= A.upper();
             upper() -= A.upper();
@@ -361,12 +377,12 @@ void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
     }
     else if (diagonal())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             upper() = -A.upper();
         }
 
-        if (A.lowerPtr_)
+        if (A.hasLower())
         {
             lower() = -A.lower();
         }
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C
index 61006d077c3..91180115b7f 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "lduAddressing.H"
-#include "demandDrivenData.H"
 #include "scalarField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -44,7 +43,7 @@ void Foam::lduAddressing::calcLosort() const
     // Scan the neighbour list to find out how many times the cell
     // appears as a neighbour of the face. Done this way to avoid guessing
     // and resizing list
-    labelList nNbrOfFace(size(), Zero);
+    labelList nNbrOfFace(size(), Foam::zero{});
 
     const labelUList& nbr = upperAddr();
 
@@ -73,9 +72,8 @@ void Foam::lduAddressing::calcLosort() const
     }
 
     // Gather the neighbours into the losort array
-    losortPtr_ = new labelList(nbr.size(), -1);
-
-    labelList& lst = *losortPtr_;
+    losortPtr_ = std::make_unique<labelList>(nbr.size(), -1);
+    auto& lst = *losortPtr_;
 
     // Set counter for losort
     label lstI = 0;
@@ -104,9 +102,8 @@ void Foam::lduAddressing::calcOwnerStart() const
 
     const labelList& own = lowerAddr();
 
-    ownerStartPtr_ = new labelList(size() + 1, own.size());
-
-    labelList& ownStart = *ownerStartPtr_;
+    ownerStartPtr_ = std::make_unique<labelList>(size() + 1, own.size());
+    auto& ownStart = *ownerStartPtr_;
 
     // Set up first lookup by hand
     ownStart[0] = 0;
@@ -139,9 +136,8 @@ void Foam::lduAddressing::calcLosortStart() const
             << abort(FatalError);
     }
 
-    losortStartPtr_ = new labelList(size() + 1, Zero);
-
-    labelList& lsrtStart = *losortStartPtr_;
+    losortStartPtr_ = std::make_unique<labelList>(size() + 1, Foam::zero{});
+    auto& lsrtStart = *losortStartPtr_;
 
     const labelList& nbr = upperAddr();
 
@@ -173,16 +169,6 @@ void Foam::lduAddressing::calcLosortStart() const
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::lduAddressing::~lduAddressing()
-{
-    deleteDemandDrivenData(losortPtr_);
-    deleteDemandDrivenData(ownerStartPtr_);
-    deleteDemandDrivenData(losortStartPtr_);
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::labelUList& Foam::lduAddressing::losortAddr() const
@@ -220,9 +206,9 @@ const Foam::labelUList& Foam::lduAddressing::losortStartAddr() const
 
 void Foam::lduAddressing::clearOut()
 {
-    deleteDemandDrivenData(losortPtr_);
-    deleteDemandDrivenData(ownerStartPtr_);
-    deleteDemandDrivenData(losortStartPtr_);
+    losortPtr_.reset(nullptr);
+    ownerStartPtr_.reset(nullptr);
+    losortStartPtr_.reset(nullptr);
 }
 
 
@@ -262,7 +248,7 @@ Foam::Tuple2<Foam::label, Foam::scalar> Foam::lduAddressing::band() const
     const labelUList& owner = lowerAddr();
     const labelUList& neighbour = upperAddr();
 
-    labelList cellBandwidth(size(), Zero);
+    labelList cellBandwidth(size(), Foam::zero{});
 
     forAll(neighbour, facei)
     {
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H
index c37cd61aec8..b610a6539d4 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -123,23 +123,17 @@ class lduAddressing
     //- Demand-driven data
 
         //- Losort addressing
-        mutable labelList* losortPtr_;
+        mutable std::unique_ptr<labelList> losortPtr_;
 
         //- Owner start addressing
-        mutable labelList* ownerStartPtr_;
+        mutable std::unique_ptr<labelList> ownerStartPtr_;
 
         //- Losort start addressing
-        mutable labelList* losortStartPtr_;
+        mutable std::unique_ptr<labelList> losortStartPtr_;
 
 
     // Private Member Functions
 
-        //- No copy construct
-        lduAddressing(const lduAddressing&) = delete;
-
-        //- No copy assignment
-        void operator=(const lduAddressing&) = delete;
-
         //- Calculate losort
         void calcLosort() const;
 
@@ -152,24 +146,32 @@ class lduAddressing
 
 public:
 
-    //- Construct with size (number of equations)
-    explicit lduAddressing(const label nEqns)
-    :
-        size_(nEqns),
-        losortPtr_(nullptr),
-        ownerStartPtr_(nullptr),
-        losortStartPtr_(nullptr)
-    {}
+    // Generated Methods
+
+        //- No copy construct
+        lduAddressing(const lduAddressing&) = delete;
+
+        //- No copy assignment
+        void operator=(const lduAddressing&) = delete;
+
+
+    // Constructors
+
+        //- Construct with size (number of equations)
+        explicit lduAddressing(const label nEqns) noexcept
+        :
+            size_(nEqns)
+        {}
 
 
     //- Destructor
-    virtual ~lduAddressing();
+    virtual ~lduAddressing() = default;
 
 
     // Member Functions
 
         //- Return number of equations
-        label size() const
+        label size() const noexcept
         {
             return size_;
         }
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
index dd323254199..e2767c6b526 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,79 +60,66 @@ Foam::lduMatrix::normTypesNames_
 
 Foam::lduMatrix::lduMatrix(const lduMesh& mesh)
 :
-    lduMesh_(mesh),
-    lowerPtr_(nullptr),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr)
+    lduMesh_(mesh)
 {}
 
 
 Foam::lduMatrix::lduMatrix(const lduMatrix& A)
 :
-    lduMesh_(A.lduMesh_),
-    lowerPtr_(nullptr),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr)
+    lduMesh_(A.lduMesh_)
 {
-    if (A.lowerPtr_)
+    if (A.diagPtr_)
     {
-        lowerPtr_ = new scalarField(*(A.lowerPtr_));
+        diagPtr_ = std::make_unique<scalarField>(*(A.diagPtr_));
     }
 
-    if (A.diagPtr_)
+    if (A.upperPtr_)
     {
-        diagPtr_ = new scalarField(*(A.diagPtr_));
+        upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
     }
 
-    if (A.upperPtr_)
+    if (A.lowerPtr_)
     {
-        upperPtr_ = new scalarField(*(A.upperPtr_));
+        lowerPtr_ = std::make_unique<scalarField>(*(A.lowerPtr_));
     }
 }
 
 
-Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse)
+Foam::lduMatrix::lduMatrix(lduMatrix&& A)
 :
     lduMesh_(A.lduMesh_),
-    lowerPtr_(nullptr),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr)
+    diagPtr_(std::move(A.diagPtr_)),
+    lowerPtr_(std::move(A.lowerPtr_)),
+    upperPtr_(std::move(A.upperPtr_))
+{}
+
+
+Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse)
+:
+    lduMesh_(A.lduMesh_)
 {
     if (reuse)
     {
-        if (A.lowerPtr_)
-        {
-            lowerPtr_ = A.lowerPtr_;
-            A.lowerPtr_ = nullptr;
-        }
-
-        if (A.diagPtr_)
-        {
-            diagPtr_ = A.diagPtr_;
-            A.diagPtr_ = nullptr;
-        }
-
-        if (A.upperPtr_)
-        {
-            upperPtr_ = A.upperPtr_;
-            A.upperPtr_ = nullptr;
-        }
+        diagPtr_ = std::move(A.diagPtr_);
+        upperPtr_ = std::move(A.upperPtr_);
+        lowerPtr_ = std::move(A.lowerPtr_);
     }
     else
     {
-        if (A.lowerPtr_)
+        // Copy assignment
+        if (A.diagPtr_)
         {
-            lowerPtr_ = new scalarField(*(A.lowerPtr_));
+            diagPtr_ = std::make_unique<scalarField>(*(A.diagPtr_));
         }
 
-        if (A.diagPtr_)
+        if (A.upperPtr_)
         {
-            diagPtr_ = new scalarField(*(A.diagPtr_));
+            upperPtr_ = std::make_unique<scalarField>(*(A.upperPtr_));
         }
 
-        if (A.upperPtr_)
+        if (A.lowerPtr_)
         {
-            upperPtr_ = new scalarField(*(A.upperPtr_));
+            lowerPtr_ = std::make_unique<scalarField>(*(A.lowerPtr_));
         }
     }
 }
@@ -140,64 +127,56 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reuse)
 
 Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is)
 :
-    lduMesh_(mesh),
-    lowerPtr_(nullptr),
-    diagPtr_(nullptr),
-    upperPtr_(nullptr)
+    lduMesh_(mesh)
 {
-    Switch hasLow(is);
-    Switch hasDiag(is);
-    Switch hasUp(is);
+    bool withLower, withDiag, withUpper;
+
+    is >> withLower >> withDiag >> withUpper;
 
-    if (hasLow)
+    if (withLower)
     {
-        lowerPtr_ = new scalarField(is);
+        lowerPtr_ = std::make_unique<scalarField>(is);
     }
-    if (hasDiag)
+    if (withDiag)
     {
-        diagPtr_ = new scalarField(is);
+        diagPtr_ = std::make_unique<scalarField>(is);
     }
-    if (hasUp)
+    if (withUpper)
     {
-        upperPtr_ = new scalarField(is);
+        upperPtr_ = std::make_unique<scalarField>(is);
     }
 }
 
 
-Foam::lduMatrix::~lduMatrix()
-{
-    if (lowerPtr_)
-    {
-        delete lowerPtr_;
-    }
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::word Foam::lduMatrix::matrixTypeName() const
+{
     if (diagPtr_)
     {
-        delete diagPtr_;
+        return
+        (
+            (!upperPtr_)
+          ? (!lowerPtr_ ? "diagonal" : "diagonal-lower")
+          : (!lowerPtr_ ? "symmetric" : "asymmetric")
+        );
     }
 
-    if (upperPtr_)
-    {
-        delete upperPtr_;
-    }
+    // is empty (or just wrong)
+    return (!upperPtr_ && !lowerPtr_ ? "empty" : "ill-defined");
 }
 
 
-Foam::scalarField& Foam::lduMatrix::lower()
+const Foam::scalarField& Foam::lduMatrix::diag() const
 {
-    if (!lowerPtr_)
+    if (!diagPtr_)
     {
-        if (upperPtr_)
-        {
-            lowerPtr_ = new scalarField(*upperPtr_);
-        }
-        else
-        {
-            lowerPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero);
-        }
+        FatalErrorInFunction
+            << "diagPtr_ unallocated"
+            << abort(FatalError);
     }
 
-    return *lowerPtr_;
+    return *diagPtr_;
 }
 
 
@@ -205,71 +184,87 @@ Foam::scalarField& Foam::lduMatrix::diag()
 {
     if (!diagPtr_)
     {
-        diagPtr_ = new scalarField(lduAddr().size(), Zero);
+        diagPtr_ =
+            std::make_unique<scalarField>(lduAddr().size(), Foam::zero{});
     }
 
     return *diagPtr_;
 }
 
 
-Foam::scalarField& Foam::lduMatrix::upper()
+Foam::scalarField& Foam::lduMatrix::diag(label size)
 {
-    if (!upperPtr_)
+    if (!diagPtr_)
     {
-        if (lowerPtr_)
-        {
-            upperPtr_ = new scalarField(*lowerPtr_);
-        }
-        else
-        {
-            upperPtr_ = new scalarField(lduAddr().lowerAddr().size(), Zero);
-        }
+        // if (size < 0)
+        // {
+        //     size = lduAddr().size();
+        // }
+        diagPtr_ = std::make_unique<scalarField>(size, Foam::zero{});
     }
 
-    return *upperPtr_;
+    return *diagPtr_;
 }
 
 
-Foam::scalarField& Foam::lduMatrix::lower(const label nCoeffs)
+const Foam::scalarField& Foam::lduMatrix::upper() const
 {
-    if (!lowerPtr_)
+    if (upperPtr_)
     {
-        if (upperPtr_)
-        {
-            lowerPtr_ = new scalarField(*upperPtr_);
-        }
-        else
+        return *upperPtr_;
+    }
+    else
+    {
+        if (!lowerPtr_)
         {
-            lowerPtr_ = new scalarField(nCoeffs, Zero);
+            FatalErrorInFunction
+                << "lowerPtr_ and upperPtr_ unallocated"
+                << abort(FatalError);
         }
-    }
 
-    return *lowerPtr_;
+        return *lowerPtr_;
+    }
 }
 
 
-Foam::scalarField& Foam::lduMatrix::diag(const label size)
+Foam::scalarField& Foam::lduMatrix::upper()
 {
-    if (!diagPtr_)
+    if (!upperPtr_)
     {
-        diagPtr_ = new scalarField(size, Zero);
+        if (lowerPtr_)
+        {
+            upperPtr_ = std::make_unique<scalarField>(*lowerPtr_);
+        }
+        else
+        {
+            upperPtr_ =
+                std::make_unique<scalarField>
+                (
+                    lduAddr().lowerAddr().size(),
+                    Foam::zero{}
+                );
+        }
     }
 
-    return *diagPtr_;
+    return *upperPtr_;
 }
 
 
-Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs)
+Foam::scalarField& Foam::lduMatrix::upper(label nCoeffs)
 {
     if (!upperPtr_)
     {
         if (lowerPtr_)
         {
-            upperPtr_ = new scalarField(*lowerPtr_);
+            upperPtr_ = std::make_unique<scalarField>(*lowerPtr_);
         }
         else
         {
-            upperPtr_ = new scalarField(nCoeffs, Zero);
+            // if (nCoeffs < 0)
+            // {
+            //     nCoeffs = lduAddr().lowerAddr().size();
+            // }
+            upperPtr_ = std::make_unique<scalarField>(nCoeffs, Foam::zero{});
         }
     }
 
@@ -279,54 +274,67 @@ Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs)
 
 const Foam::scalarField& Foam::lduMatrix::lower() const
 {
-    if (!lowerPtr_ && !upperPtr_)
-    {
-        FatalErrorInFunction
-            << "lowerPtr_ or upperPtr_ unallocated"
-            << abort(FatalError);
-    }
-
     if (lowerPtr_)
     {
         return *lowerPtr_;
     }
     else
     {
+        if (!upperPtr_)
+        {
+            FatalErrorInFunction
+                << "lowerPtr_ and upperPtr_ unallocated"
+                << abort(FatalError);
+        }
+
         return *upperPtr_;
     }
 }
 
 
-const Foam::scalarField& Foam::lduMatrix::diag() const
+Foam::scalarField& Foam::lduMatrix::lower()
 {
-    if (!diagPtr_)
+    if (!lowerPtr_)
     {
-        FatalErrorInFunction
-            << "diagPtr_ unallocated"
-            << abort(FatalError);
+        if (upperPtr_)
+        {
+            lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
+        }
+        else
+        {
+            lowerPtr_ =
+                std::make_unique<scalarField>
+                (
+                    lduAddr().lowerAddr().size(),
+                    Foam::zero{}
+                );
+        }
     }
 
-    return *diagPtr_;
+    return *lowerPtr_;
 }
 
 
-const Foam::scalarField& Foam::lduMatrix::upper() const
+Foam::scalarField& Foam::lduMatrix::lower(label nCoeffs)
 {
-    if (!lowerPtr_ && !upperPtr_)
+    if (!lowerPtr_)
     {
-        FatalErrorInFunction
-            << "lowerPtr_ or upperPtr_ unallocated"
-            << abort(FatalError);
+        if (upperPtr_)
+        {
+            lowerPtr_ = std::make_unique<scalarField>(*upperPtr_);
+        }
+        else
+        {
+            // if (nCoeffs < 0)
+            // {
+            //     nCoeffs = lduAddr().lowerAddr().size();
+            // }
+            lowerPtr_ =
+                std::make_unique<scalarField>(nCoeffs, Foam::zero{});
+        }
     }
 
-    if (upperPtr_)
-    {
-        return *upperPtr_;
-    }
-    else
-    {
-        return *lowerPtr_;
-    }
+    return *lowerPtr_;
 }
 
 
@@ -377,28 +385,25 @@ void Foam::lduMatrix::setResidualField
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum)
+Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& mat)
 {
-    Switch hasLow = ldum.hasLower();
-    Switch hasDiag = ldum.hasDiag();
-    Switch hasUp = ldum.hasUpper();
+    os  << mat.hasLower() << token::SPACE
+        << mat.hasDiag() << token::SPACE
+        << mat.hasUpper() << token::SPACE;
 
-    os  << hasLow << token::SPACE << hasDiag << token::SPACE
-        << hasUp << token::SPACE;
-
-    if (hasLow)
+    if (mat.hasLower())
     {
-        os  << ldum.lower();
+        os  << mat.lower();
     }
 
-    if (hasDiag)
+    if (mat.hasDiag())
     {
-        os  << ldum.diag();
+        os  << mat.diag();
     }
 
-    if (hasUp)
+    if (mat.hasUpper())
     {
-        os  << ldum.upper();
+        os  << mat.upper();
     }
 
     os.check(FUNCTION_NAME);
@@ -413,54 +418,50 @@ Foam::Ostream& Foam::operator<<
     const InfoProxy<lduMatrix>& iproxy
 )
 {
-    const auto& ldum = *iproxy;
-
-    Switch hasLow = ldum.hasLower();
-    Switch hasDiag = ldum.hasDiag();
-    Switch hasUp = ldum.hasUpper();
+    const auto& mat = *iproxy;
 
-    os  << "Lower:" << hasLow
-        << " Diag:" << hasDiag
-        << " Upper:" << hasUp << endl;
+    os  << "Lower:" << Switch::name(mat.hasLower())
+        << " Diag:" << Switch::name(mat.hasDiag())
+        << " Upper:" << Switch::name(mat.hasUpper()) << endl;
 
-    if (hasLow)
+    if (mat.hasLower())
     {
-        os  << "lower:" << ldum.lower().size() << endl;
+        os  << "lower:" << mat.lower().size() << endl;
     }
-    if (hasDiag)
+    if (mat.hasDiag())
     {
-        os  << "diag :" << ldum.diag().size() << endl;
+        os  << "diag :" << mat.diag().size() << endl;
     }
-    if (hasUp)
+    if (mat.hasUpper())
     {
-        os  << "upper:" << ldum.upper().size() << endl;
+        os  << "upper:" << mat.upper().size() << endl;
     }
 
 
-    //if (hasLow)
+    //if (hasLower)
     //{
     //    os  << "lower contents:" << endl;
-    //    forAll(ldum.lower(), i)
+    //    forAll(mat.lower(), i)
     //    {
-    //        os  << "i:" << i << "\t" << ldum.lower()[i] << endl;
+    //        os  << "i:" << i << "\t" << mat.lower()[i] << endl;
     //    }
     //    os  << endl;
     //}
     //if (hasDiag)
     //{
     //    os  << "diag contents:" << endl;
-    //    forAll(ldum.diag(), i)
+    //    forAll(mat.diag(), i)
     //    {
-    //        os  << "i:" << i << "\t" << ldum.diag()[i] << endl;
+    //        os  << "i:" << i << "\t" << mat.diag()[i] << endl;
     //    }
     //    os  << endl;
     //}
-    //if (hasUp)
+    //if (hasUpper)
     //{
     //    os  << "upper contents:" << endl;
-    //    forAll(ldum.upper(), i)
+    //    forAll(mat.upper(), i)
     //    {
-    //        os  << "i:" << i << "\t" << ldum.upper()[i] << endl;
+    //        os  << "i:" << i << "\t" << mat.upper()[i] << endl;
     //    }
     //    os  << endl;
     //}
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
index e13bb4cd057..a4dd1b06017 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,13 +57,14 @@ SourceFiles
 #include "primitiveFieldsFwd.H"
 #include "FieldField.H"
 #include "lduInterfaceFieldPtrsList.H"
+#include "solverPerformance.H"
 #include "typeInfo.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
-#include "solverPerformance.H"
 #include "InfoProxy.H"
 #include "Enum.H"
 #include "profilingTrigger.H"
+#include <functional>  // For reference_wrapper
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,8 +90,14 @@ class lduMatrix
         //const lduMesh& lduMesh_;
         std::reference_wrapper<const lduMesh> lduMesh_;
 
-        //- Coefficients (not including interfaces)
-        scalarField *lowerPtr_, *diagPtr_, *upperPtr_;
+        //- Diagonal coefficients
+        std::unique_ptr<scalarField> diagPtr_;
+
+        //- Off-diagonal coefficients (not including interfaces)
+        std::unique_ptr<scalarField> lowerPtr_;
+
+        //- Off-diagonal coefficients (not including interfaces)
+        std::unique_ptr<scalarField> upperPtr_;
 
 
 public:
@@ -115,6 +122,7 @@ public:
         static const scalar defaultTolerance;
 
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for lduMatrix solvers
     class solver
     {
@@ -331,6 +339,7 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for lduMatrix smoothers
     class smoother
     {
@@ -478,6 +487,7 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
     //- Abstract base-class for lduMatrix preconditioners
     class preconditioner
     {
@@ -582,6 +592,8 @@ public:
     };
 
 
+    // -----------------------------------------------------------------------
+
     // Static Data
 
         // Declare name of the class and its debug switch
@@ -590,101 +602,101 @@ public:
 
     // Constructors
 
-        //- Construct given an LDU addressed mesh.
-        //  The coefficients are initially empty for subsequent setting.
-        lduMatrix(const lduMesh&);
+        //- Construct (without coefficients) for an LDU addressed mesh.
+        //  Not yet 'explicit' (legacy code may rely on implicit construct)
+        lduMatrix(const lduMesh& mesh);
 
-        //- Construct as copy
+        //- Copy construct
         lduMatrix(const lduMatrix&);
 
+        //- Move construct
+        lduMatrix(lduMatrix&&);
+
         //- Construct as copy or re-use as specified.
         lduMatrix(lduMatrix&, bool reuse);
 
         //- Construct given an LDU addressed mesh and an Istream
         //- from which the coefficients are read
-        lduMatrix(const lduMesh&, Istream&);
+        lduMatrix(const lduMesh& mesh, Istream& is);
 
 
     //- Destructor
-    ~lduMatrix();
+    ~lduMatrix() = default;
 
 
     // Member Functions
 
-        // Access to addressing
+    // Addressing
 
-            //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const noexcept
-            {
-                return lduMesh_;
-            }
+        //- Return the LDU mesh from which the addressing is obtained
+        const lduMesh& mesh() const noexcept
+        {
+            return lduMesh_;
+        }
 
-            //- Set the LDU mesh containing the addressing is obtained
-            void setLduMesh(const lduMesh& m)
-            {
-                lduMesh_ = m;
-            }
+        //- Set the LDU mesh containing the addressing
+        void setLduMesh(const lduMesh& m)
+        {
+            lduMesh_ = m;
+        }
 
-            //- Return the LDU addressing
-            const lduAddressing& lduAddr() const
-            {
-                return mesh().lduAddr();
-            }
+        //- Return the LDU addressing
+        const lduAddressing& lduAddr() const
+        {
+            return mesh().lduAddr();
+        }
 
-            //- Return the patch evaluation schedule
-            const lduSchedule& patchSchedule() const
-            {
-                return lduAddr().patchSchedule();
-            }
+        //- Return the patch evaluation schedule
+        const lduSchedule& patchSchedule() const
+        {
+            return mesh().lduAddr().patchSchedule();
+        }
 
 
-        // Access to coefficients
+    // Coefficients
 
-            scalarField& lower();
-            scalarField& diag();
-            scalarField& upper();
+        const scalarField& diag() const;
+        const scalarField& upper() const;
+        const scalarField& lower() const;
 
-            // Size with externally provided sizes (for constructing with 'fake'
-            // mesh in GAMG)
+        scalarField& diag();
+        scalarField& upper();
+        scalarField& lower();
 
-                scalarField& lower(const label size);
-                scalarField& diag(const label nCoeffs);
-                scalarField& upper(const label nCoeffs);
+        // Size with externally provided sizes
+        // (for constructing with 'fake' mesh in GAMG)
 
+        scalarField& diag(label size);
+        scalarField& upper(label nCoeffs);
+        scalarField& lower(label nCoeffs);
 
-            const scalarField& lower() const;
-            const scalarField& diag() const;
-            const scalarField& upper() const;
 
-            bool hasDiag() const noexcept
-            {
-                return (diagPtr_);
-            }
+    // Characteristics
 
-            bool hasUpper() const noexcept
-            {
-                return (upperPtr_);
-            }
+        //- The matrix type (empty, diagonal, symmetric, ...)
+        word matrixTypeName() const;
 
-            bool hasLower() const noexcept
-            {
-                return (lowerPtr_);
-            }
+        bool hasDiag() const noexcept { return bool(diagPtr_); }
+        bool hasUpper() const noexcept { return bool(upperPtr_); }
+        bool hasLower() const noexcept { return bool(lowerPtr_); }
 
-            bool diagonal() const noexcept
-            {
-                return (diagPtr_ && !lowerPtr_ && !upperPtr_);
-            }
+        //- Matrix has diagonal only
+        bool diagonal() const noexcept
+        {
+            return (diagPtr_ && !lowerPtr_ && !upperPtr_);
+        }
 
-            bool symmetric() const noexcept
-            {
-                return (diagPtr_ && (!lowerPtr_ && upperPtr_));
-            }
+        //- Matrix is symmetric
+        bool symmetric() const noexcept
+        {
+            return (diagPtr_ && !lowerPtr_ && upperPtr_);
+        }
 
-            bool asymmetric() const noexcept
-            {
-                return (diagPtr_ && lowerPtr_ && upperPtr_);
-            }
+        //- Matrix is asymmetric (ie, full)
+        bool asymmetric() const noexcept
+        {
+            return (diagPtr_ && lowerPtr_ && upperPtr_);
+        }
 
 
         // Operations
@@ -801,8 +813,12 @@ public:
 
     // Member Operators
 
+        //- Copy assignment
         void operator=(const lduMatrix&);
 
+        //- Move assignment
+        void operator=(lduMatrix&&);
+
         void negate();
 
         void operator+=(const lduMatrix&);
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
index 4d307f57100..1e83a03c75b 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,7 +86,7 @@ void Foam::lduMatrix::sumMagOffDiag
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::lduMatrix::operator=(const lduMatrix& A)
 {
@@ -95,38 +95,49 @@ void Foam::lduMatrix::operator=(const lduMatrix& A)
         return;  // Self-assignment is a no-op
     }
 
-    if (A.lowerPtr_)
+    if (A.hasLower())
     {
         lower() = A.lower();
     }
-    else if (lowerPtr_)
+    else
     {
-        delete lowerPtr_;
-        lowerPtr_ = nullptr;
+        lowerPtr_.reset(nullptr);
     }
 
-    if (A.upperPtr_)
+    if (A.hasUpper())
     {
         upper() = A.upper();
     }
-    else if (upperPtr_)
+    else
     {
-        delete upperPtr_;
-        upperPtr_ = nullptr;
+        upperPtr_.reset(nullptr);
     }
 
-    if (A.diagPtr_)
+    if (A.hasDiag())
     {
         diag() = A.diag();
     }
 }
 
 
+void Foam::lduMatrix::operator=(lduMatrix&& A)
+{
+    if (this == &A)
+    {
+        return;  // Self-assignment is a no-op
+    }
+
+    diagPtr_ = std::move(A.diagPtr_);
+    upperPtr_ = std::move(A.upperPtr_);
+    lowerPtr_ = std::move(A.lowerPtr_);
+}
+
+
 void Foam::lduMatrix::negate()
 {
-    if (lowerPtr_)
+    if (diagPtr_)
     {
-        lowerPtr_->negate();
+        diagPtr_->negate();
     }
 
     if (upperPtr_)
@@ -134,16 +145,16 @@ void Foam::lduMatrix::negate()
         upperPtr_->negate();
     }
 
-    if (diagPtr_)
+    if (lowerPtr_)
     {
-        diagPtr_->negate();
+        lowerPtr_->negate();
     }
 }
 
 
 void Foam::lduMatrix::operator+=(const lduMatrix& A)
 {
-    if (A.diagPtr_)
+    if (A.hasDiag())
     {
         diag() += A.diag();
     }
@@ -168,7 +179,7 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
     }
     else if (asymmetric() && A.symmetric())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             lower() += A.upper();
             upper() += A.upper();
@@ -187,12 +198,12 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
     }
     else if (diagonal())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             upper() = A.upper();
         }
 
-        if (A.lowerPtr_)
+        if (A.hasLower())
         {
             lower() = A.lower();
         }
@@ -206,15 +217,8 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
         {
             WarningInFunction
                 << "Unknown matrix type combination" << nl
-                << "    this :"
-                << " diagonal:" << diagonal()
-                << " symmetric:" << symmetric()
-                << " asymmetric:" << asymmetric() << nl
-                << "    A    :"
-                << " diagonal:" << A.diagonal()
-                << " symmetric:" << A.symmetric()
-                << " asymmetric:" << A.asymmetric()
-                << endl;
+                << "    this : " << this->matrixTypeName()
+                << "    A    : " << A.matrixTypeName() << endl;
         }
     }
 }
@@ -247,7 +251,7 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A)
     }
     else if (asymmetric() && A.symmetric())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             lower() -= A.upper();
             upper() -= A.upper();
@@ -266,12 +270,12 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A)
     }
     else if (diagonal())
     {
-        if (A.upperPtr_)
+        if (A.hasUpper())
         {
             upper() = -A.upper();
         }
 
-        if (A.lowerPtr_)
+        if (A.hasLower())
         {
             lower() = -A.lower();
         }
@@ -285,15 +289,8 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A)
         {
             WarningInFunction
                 << "Unknown matrix type combination" << nl
-                << "    this :"
-                << " diagonal:" << diagonal()
-                << " symmetric:" << symmetric()
-                << " asymmetric:" << asymmetric() << nl
-                << "    A    :"
-                << " diagonal:" << A.diagonal()
-                << " symmetric:" << A.symmetric()
-                << " asymmetric:" << A.asymmetric()
-                << endl;
+                << "    this : " << this->matrixTypeName()
+                << "    A    : " << A.matrixTypeName() << endl;
         }
     }
 }
-- 
GitLab