Feature single precision solve type

applications/test/field1/Test-field1.C

@@ -24,6 +24,10 @@ License
 Description
     Simple field tests
 
+    Test use of Kahan/Neumaier to extend precision for when running SPDP
+    mode. Conclusion is that it is easier/quicker to run these summation
+    loops as double precision (i.e. solveScalar).
+
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveFields.H"
@@ -32,6 +36,116 @@ Description
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class CombineOp, class ResultType>
+void sumNeumaier
+(
+    const UList<Type>& vals,
+    const CombineOp& cop,
+    ResultType& result
+)
+{
+    // Neumaier version of Kahan
+    ResultType sum = Zero;
+    ResultType c = Zero;
+    for (const Type& vali : vals)
+    {
+        ResultType val;
+        cop(val, vali);
+
+        const ResultType t = sum + val;
+
+        for
+        (
+            direction cmpt = 0;
+            cmpt < pTraits<ResultType>::nComponents;
+            cmpt++
+        )
+        {
+            if (mag(sum[cmpt]) >= mag(val[cmpt]))
+            {
+                // If sum is bigger, low-order digits of input[i] are lost.
+                c[cmpt] += (sum[cmpt] - t[cmpt]) + val[cmpt];
+            }
+            else
+            {
+                // Else low-order digits of sum are lost.
+                c[cmpt] += (val[cmpt] - t[cmpt]) + sum[cmpt];
+            }
+        }
+        sum = t;
+    }
+    result = sum + c;
+}
+
+
+template<class CombineOp, class ResultType>
+void sumNeumaier
+(
+    const UList<scalar>& vals,
+    const CombineOp& cop,
+    ResultType& result
+)
+{
+    // Neumaier version of Kahan
+    ResultType sum = Zero;
+    ResultType c = Zero;
+    for (const scalar vali : vals)
+    {
+        ResultType val;
+        cop(val, vali);
+
+        const ResultType t = sum + val;
+        if (mag(sum) >= mag(val))
+        {
+            // If sum is bigger, low-order digits of input[i] are lost.
+            c += (sum - t) + val;
+        }
+        else
+        {
+            // Else low-order digits of sum are lost.
+            c += (val - t) + sum;
+        }
+        sum = t;
+    }
+    result = sum + c;
+}
+
+
+template<class Type>
+Type mySum(const UList<Type>& f)
+{
+    typedef typename Foam::typeOfSolve<Type>::type solveType;
+
+    solveType Sum = Zero;
+
+    if (f.size())
+    {
+        sumNeumaier(f, eqOp<solveType>(), Sum);
+    }
+
+    return Type(Sum);
+}
+
+
+//- The sumSqr always adds only positive numbers. Here there can never be any
+//  cancellation of truncation errors.
+template<class Type>
+typename outerProduct1<Type>::type mySumSqr(const UList<Type>& f)
+{
+    typedef typename outerProduct1<solveScalar>::type prodType;
+
+    prodType result = Zero;
+
+    if (f.size())
+    {
+        sumNeumaier(f, eqSqrOp<prodType>(), result);
+    }
+
+    return result;
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -63,6 +177,71 @@ int main(int argc, char *argv[])
 
     Info<< "negated: " << lfield << nl;
 
+
+    // Summation (compile in SPDP)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Pout.precision(16);
+    Sout.precision(16);
+
+    const scalar SMALLS(1e-6);
+    const vector SMALLV(SMALLS, SMALLS, SMALLS);
+    const scalar GREATS(1e6);
+    const vector GREATV(GREATS, GREATS, GREATS);
+
+    // scalarField summation
+    {
+        scalarField sfield(10, SMALLS);
+
+        sfield[8] = GREATS;
+        sfield[9] = -sfield[8];
+        Info<< "scalarField:" << sfield.size() << nl
+            << "    sum      :" << sum(sfield) << nl
+            << "    corrected:" << mySum(sfield) << endl;
+    }
+    // vectorField summation
+    {
+        vectorField vfield(10, SMALLV);
+
+        vfield[8] = GREATV;
+        vfield[9] = -vfield[8];
+        Info<< "vectorField:" << vfield.size() << nl
+            << "    sum      :" << sum(vfield) << nl
+            << "    corrected:" << mySum(vfield) << endl;
+    }
+    // sphericalTensorField summation
+    {
+        sphericalTensorField tfield(10, SMALLS);
+
+        tfield[8] = GREATS;
+        tfield[9] = -tfield[8];
+        Info<< "sphericalTensorField:" << tfield.size() << nl
+            << "    sum      :" << sum(tfield) << nl
+            << "    corrected:" << mySum(tfield) << endl;
+    }
+    // symmTensorField summation
+    {
+        symmTensorField tfield(10, SMALLS*symmTensor::I);
+
+        tfield[8] = GREATS*symmTensor::I;
+        tfield[9] = -tfield[8];
+        Info<< "symmTensorField:" << tfield.size() << nl
+            << "    sum      :" << sum(tfield) << nl
+            << "    corrected:" << mySum(tfield) << endl;
+    }
+    // tensorField summation
+    {
+        tensorField tfield(10, SMALLS*tensor::I);
+
+        tfield[8] = GREATS*tensor::I;
+        tfield[9] = -tfield[8];
+        Info<< "tensorField:" << tfield.size() << nl
+            << "    sum      :" << sum(tfield) << nl
+            << "    corrected:" << mySum(tfield) << endl;
+    }
+
+
+
     return 0;
 }
 

