ENH: add primitives support for mixed precision (#1086)

- add vsmall pTraits for scalars
- report the solve scalar in buildArch information

etc/bashrc

@@ -81,7 +81,7 @@ export WM_COMPILER=Gcc
 export WM_ARCH_OPTION=64
 
 # [WM_PRECISION_OPTION] - Floating-point precision:
-# = DP | SP
+# = DP | SP | SPDP
 export WM_PRECISION_OPTION=DP
 
 # [WM_LABEL_SIZE] - Label size in bits:

etc/config.csh/aliases

@@ -25,8 +25,9 @@
 alias wmSet 'source $WM_PROJECT_DIR/etc/cshrc'
 alias wmInt32 'wmSet WM_LABEL_SIZE=32'
 alias wmInt64 'wmSet WM_LABEL_SIZE=64'
-alias wmSP 'wmSet WM_PRECISION_OPTION=SP'
 alias wmDP 'wmSet WM_PRECISION_OPTION=DP'
+alias wmSP 'wmSet WM_PRECISION_OPTION=SP'
+alias wmSPDP 'wmSet WM_PRECISION_OPTION=SPDP'
 
 # Clear env
 alias wmUnset 'source $WM_PROJECT_DIR/etc/config.csh/unset'

etc/config.csh/unset

@@ -162,8 +162,9 @@ endif
 unalias wmSet
 unalias wmInt32
 unalias wmInt64
-unalias wmSP
 unalias wmDP
+unalias wmSP
+unalias wmSPDP
 
 unalias wmUnset
 

etc/config.sh/aliases

@@ -25,8 +25,9 @@
 alias wmSet='. $WM_PROJECT_DIR/etc/bashrc'
 alias wmInt32='wmSet WM_LABEL_SIZE=32'
 alias wmInt64='wmSet WM_LABEL_SIZE=64'
-alias wmSP='wmSet WM_PRECISION_OPTION=SP'
 alias wmDP='wmSet WM_PRECISION_OPTION=DP'
+alias wmSP='wmSet WM_PRECISION_OPTION=SP'
+alias wmSPDP='wmSet WM_PRECISION_OPTION=SPDP'
 
 # Clear env
 alias wmUnset='. $WM_PROJECT_DIR/etc/config.sh/unset'

etc/config.sh/unset

@@ -156,8 +156,9 @@ fi
 unalias wmSet       2>/dev/null
 unalias wmInt32     2>/dev/null
 unalias wmInt64     2>/dev/null
-unalias wmSP        2>/dev/null
 unalias wmDP        2>/dev/null
+unalias wmSP        2>/dev/null
+unalias wmSPDP      2>/dev/null
 
 unalias wmUnset     2>/dev/null
 

etc/cshrc

@@ -83,7 +83,7 @@ setenv WM_COMPILER Gcc
 setenv WM_ARCH_OPTION 64
 
 # [WM_PRECISION_OPTION] - Floating-point precision:
-# = DP | SP
+# = DP | SP | SPDP
 setenv WM_PRECISION_OPTION DP
 
 # [WM_LABEL_SIZE] - Label size in bits:

src/OpenFOAM/Make/files

@@ -79,7 +79,7 @@ primitives/Vector/complexVector/complexVector.C
 primitives/Vector/doubleVector/doubleVector.C
 primitives/Tensor/doubleTensor/doubleTensor.C
 #endif
-#if !defined(WM_SP)
+#if !defined(WM_SP) && !defined(WM_SPDP)
 primitives/Vector/floatVector/floatVector.C
 primitives/Tensor/floatTensor/floatTensor.C
 #endif

src/OpenFOAM/fields/Fields/Field/PrecisionAdaptor/PrecisionAdaptor.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PrecisionAdaptor
+
+Description
+    Conversion adaptor for Field that either wraps as a tmp reference
+    or creates the necessary tmp and copies the values on construction
+    and destruction.
+    This provides automatic conversion between (scalar) types for use
+    with linear solvers able to run mixed precision.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PrecisionAdaptor_H
+#define PrecisionAdaptor_H
+
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- A const Field wrapper with possible data conversion
+template<class Type, class InputType>
+class ConstPrecisionAdaptor
+:
+    public tmp<Field<Type>>
+{
+public:
+
+    // Constructors
+
+        //- Construct from InputType
+        ConstPrecisionAdaptor(const Field<InputType>& input)
+        :
+            tmp<Field<Type>>()
+        {
+            if (std::is_same<Type, InputType>::value)
+            {
+                // Set reference - cast for compiler to handle different types
+                this->cref(reinterpret_cast<const Field<Type>&>(input));
+            }
+            else
+            {
+                this->reset(new Field<Type>(input.size()));
+                std::copy(input.cbegin(), input.cend(), this->ref().begin());
+            }
+        }
+
+
+    // Member Functions
+
+        //- Return the field
+        static const Field<Type>& get
+        (
+            const Field<InputType>& input,
+            Field<Type>& dst
+        )
+        {
+            if (std::is_same<Type, InputType>::value)
+            {
+                return reinterpret_cast<const Field<Type>&>(input);
+            }
+            else
+            {
+                dst.resize(input.size());
+                std::copy(input.cbegin(), input.cend(), dst.begin());
+                return dst;
+            }
+        }
+};
+
+
+//- A Field wrapper with possible data conversion
+template<class Type, class InputType>
+class PrecisionAdaptor
+:
+    public tmp<Field<Type>>
+{
+    // Private Data
+
+        //- Reference to underlying field
+        Field<InputType>& ref_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from Field<InputType>, copying on input as required
+        PrecisionAdaptor(Field<InputType>& input)
+        :
+            tmp<Field<Type>>(),
+            ref_(input)
+        {
+            if (std::is_same<Type, InputType>::value)
+            {
+                // Set reference - cast for compiler to handle different types
+                this->cref(reinterpret_cast<const Field<Type>&>(input));
+            }
+            else
+            {
+                this->reset(new Field<Type>(input.size()));
+                std::copy(input.cbegin(), input.cend(), this->ref().begin());
+            }
+        }
+
+
+    //- Destructor, copying on destroy
+    ~PrecisionAdaptor()
+    {
+        if (this->isTmp())
+        {
+            const Field<Type>& store = this->cref();
+            ref_.resize(store.size());
+            std::copy(store.cbegin(), store.cend(), ref_.begin());
+        }
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/primitiveFieldsFwd.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
@@ -48,6 +48,7 @@ template<class Type> class Field;
 
 typedef Field<label> labelField;
 typedef Field<scalar> scalarField;
+typedef Field<solveScalar> solveScalarField;
 typedef Field<vector> vectorField;
 typedef Field<sphericalTensor> sphericalTensorField;
 typedef Field<symmTensor> symmTensorField;

src/OpenFOAM/global/global.Cver

@@ -98,6 +98,7 @@ const std::string Foam::foamVersion::buildArch
 #endif
     ";label="  + std::to_string(8*sizeof(Foam::label))
   + ";scalar=" + std::to_string(8*sizeof(Foam::scalar))
+  + ";solve=" + std::to_string(8*sizeof(Foam::solveScalar))
 );
 
 

src/OpenFOAM/primitives/Scalar/Scalar.C

@@ -39,6 +39,7 @@ const Scalar pTraits<Scalar>::min = -ScalarVGREAT;
 const Scalar pTraits<Scalar>::max = ScalarVGREAT;
 const Scalar pTraits<Scalar>::rootMin = -ScalarROOTVGREAT;
 const Scalar pTraits<Scalar>::rootMax = ScalarROOTVGREAT;
+const Scalar pTraits<Scalar>::vsmall = ScalarVSMALL;
 
 const char* const pTraits<Scalar>::componentNames[] = { "" };
 

src/OpenFOAM/primitives/Scalar/Scalar.H

@@ -78,6 +78,7 @@ public:
         static const Scalar min;
         static const Scalar rootMax;
         static const Scalar rootMin;
+        static const Scalar vsmall;
 
 
     // Constructors

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

@@ -28,7 +28,7 @@ Typedef
 
 Description
     A floating-point number identical to float or double depending on
-    whether WM_SP or WM_DP is defined.
+    whether WM_SP, WM_SPDP or WM_DP is defined.
 
 SourceFiles
     scalar.C
@@ -43,13 +43,18 @@ SourceFiles
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#if defined(WM_SP)
+#if defined(WM_SP) || defined(WM_SPDP)
 
 // Define scalar as a float
 
 namespace Foam
 {
     typedef floatScalar scalar;
+    #if defined(WM_SPDP)
+    typedef doubleScalar solveScalar;
+    #else
+    typedef floatScalar solveScalar;
+    #endif
 
     constexpr scalar GREAT = floatScalarGREAT;
     constexpr scalar VGREAT = floatScalarVGREAT;
@@ -89,6 +94,7 @@ namespace Foam
 namespace Foam
 {
     typedef doubleScalar scalar;
+    typedef doubleScalar solveScalar;
 
     constexpr scalar GREAT = doubleScalarGREAT;
     constexpr scalar VGREAT = doubleScalarVGREAT;

src/OpenFOAM/primitives/Tensor/floatTensor/floatTensor.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) 2016 OpenFOAM Foundation
@@ -50,7 +50,7 @@ namespace Foam
 typedef Tensor<float> floatTensor;
 
 //- Data associated with floatTensor type are contiguous
-#if !defined(WM_SP)
+#if !defined(WM_SP) && !defined(WM_SPDP)
 template<>
 inline bool contiguous<floatTensor>() {return true;}
 #endif

src/OpenFOAM/primitives/Vector/floatVector/floatVector.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
@@ -50,7 +50,7 @@ namespace Foam
 typedef Vector<float> floatVector;
 
 //- Data associated with floatVector type are contiguous
-#if !defined(WM_SP)
+#if !defined(WM_SP) && !defined(WM_SPDP)
 template<>
 inline bool contiguous<floatVector>() {return true;}
 #endif

src/Pstream/mpi/UPstream.C

@@ -39,7 +39,7 @@ License
 #include <cstdlib>
 #include <csignal>
 
-#if defined(WM_SP)
+#if defined(WM_SP) || defined(WM_SPDP)
     #define MPI_SCALAR MPI_FLOAT
 #elif defined(WM_DP)
     #define MPI_SCALAR MPI_DOUBLE

src/dummyThirdParty/metisDecomp/metis.h

@@ -27,14 +27,14 @@
 
 
 // Float type: OpenFOAM uses WM_SP, WM_DP, metis.h uses REALTYPEWIDTH
-#if defined(WM_SP)
+#if defined(WM_SP) || defined(WM_SPDP)
   typedef float real_t;
   #define REALTYPEWIDTH 32
 #elif defined(WM_DP)
   typedef double real_t;
   #define REALTYPEWIDTH 64
 #else
-  #error "Incorrect user-supplied value for WM_SP / WM_DP  (metis REALTYPEWIDTH)"
+  #error "Incorrect user-supplied value for WM_SP (WM_SPDP) / WM_DP  (metis REALTYPEWIDTH)"
 #endif