diff --git a/applications/test/field1/Test-field1.C b/applications/test/field1/Test-field1.C
index c2b739e398676ebcb43dcd497366c4e44d514eec..3fdfcb6a7ff429d21435a5615f378e53071f8d92 100644
--- a/applications/test/field1/Test-field1.C
+++ b/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;
}
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
index 21badff3c815d883ef13ea5c528f4f5742bff5bc..0c5091fab6eb8e04a7671b56868cecb95660b019 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
+++ b/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)
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C
index d71ae13909507ed55167b2ba5edf72891bcae738..c7d4b1bd05547a89d9a88fd20c65d78d811de989 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C
+++ b/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,
- scalarField& cellVols
+ vectorField& cellCtrs_s,
+ 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 =
- fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
+ solveScalar pyr3Vol =
+ 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 =
- fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
+ solveScalar pyr3Vol =
+ 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;
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 05a56e62d4c8db12bb324a89aa5123a3adce9736..4f4a1445f37c719769a6bf878fd90283dfea6751 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/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);
- scalar sumMagClosedBoundary = 0;
+ Vector<solveScalar> sumClosed(Zero);
+ 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(sumMagClosedBoundary, sumOp<scalar>());
+ reduce(sumClosed, sumOp<Vector<solveScalar>>());
+ reduce(sumMagClosedBoundary, sumOp<solveScalar>());
- vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
+ Vector<solveScalar> openness = sumClosed/(sumMagClosedBoundary + VSMALL);
if (cmptMax(cmptMag(openness)) > closedThreshold_)
{
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
index 95db7415385f9d71cfef1f8742ad45c5f03847f5..98a8b032b7332216ce0bfd939ba8b8fa4bcf9157 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
+++ b/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 point& nextPoint = p[f.nextLabel(fp)];
+ const solveVector thisPoint = p[f[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);
}
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
index dfd90c25ebe01b562d1db3116edc75150ad28f8d..ac906072acf557caf7e572221ffd844466b41904 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
+++ b/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;
- scalar sumA = 0.0;
- vector sumAc = Zero;
+ typedef Vector<solveScalar> solveVector;
- 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]];
-
- vector c = p[f[pi]] + nextPoint + fCentre;
- vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
- scalar a = mag(n);
+ const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+ const solveVector nextPoint(p[f[nextPi]]);
+ const solveVector thisPoint(p[f[pi]]);
+ solveVector c = thisPoint + nextPoint + fCentre;
+ solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
+ solveScalar a = mag(n);
sumN += n;
sumA += a;
sumAc += a*c;
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H b/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
index e6888e71eb78ccb835c3ed9f563ffe9164f2b00b..25359893165017312b98436491623de9959d1c98 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
+++ b/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;
+ };
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H b/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H
index 47d895d54b4e264a4c19a1606d9106ee4069bf8f..3ac5117a5d84711b58b8e74aa108eac5e63cb2e5 100644
--- a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H
+++ b/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
diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
index 37e8bf8c98b02a3ddd34f0c9605a9d5aa8a4a220..998c6b5b3394673975cfb44d31e27bbbb09daa54 100644
--- a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
+++ b/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
diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H
index 561232f630bb63f2b75401e7ad26d6cc604d0f77..253bcc00d587b8cb44d418ce3eac4dc7466daabb 100644
--- a/src/OpenFOAM/primitives/Tensor/Tensor.H
+++ b/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
diff --git a/src/OpenFOAM/primitives/Vector/Vector.H b/src/OpenFOAM/primitives/Vector/Vector.H
index f4fc10316aa7bde5f1abc82ae949968231083dac..100a6f704da71733fb055bc5dc285db06d6ca397 100644
--- a/src/OpenFOAM/primitives/Vector/Vector.H
+++ b/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
diff --git a/src/OpenFOAM/primitives/VectorSpace/products.H b/src/OpenFOAM/primitives/VectorSpace/products.H
index e0643885cce7586f038ba1409242c24c7e0d5214..fcf9f4a46494ee17fbb845918bf0493e2856b8b5 100644
--- a/src/OpenFOAM/primitives/VectorSpace/products.H
+++ b/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