src/OpenFOAM/fields/Fields/Field/FieldFunctions.C

@@ -390,14 +390,16 @@ TMP_UNARY_FUNCTION(Type, min)
 template<class Type>
 Type sum(const UList<Type>& f)
 {
-    Type Sum = Zero;
+    typedef typename Foam::typeOfSolve<Type>::type solveType;
+
+    solveType Sum = Zero;
 
     if (f.size())
     {
-        TFOR_ALL_S_OP_F(Type, Sum, +=, Type, f)
+        TFOR_ALL_S_OP_FUNC_F(solveType, Sum, +=, solveType, Type, f)
     }
 
-    return Sum;
+    return Type(Sum);
 }
 
 TMP_UNARY_FUNCTION(Type, sum)

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
+#include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -75,10 +76,17 @@ void Foam::primitiveMesh::makeCellCentresAndVols
 (
     const vectorField& fCtrs,
     const vectorField& fAreas,
-    vectorField& cellCtrs,
+    vectorField& cellCtrs_s,
-    scalarField& cellVols
+    scalarField& cellVols_s
 ) const
 {
+    typedef Vector<solveScalar> solveVector;
+
+    PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s);
+    Field<solveVector>& cellCtrs = tcellCtrs.ref();
+    PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s);
+    Field<solveScalar>& cellVols = tcellVols.ref();
+
     // Clear the fields for accumulation
     cellCtrs = Zero;
     cellVols = 0.0;
@@ -89,18 +97,18 @@ void Foam::primitiveMesh::makeCellCentresAndVols
     // first estimate the approximate cell centre as the average of
     // face centres
 
-    vectorField cEst(nCells(), Zero);
+    Field<solveVector> cEst(nCells(), Zero);
     labelField nCellFaces(nCells(), Zero);
 
     forAll(own, facei)
     {
-        cEst[own[facei]] += fCtrs[facei];
+        cEst[own[facei]] += solveVector(fCtrs[facei]);
         ++nCellFaces[own[facei]];
     }
 
     forAll(nei, facei)
     {
-        cEst[nei[facei]] += fCtrs[facei];
+        cEst[nei[facei]] += solveVector(fCtrs[facei]);
         ++nCellFaces[nei[facei]];
     }
 
