ENH: mixed precision solver infrastructure (#1086)

- lduSolver: solveScalar for coarse-level solve; normFactor etc

- lduMatrix: have smoother natively use solveScalarField

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2012 OpenFOAM Foundation
@@ -46,7 +46,7 @@ Foam::cyclicLduInterfaceField::~cyclicLduInterfaceField()
 
 void Foam::cyclicLduInterfaceField::transformCoupleField
 (
-    scalarField& f,
+    solveScalarField& f,
     const direction cmpt
 ) const
 {

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2012 OpenFOAM Foundation
@@ -93,7 +93,7 @@ public:
         //- Transform given patch component field
         void transformCoupleField
         (
-            scalarField& f,
+            solveScalarField& f,
             const direction cmpt
         ) const;
 };

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.C

@@ -35,10 +35,4 @@ defineTypeNameAndDebug(lduInterfaceField, 0);
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::lduInterfaceField::~lduInterfaceField()
-{}
-
-
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2013 OpenFOAM Foundation
@@ -93,7 +93,7 @@ public:
 
 
     //- Destructor
-    virtual ~lduInterfaceField();
+    virtual ~lduInterfaceField() = default;
 
 
     // Member Functions
@@ -137,11 +137,11 @@ public:
             //  or subtract coupled contributions to matrix
             virtual void initInterfaceMatrixUpdate
             (
-                scalarField&,
+                solveScalarField& result,
                 const bool add,
-                const scalarField&,
-                const scalarField&,
-                const direction,
+                const solveScalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const
             {}
@@ -150,11 +150,11 @@ public:
             //  or subtract coupled contributions to matrix
             virtual void updateInterfaceMatrix
             (
-                scalarField&,
+                solveScalarField& result,
                 const bool add,
-                const scalarField&,
-                const scalarField&,
-                const direction,
+                const solveScalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const = 0;
 

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/processorLduInterfaceField/processorLduInterfaceField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2012 OpenFOAM Foundation
@@ -46,7 +46,7 @@ Foam::processorLduInterfaceField::~processorLduInterfaceField()
 
 void Foam::processorLduInterfaceField::transformCoupleField
 (
-    scalarField& f,
+    solveScalarField& f,
     const direction cmpt
 ) const
 {

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/processorLduInterfaceField/processorLduInterfaceField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2014 OpenFOAM Foundation
@@ -99,7 +99,7 @@ public:
         //- Transform given patch component field
         void transformCoupleField
         (
-            scalarField& f,
+            solveScalarField& f,
             const direction cmpt
         ) const;
 };

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -29,7 +29,7 @@ License
 #include "IOstreams.H"
 #include "Switch.H"
 #include "objectRegistry.H"
-#include "IOField.H"
+#include "scalarIOField.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -319,7 +319,7 @@ const Foam::scalarField& Foam::lduMatrix::upper() const
 
 void Foam::lduMatrix::setResidualField
 (
-    const Field<scalar>& residual,
+    const scalarField& residual,
     const word& fieldName,
     const bool initial
 ) const
@@ -339,8 +339,8 @@ void Foam::lduMatrix::setResidualField
         lookupName = word("residual:" + fieldName);
     }
 
-    IOField<scalar>* residualPtr =
-        lduMesh_.thisDb().getObjectPtr<IOField<scalar>>(lookupName);
+    scalarIOField* residualPtr =
+        lduMesh_.thisDb().getObjectPtr<scalarIOField>(lookupName);
 
     if (residualPtr)
     {

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -262,12 +262,12 @@ public:
 
             //- Return the matrix norm used to normalise the residual for the
             //- stopping criterion
-            scalar normFactor
+            solveScalarField::cmptType normFactor
             (
-                const scalarField& psi,
-                const scalarField& source,
-                const scalarField& Apsi,
-                scalarField& tmpField
+                const solveScalarField& psi,
+                const solveScalarField& source,
+                const solveScalarField& Apsi,
+                solveScalarField& tmpField
             ) const;
     };
 
@@ -404,11 +404,20 @@ public:
             //- Smooth the solution for a given number of sweeps
             virtual void smooth
             (
-                scalarField& psi,
+                solveScalarField& psi,
                 const scalarField& source,
                 const direction cmpt,
                 const label nSweeps
             ) const = 0;
+
+            //- Smooth the solution for a given number of sweeps
+            virtual void scalarSmooth
+            (
+                solveScalarField& psi,
+                const solveScalarField& source,
+                const direction cmpt,
+                const label nSweeps
+            ) const = 0;
     };
 
 
@@ -495,8 +504,8 @@ public:
             //- Return wA the preconditioned form of residual rA
             virtual void precondition
             (
-                scalarField& wA,
-                const scalarField& rA,
+                solveScalarField& wA,
+                const solveScalarField& rA,
                 const direction cmpt=0
             ) const = 0;
 
@@ -505,8 +514,8 @@ public:
             //  This is only required for preconditioning asymmetric matrices.
             virtual void preconditionT
             (
-                scalarField& wT,
-                const scalarField& rT,
+                solveScalarField& wT,
+                const solveScalarField& rT,
                 const direction cmpt=0
             ) const
             {
@@ -624,8 +633,8 @@ public:
             //- Matrix multiplication with updated interfaces.
             void Amul
             (
-                scalarField&,
-                const tmp<scalarField>&,
+                solveScalarField&,
+                const tmp<solveScalarField>&,
                 const FieldField<Field, scalar>&,
                 const lduInterfaceFieldPtrsList&,
                 const direction cmpt
@@ -634,8 +643,8 @@ public:
             //- Matrix transpose multiplication with updated interfaces.
             void Tmul
             (
-                scalarField&,
-                const tmp<scalarField>&,
+                solveScalarField&,
+                const tmp<solveScalarField>&,
                 const FieldField<Field, scalar>&,
                 const lduInterfaceFieldPtrsList&,
                 const direction cmpt
@@ -645,7 +654,7 @@ public:
             //- Sum the coefficients on each row of the matrix
             void sumA
             (
-                scalarField&,
+                solveScalarField&,
                 const FieldField<Field, scalar>&,
                 const lduInterfaceFieldPtrsList&
             ) const;
@@ -653,17 +662,17 @@ public:
 
             void residual
             (
-                scalarField& rA,
-                const scalarField& psi,
+                solveScalarField& rA,
+                const solveScalarField& psi,
                 const scalarField& source,
                 const FieldField<Field, scalar>& interfaceBouCoeffs,
                 const lduInterfaceFieldPtrsList& interfaces,
                 const direction cmpt
             ) const;
 
-            tmp<scalarField> residual
+            tmp<solveScalarField> residual
             (
-                const scalarField& psi,
+                const solveScalarField& psi,
                 const scalarField& source,
                 const FieldField<Field, scalar>& interfaceBouCoeffs,
                 const lduInterfaceFieldPtrsList& interfaces,
@@ -678,8 +687,8 @@ public:
                 const bool add,
                 const FieldField<Field, scalar>& interfaceCoeffs,
                 const lduInterfaceFieldPtrsList& interfaces,
-                const scalarField& psiif,
-                scalarField& result,
+                const solveScalarField& psiif,
+                solveScalarField& result,
                 const direction cmpt
             ) const;
 
@@ -689,8 +698,8 @@ public:
                 const bool add,
                 const FieldField<Field, scalar>& interfaceCoeffs,
                 const lduInterfaceFieldPtrsList& interfaces,
-                const scalarField& psiif,
-                scalarField& result,
+                const solveScalarField& psiif,
+                solveScalarField& result,
                 const direction cmpt
             ) const;
 
@@ -698,7 +707,7 @@ public:
             //- if it exists
             void setResidualField
             (
-                const Field<scalar>& residual,
+                const scalarField& residual,
                 const word& fieldName,
                 const bool initial
             ) const;

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -35,17 +35,17 @@ Description
 
 void Foam::lduMatrix::Amul
 (
-    scalarField& Apsi,
-    const tmp<scalarField>& tpsi,
+    solveScalarField& Apsi,
+    const tmp<solveScalarField>& tpsi,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
     const direction cmpt
 ) const
 {
-    scalar* __restrict__ ApsiPtr = Apsi.begin();
+    solveScalar* __restrict__ ApsiPtr = Apsi.begin();
 
-    const scalarField& psi = tpsi();
-    const scalar* const __restrict__ psiPtr = psi.begin();
+    const solveScalarField& psi = tpsi();
+    const solveScalar* const __restrict__ psiPtr = psi.begin();
 
     const scalar* const __restrict__ diagPtr = diag().begin();
 
@@ -98,17 +98,17 @@ void Foam::lduMatrix::Amul
 
 void Foam::lduMatrix::Tmul
 (
-    scalarField& Tpsi,
-    const tmp<scalarField>& tpsi,
+    solveScalarField& Tpsi,
+    const tmp<solveScalarField>& tpsi,
     const FieldField<Field, scalar>& interfaceIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
     const direction cmpt
 ) const
 {
-    scalar* __restrict__ TpsiPtr = Tpsi.begin();
+    solveScalar* __restrict__ TpsiPtr = Tpsi.begin();
 
-    const scalarField& psi = tpsi();
-    const scalar* const __restrict__ psiPtr = psi.begin();
+    const solveScalarField& psi = tpsi();
+    const solveScalar* const __restrict__ psiPtr = psi.begin();
 
     const scalar* const __restrict__ diagPtr = diag().begin();
 
@@ -159,12 +159,12 @@ void Foam::lduMatrix::Tmul
 
 void Foam::lduMatrix::sumA
 (
-    scalarField& sumA,
+    solveScalarField& sumA,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces
 ) const
 {
-    scalar* __restrict__ sumAPtr = sumA.begin();
+    solveScalar* __restrict__ sumAPtr = sumA.begin();
 
     const scalar* __restrict__ diagPtr = diag().begin();
 
@@ -208,17 +208,17 @@ void Foam::lduMatrix::sumA
 
 void Foam::lduMatrix::residual
 (
-    scalarField& rA,
-    const scalarField& psi,
+    solveScalarField& rA,
+    const solveScalarField& psi,
     const scalarField& source,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
     const direction cmpt
 ) const
 {
-    scalar* __restrict__ rAPtr = rA.begin();
+    solveScalar* __restrict__ rAPtr = rA.begin();
 
-    const scalar* const __restrict__ psiPtr = psi.begin();
+    const solveScalar* const __restrict__ psiPtr = psi.begin();
     const scalar* const __restrict__ diagPtr = diag().begin();
     const scalar* const __restrict__ sourcePtr = source.begin();
 
@@ -278,16 +278,16 @@ void Foam::lduMatrix::residual
 }
 
 
-Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual
+Foam::tmp<Foam::Field<Foam::solveScalar>> Foam::lduMatrix::residual
 (
-    const scalarField& psi,
+    const solveScalarField& psi,
     const scalarField& source,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
     const direction cmpt
 ) const
 {
-    tmp<scalarField> trA(new scalarField(psi.size()));
+    tmp<solveScalarField> trA(new solveScalarField(psi.size()));
     residual(trA.ref(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
     return trA;
 }

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -170,12 +170,12 @@ void Foam::lduMatrix::solver::read(const dictionary& solverControls)
 }
 
 
-Foam::scalar Foam::lduMatrix::solver::normFactor
+Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
 (
-    const scalarField& psi,
-    const scalarField& source,
-    const scalarField& Apsi,
-    scalarField& tmpField
+    const solveScalarField& psi,
+    const solveScalarField& source,
+    const solveScalarField& Apsi,
+    solveScalarField& tmpField
 ) const
 {
     // --- Calculate A dot reference value of psi

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -34,8 +34,8 @@ void Foam::lduMatrix::initMatrixInterfaces
     const bool add,
     const FieldField<Field, scalar>& coupleCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const scalarField& psiif,
-    scalarField& result,
+    const solveScalarField& psiif,
+    solveScalarField& result,
     const direction cmpt
 ) const
 {
@@ -103,8 +103,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
     const bool add,
     const FieldField<Field, scalar>& coupleCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const scalarField& psiif,
-    scalarField& result,
+    const solveScalarField& psiif,
+    solveScalarField& result,
     const direction cmpt
 ) const
 {

src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2015 OpenFOAM Foundation
@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "DICPreconditioner.H"
+#include <algorithm>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -48,8 +49,11 @@ Foam::DICPreconditioner::DICPreconditioner
 )
 :
     lduMatrix::preconditioner(sol),
-    rD_(sol.matrix().diag())
+    rD_(sol.matrix().diag().size())
 {
+    const scalarField& diag = sol.matrix().diag();
+    std::copy(diag.begin(), diag.end(), rD_.begin());
+
     calcReciprocalD(rD_, sol.matrix());
 }
 
@@ -58,11 +62,11 @@ Foam::DICPreconditioner::DICPreconditioner
 
 void Foam::DICPreconditioner::calcReciprocalD
 (
-    scalarField& rD,
+    solveScalarField& rD,
     const lduMatrix& matrix
 )
 {
-    scalar* __restrict__ rDPtr = rD.begin();
+    solveScalar* __restrict__ rDPtr = rD.begin();
 
     const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
     const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
@@ -88,14 +92,14 @@ void Foam::DICPreconditioner::calcReciprocalD
 
 void Foam::DICPreconditioner::precondition
 (
-    scalarField& wA,
-    const scalarField& rA,
+    solveScalarField& wA,
+    const solveScalarField& rA,
     const direction
 ) const
 {
-    scalar* __restrict__ wAPtr = wA.begin();
-    const scalar* __restrict__ rAPtr = rA.begin();
-    const scalar* __restrict__ rDPtr = rD_.begin();
+    solveScalar* __restrict__ wAPtr = wA.begin();
+    const solveScalar* __restrict__ rAPtr = rA.begin();
+    const solveScalar* __restrict__ rDPtr = rD_.begin();
 
     const label* const __restrict__ uPtr =
         solver_.matrix().lduAddr().upperAddr().begin();
@@ -104,9 +108,9 @@ void Foam::DICPreconditioner::precondition
     const scalar* const __restrict__ upperPtr =
         solver_.matrix().upper().begin();
 
-    label nCells = wA.size();
-    label nFaces = solver_.matrix().upper().size();
-    label nFacesM1 = nFaces - 1;
+    const label nCells = wA.size();
+    const label nFaces = solver_.matrix().upper().size();
+    const label nFacesM1 = nFaces - 1;
 
     for (label cell=0; cell<nCells; cell++)
     {

src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011 OpenFOAM Foundation
@@ -60,7 +60,7 @@ class DICPreconditioner
     // Private data
 
         //- The reciprocal preconditioned diagonal
-        scalarField rD_;
+        solveScalarField rD_;
 
 
 public:
@@ -87,13 +87,13 @@ public:
     // Member Functions