@@ -111,12 +119,15 @@ void Foam::primitiveMesh::makeCellCentresAndVols
 
     forAll(own, facei)
     {
+        const solveVector fc(fCtrs[facei]);
+        const solveVector fA(fAreas[facei]);
+
         // Calculate 3*face-pyramid volume
-        scalar pyr3Vol =
+        solveScalar pyr3Vol =
-            fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
+            fA & (fc - cEst[own[facei]]);
 
         // Calculate face-pyramid centre
-        vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
+        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
 
         // Accumulate volume-weighted face-pyramid centre
         cellCtrs[own[facei]] += pyr3Vol*pc;
@@ -127,12 +138,15 @@ void Foam::primitiveMesh::makeCellCentresAndVols
 
     forAll(nei, facei)
     {
+        const solveVector fc(fCtrs[facei]);
+        const solveVector fA(fAreas[facei]);
+
         // Calculate 3*face-pyramid volume
-        scalar pyr3Vol =
+        solveScalar pyr3Vol =
-            fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
+            fA & (cEst[nei[facei]] - fc);
 
         // Calculate face-pyramid centre
-        vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
+        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
 
         // Accumulate volume-weighted face-pyramid centre
         cellCtrs[nei[facei]] += pyr3Vol*pc;

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -61,22 +61,22 @@ bool Foam::primitiveMesh::checkClosedBoundary
     // Loop through all boundary faces and sum up the face area vectors.
     // For a closed boundary, this should be zero in all vector components
 
-    vector sumClosed(Zero);
+    Vector<solveScalar> sumClosed(Zero);
-    scalar sumMagClosedBoundary = 0;
+    solveScalar sumMagClosedBoundary = 0;
 
     for (label facei = nInternalFaces(); facei < areas.size(); facei++)
     {
         if (!internalOrCoupledFaces.size() || !internalOrCoupledFaces[facei])
         {
-            sumClosed += areas[facei];
+            sumClosed += Vector<solveScalar>(areas[facei]);
             sumMagClosedBoundary += mag(areas[facei]);
         }
     }
 
-    reduce(sumClosed, sumOp<vector>());
+    reduce(sumClosed, sumOp<Vector<solveScalar>>());
-    reduce(sumMagClosedBoundary, sumOp<scalar>());
+    reduce(sumMagClosedBoundary, sumOp<solveScalar>());
 
-    vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
+    Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
 
     if (cmptMax(cmptMag(openness)) > closedThreshold_)
     {

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -420,6 +420,7 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
     tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
     scalarField& faceFlatness = tfaceFlatness.ref();
 
+    typedef Vector<solveScalar> solveVector;
 
     forAll(fcs, facei)
     {
@@ -427,20 +428,20 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
 
         if (f.size() > 3 && magAreas[facei] > ROOTVSMALL)
         {
-            const point& fc = fCtrs[facei];
+            const solveVector fc = fCtrs[facei];
 
             // Calculate the sum of magnitude of areas and compare to magnitude
             // of sum of areas.
 
-            scalar sumA = 0.0;
+            solveScalar sumA = 0.0;
 
             forAll(f, fp)
             {
-                const point& thisPoint = p[f[fp]];
+                const solveVector thisPoint = p[f[fp]];
-                const point& nextPoint = p[f.nextLabel(fp)];
+                const solveVector nextPoint = p[f.nextLabel(fp)];
 
                 // Triangle around fc.
-                vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
+                solveVector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
                 sumA += mag(n);
             }
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C

@@ -84,7 +84,7 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
     forAll(fs, facei)
     {
         const labelList& f = fs[facei];
-        label nPoints = f.size();
+        const label nPoints = f.size();
 
         // If the face is a triangle, do a direct calculation for efficiency
         // and to avoid round-off error-related problems
@@ -95,26 +95,29 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
         }
         else
         {
-            vector sumN = Zero;
+            typedef Vector<solveScalar> solveVector;
-            scalar sumA = 0.0;
-            vector sumAc = Zero;
 
-            point fCentre = p[f[0]];
+            solveVector sumN = Zero;
+            solveScalar sumA = 0.0;
+            solveVector sumAc = Zero;
+
+            solveVector fCentre = p[f[0]];
             for (label pi = 1; pi < nPoints; pi++)
             {
-                fCentre += p[f[pi]];
+                fCentre += solveVector(p[f[pi]]);
             }
 
             fCentre /= nPoints;
 
             for (label pi = 0; pi < nPoints; pi++)
             {
-                const point& nextPoint = p[f[(pi + 1) % nPoints]];
+                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
-
+                const solveVector nextPoint(p[f[nextPi]]);
-                vector c = p[f[pi]] + nextPoint + fCentre;
+                const solveVector thisPoint(p[f[pi]]);
-                vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
-                scalar a = mag(n);
 
+                solveVector c = thisPoint + nextPoint + fCentre;
+                solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
+                solveScalar a = mag(n);
                 sumN += n;
                 sumA += a;
                 sumAc += a*c;

src/OpenFOAM/primitives/Scalar/scalar/scalar.H

@@ -154,7 +154,7 @@ namespace Foam
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Additional transcendental functions
+// Additional transcendental functions and specialisations
 
 namespace Foam
 {
@@ -172,6 +172,15 @@ namespace Foam
 
     //- Lower incomplete gamma function
     scalar incGamma_P(const scalar a, const scalar x);
+
+    //- Type to use for extended precision
+    template<>
+    class typeOfSolve<scalar>
+    {
+    public:
+
+        typedef solveScalar type;
+    };
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H

@@ -137,6 +137,15 @@ struct is_contiguous_scalar<SphericalTensor<Cmpt>>
 {};
 
 
+template<class Cmpt>
+class typeOfSolve<SphericalTensor<Cmpt>>
+{
+public:
+
+    typedef SphericalTensor<solveScalar> type;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/SymmTensor/SymmTensor.H

@@ -176,6 +176,15 @@ public:
 };
 
 
+template<class Cmpt>
+class typeOfSolve<SymmTensor<Cmpt>>
+{
+public:
+
+    typedef SymmTensor<solveScalar> type;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/Tensor/Tensor.H

@@ -344,6 +344,15 @@ public:
 };
 
 
+template<class Cmpt>
+class typeOfSolve<Tensor<Cmpt>>
+{
+public:
+
+    typedef Tensor<solveScalar> type;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/Vector/Vector.H

@@ -167,6 +167,15 @@ public:
 };
 
 
+template<class Cmpt>
+class typeOfSolve<Vector<Cmpt>>
+{
+public:
+
+    typedef Vector<solveScalar> type;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/VectorSpace/products.H

@@ -35,6 +35,7 @@ Description
 #ifndef products_H
 #define products_H
 
+#include "direction.H"
 #include "pTraits.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,6 +73,16 @@ class symmTypeOfRank
 {};
 
 
+//- The extended precision type (solveScalar for float)
+template<class Type>
+class typeOfSolve
+{
+public:
+
+    typedef Type type;
+};
+
+
 //- The magnitude type for given argument.
 template<class arg1>
 class typeOfMag