diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C
index cc0240be0b..0edc4f1273 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C
@@ -614,9 +614,12 @@ Foam::sqr(const dimensioned<Type>& dt)
 }
 
 template<class Type>
-Foam::dimensioned<Foam::scalar> Foam::magSqr(const dimensioned<Type>& dt)
+Foam::dimensioned<typename Foam::typeOfMag<Type>::type>
+Foam::magSqr(const dimensioned<Type>& dt)
 {
-    return dimensioned<scalar>
+    typedef typename typeOfMag<Type>::type magType;
+
+    return dimensioned<magType>
     (
         "magSqr(" + dt.name() + ')',
         magSqr(dt.dimensions()),
@@ -625,9 +628,12 @@ Foam::dimensioned<Foam::scalar> Foam::magSqr(const dimensioned<Type>& dt)
 }
 
 template<class Type>
-Foam::dimensioned<Foam::scalar> Foam::mag(const dimensioned<Type>& dt)
+Foam::dimensioned<typename Foam::typeOfMag<Type>::type>
+Foam::mag(const dimensioned<Type>& dt)
 {
-    return dimensioned<scalar>
+    typedef typename typeOfMag<Type>::type magType;
+
+    return dimensioned<magType>
     (
         "mag(" + dt.name() + ')',
         dt.dimensions(),
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
index 2e768cbda5..b7bd77bccb 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
@@ -437,10 +437,12 @@ dimensioned<typename outerProduct<Type, Type>::type>
 sqr(const dimensioned<Type>&);
 
 template<class Type>
-dimensioned<scalar> magSqr(const dimensioned<Type>&);
+dimensioned<typename typeOfMag<Type>::type>
+magSqr(const dimensioned<Type>& dt);
 
 template<class Type>
-dimensioned<scalar> mag(const dimensioned<Type>&);
+dimensioned<typename typeOfMag<Type>::type>
+mag(const dimensioned<Type>& dt);
 
 template<class Type>
 dimensioned<Type> cmptMultiply
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
index 2573a39afb..ae637a7417 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
@@ -141,13 +141,13 @@ sqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
 
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> magSqr
-(
-    const DimensionedField<Type, GeoMesh>& df
-)
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+magSqr(const DimensionedField<Type, GeoMesh>& df)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres =
-        tmp<DimensionedField<scalar, GeoMesh>>::New
+        tmp<DimensionedField<magType, GeoMesh>>::New
         (
             IOobject
             (
@@ -165,15 +165,15 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
 }
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> magSqr
-(
-    const tmp<DimensionedField<Type, GeoMesh>>& tdf
-)
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+magSqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     const DimensionedField<Type, GeoMesh>& df = tdf();
 
     auto tres =
-        reuseTmpDimensionedField<scalar, Type, GeoMesh>::New
+        reuseTmpDimensionedField<magType, Type, GeoMesh>::New
         (
             tdf,
             "magSqr(" + df.name() + ')',
@@ -188,13 +188,13 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
 
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> mag
-(
-    const DimensionedField<Type, GeoMesh>& df
-)
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+mag(const DimensionedField<Type, GeoMesh>& df)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres =
-        tmp<DimensionedField<scalar, GeoMesh>>::New
+        tmp<DimensionedField<magType, GeoMesh>>::New
         (
             IOobject
             (
@@ -212,15 +212,15 @@ tmp<DimensionedField<scalar, GeoMesh>> mag
 }
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> mag
-(
-    const tmp<DimensionedField<Type, GeoMesh>>& tdf
-)
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+mag(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     const DimensionedField<Type, GeoMesh>& df = tdf();
 
     auto tres =
-        reuseTmpDimensionedField<scalar, Type, GeoMesh>::New
+        reuseTmpDimensionedField<magType, Type, GeoMesh>::New
         (
             tdf,
             "mag(" + df.name() + ')',
@@ -325,11 +325,10 @@ UNARY_REDUCTION_FUNCTION(Type, max, gMax)
 UNARY_REDUCTION_FUNCTION(Type, min, gMin)
 UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
 UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
-
 UNARY_REDUCTION_FUNCTION(MinMax<Type>, minMax, gMinMax)
 UNARY_REDUCTION_FUNCTION(scalarMinMax, minMaxMag, gMinMaxMag)
 
-UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
+UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
 
 #undef UNARY_REDUCTION_FUNCTION
 
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.H
index 198aa981eb..1b055d94bb 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.H
@@ -62,34 +62,29 @@ tmp<DimensionedField<typename outerProduct<Type, Type>::type, GeoMesh>>
 sqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> magSqr
-(
-    const DimensionedField<Type, GeoMesh>& df
-);
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+magSqr(const DimensionedField<Type, GeoMesh>& df);
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> magSqr
-(
-    const tmp<DimensionedField<Type, GeoMesh>>& tdf
-);
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+magSqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> mag
-(
-    const DimensionedField<Type, GeoMesh>& df
-);
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+mag(const DimensionedField<Type, GeoMesh>& df);
 
 template<class Type, class GeoMesh>
-tmp<DimensionedField<scalar, GeoMesh>> mag
-(
-    const tmp<DimensionedField<Type, GeoMesh>>& tdf
-);
+tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
+mag(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 
 template<class Type, class GeoMesh>
 tmp
 <
     DimensionedField
-        <typename DimensionedField<Type, GeoMesh>::cmptType, GeoMesh>
+    <
+        typename DimensionedField<Type, GeoMesh>::cmptType,
+        GeoMesh
+    >
 >
 cmptAv(const DimensionedField<Type, GeoMesh>& df);
 
@@ -97,7 +92,10 @@ template<class Type, class GeoMesh>
 tmp
 <
     DimensionedField
-        <typename DimensionedField<Type, GeoMesh>::cmptType, GeoMesh>
+    <
+        typename DimensionedField<Type, GeoMesh>::cmptType,
+        GeoMesh
+    >
 >
 cmptAv(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 
@@ -123,7 +121,7 @@ UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
 UNARY_REDUCTION_FUNCTION(MinMax<Type>, minMax, gMinMax)
 UNARY_REDUCTION_FUNCTION(scalarMinMax, minMaxMag, gMinMaxMag)
 
-UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
+UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
 
 #undef UNARY_REDUCTION_FUNCTION
 
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
index 159c223198..2a5abbb551 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
@@ -158,7 +158,11 @@ sqr(const tmp<FieldField<Field, Type>>& tf)
 
 
 template<template<class> class Field, class Type>
-void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
+void magSqr
+(
+    FieldField<Field, typename typeOfMag<Type>::type>& sf,
+    const FieldField<Field, Type>& f
+)
 {
     forAll(sf, i)
     {
@@ -167,11 +171,14 @@ void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
 }
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f)
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+magSqr(const FieldField<Field, Type>& f)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres
     (
-        FieldField<Field, scalar>::NewCalculatedType(f)
+        FieldField<Field, magType>::NewCalculatedType(f)
     );
 
     magSqr(tres.ref(), f);
@@ -179,11 +186,14 @@ tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f)
 }
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf)
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+magSqr(const tmp<FieldField<Field, Type>>& tf)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres
     (
-        reuseTmpFieldField<Field, scalar, Type>::New(tf)
+        reuseTmpFieldField<Field, magType, Type>::New(tf)
     );
 
     magSqr(tres.ref(), tf());
@@ -193,7 +203,11 @@ tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf)
 
 
 template<template<class> class Field, class Type>
-void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
+void mag
+(
+    FieldField<Field, typename typeOfMag<Type>::type>& sf,
+    const FieldField<Field, Type>& f
+)
 {
     forAll(sf, i)
     {
@@ -202,11 +216,14 @@ void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
 }
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f)
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+mag(const FieldField<Field, Type>& f)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres
     (
-        FieldField<Field, scalar>::NewCalculatedType(f)
+        FieldField<Field, magType>::NewCalculatedType(f)
     );
 
     mag(tres.ref(), f);
@@ -214,11 +231,14 @@ tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f)
 }
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> mag(const tmp<FieldField<Field, Type>>& tf)
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+mag(const tmp<FieldField<Field, Type>>& tf)
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres
     (
-        reuseTmpFieldField<Field, scalar, Type>::New(tf)
+        reuseTmpFieldField<Field, magType, Type>::New(tf)
     );
 
     mag(tres.ref(), tf());
@@ -480,19 +500,21 @@ Type sum(const FieldField<Field, Type>& f)
 TMP_UNARY_FUNCTION(Type, sum)
 
 template<template<class> class Field, class Type>
-scalar sumMag(const FieldField<Field, Type>& f)
+typename typeOfMag<Type>::type sumMag(const FieldField<Field, Type>& f)
 {
-    scalar SumMag = 0.0;
+    typedef typename typeOfMag<Type>::type magType;
+
+    magType result = Zero;
 
     forAll(f, i)
     {
-        SumMag += sumMag(f[i]);
+        result += sumMag(f[i]);
     }
 
-    return SumMag;
+    return result;
 }
 
-TMP_UNARY_FUNCTION(scalar, sumMag)
+TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
 
 template<template<class> class Field, class Type>
 Type average(const FieldField<Field, Type>& f)
@@ -555,25 +577,24 @@ TMP_UNARY_FUNCTION(scalarMinMax, minMaxMag)
 
 
 // With reduction on ReturnType
-#define G_UNARY_FUNCTION(returnType, gFunc, func, rFunc)                       \
+#define G_UNARY_FUNCTION(ReturnType, gFunc, func, rFunc)                       \
                                                                                \
 template<template<class> class Field, class Type>                              \
-returnType gFunc(const FieldField<Field, Type>& f)                             \
+ReturnType gFunc(const FieldField<Field, Type>& f)                             \
 {                                                                              \
-    returnType res = func(f);                                                  \
-    reduce(res, rFunc##Op<returnType>());                                      \
+    ReturnType res = func(f);                                                  \
+    reduce(res, rFunc##Op<ReturnType>());                                      \
     return res;                                                                \
 }                                                                              \
-TMP_UNARY_FUNCTION(returnType, gFunc)
+TMP_UNARY_FUNCTION(ReturnType, gFunc)
 
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
-
 G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
 G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
 
-G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
+G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
 
 #undef G_UNARY_FUNCTION
 
diff --git a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.H b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.H
index 18be310352..9b1dd6b459 100644
--- a/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.H
+++ b/src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.H
@@ -90,23 +90,35 @@ sqr(const tmp<FieldField<Field, Type>>& tf);
 
 
 template<template<class> class Field, class Type>
-void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f);
+void magSqr
+(
+    FieldField<Field, typename typeOfMag<Type>::type>& sf,
+    const FieldField<Field, Type>& f
+);
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f);
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+magSqr(const FieldField<Field, Type>& f);
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf);
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+magSqr(const tmp<FieldField<Field, Type>>& tf);
 
 
 template<template<class> class Field, class Type>
-void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f);
+void mag
+(
+    FieldField<Field, typename typeOfMag<Type>::type>& res,
+    const FieldField<Field, Type>& f
+);
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f);
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+mag(const FieldField<Field, Type>& f);
 
 template<template<class> class Field, class Type>
-tmp<FieldField<Field, scalar>> mag(const tmp<FieldField<Field, Type>>& tf);
+tmp<FieldField<Field, typename typeOfMag<Type>::type>>
+mag(const tmp<FieldField<Field, Type>>& tf);
 
 
 template<template<class> class Field, class Type>
@@ -214,9 +226,9 @@ TMP_UNARY_FUNCTION(Type, sum)
 
 
 template<template<class> class Field, class Type>
-scalar sumMag(const FieldField<Field, Type>& f);
+typename typeOfMag<Type>::type sumMag(const FieldField<Field, Type>& f);
 
-TMP_UNARY_FUNCTION(scalar, sumMag)
+TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
 
 
 template<template<class> class Field, class Type>
@@ -239,20 +251,19 @@ TMP_UNARY_FUNCTION(scalarMinMax, minMaxMag)
 
 
 // With reduction on ReturnType
-#define G_UNARY_FUNCTION(returnType, gFunc, func, rFunc)                       \
+#define G_UNARY_FUNCTION(ReturnType, gFunc, func, rFunc)                       \
                                                                                \
 template<template<class> class Field, class Type>                              \
-returnType gFunc(const FieldField<Field, Type>& f);                            \
-TMP_UNARY_FUNCTION(returnType, gFunc)
+ReturnType gFunc(const FieldField<Field, Type>& f);                            \
+TMP_UNARY_FUNCTION(ReturnType, gFunc)
 
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
-
 G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
 G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
 
-G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
+G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
 
 #undef G_UNARY_FUNCTION
 
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
index 9404946e07..b6a4986407 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
@@ -146,23 +146,35 @@ sqr(const tmp<Field<Type>>& tf)
 
 
 template<class Type>
-void magSqr(Field<scalar>& res, const UList<Type>& f)
+void magSqr
+(
+    Field<typename typeOfMag<Type>::type>& res,
+    const UList<Type>& f
+)
 {
-    TFOR_ALL_F_OP_FUNC_F(scalar, res, =, magSqr, Type, f)
+    typedef typename typeOfMag<Type>::type magType;
+
+    TFOR_ALL_F_OP_FUNC_F(magType, res, =, magSqr, Type, f)
 }
 
 template<class Type>
-tmp<Field<scalar>> magSqr(const UList<Type>& f)
+tmp<Field<typename typeOfMag<Type>::type>>
+magSqr(const UList<Type>& f)
 {
-    auto tres = tmp<Field<scalar>>::New(f.size());
+    typedef typename typeOfMag<Type>::type magType;
+
+    auto tres = tmp<Field<magType>>::New(f.size());
     magSqr(tres.ref(), f);
     return tres;
 }
 
 template<class Type>
-tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf)
+tmp<Field<typename typeOfMag<Type>::type>>
+magSqr(const tmp<Field<Type>>& tf)
 {
-    auto tres = reuseTmp<scalar, Type>::New(tf);
+    typedef typename typeOfMag<Type>::type magType;
+
+    auto tres = reuseTmp<magType, Type>::New(tf);
     magSqr(tres.ref(), tf());
     tf.clear();
     return tres;
@@ -170,23 +182,35 @@ tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf)
 
 
 template<class Type>
-void mag(Field<scalar>& res, const UList<Type>& f)
+void mag
+(
+    Field<typename typeOfMag<Type>::type>& res,
+    const UList<Type>& f
+)
 {
-    TFOR_ALL_F_OP_FUNC_F(scalar, res, =, mag, Type, f)
+    typedef typename typeOfMag<Type>::type magType;
+
+    TFOR_ALL_F_OP_FUNC_F(magType, res, =, mag, Type, f)
 }
 
 template<class Type>
-tmp<Field<scalar>> mag(const UList<Type>& f)
+tmp<Field<typename typeOfMag<Type>::type>>
+mag(const UList<Type>& f)
 {
-    auto tres = tmp<Field<scalar>>::New(f.size());
+    typedef typename typeOfMag<Type>::type magType;
+
+    auto tres = tmp<Field<magType>>::New(f.size());
     mag(tres.ref(), f);
     return tres;
 }
 
 template<class Type>
-tmp<Field<scalar>> mag(const tmp<Field<Type>>& tf)
+tmp<Field<typename typeOfMag<Type>::type>>
+mag(const tmp<Field<Type>>& tf)
 {
-    auto tres = reuseTmp<scalar, Type>::New(tf);
+    typedef typename typeOfMag<Type>::type magType;
+
+    auto tres = reuseTmp<magType, Type>::New(tf);
     mag(tres.ref(), tf());
     tf.clear();
     return tres;
@@ -341,14 +365,14 @@ TMP_UNARY_FUNCTION(Type, min)
 template<class Type>
 Type sum(const UList<Type>& f)
 {
+    Type Sum = Zero;
+
     if (f.size())
     {
-        Type Sum = Zero;
         TFOR_ALL_S_OP_F(Type, Sum, +=, Type, f)
-        return Sum;
     }
 
-    return Zero;
+    return Sum;
 }
 
 TMP_UNARY_FUNCTION(Type, sum)
@@ -413,15 +437,15 @@ Type minMagSqr(const UList<Type>& f)
 TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 template<class Type>
-typename pTraits<Type>::cmptType
+typename scalarProduct<Type, Type>::type
 sumProd(const UList<Type>& f1, const UList<Type>& f2)
 {
-    typedef typename pTraits<Type>::cmptType outType;
+    typedef typename scalarProduct<Type, Type>::type prodType;
 
-    outType result = Zero;
+    prodType result = Zero;
     if (f1.size() && (f1.size() == f2.size()))
     {
-        TFOR_ALL_S_OP_F_OP_F(outType, result, +=, Type, f1, &&, Type, f2)
+        TFOR_ALL_S_OP_F_OP_F(prodType, result, +=, Type, f1, &&, Type, f2)
     }
     return result;
 }
@@ -450,41 +474,54 @@ Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2)
 
 
 template<class Type>
-scalar sumSqr(const UList<Type>& f)
+typename outerProduct1<Type>::type
+sumSqr(const UList<Type>& f)
 {
-    scalar SumSqr = 0;
+    typedef typename outerProduct1<Type>::type prodType;
+    prodType result = Zero;
     if (f.size())
     {
-        TFOR_ALL_S_OP_FUNC_F(scalar, SumSqr, +=, sqr, Type, f)
+        TFOR_ALL_S_OP_FUNC_F(prodType, result, +=, sqr, Type, f)
     }
-    return SumSqr;
+    return result;
 }
 
-TMP_UNARY_FUNCTION(scalar, sumSqr)
-
 template<class Type>
-scalar sumMag(const UList<Type>& f)
+typename outerProduct1<Type>::type
+sumSqr(const tmp<Field<Type>>& tf)
 {
-    scalar SumMag = 0;
-    if (f.size())
-    {
-        TFOR_ALL_S_OP_FUNC_F(scalar, SumMag, +=, mag, Type, f)
-    }
-    return SumMag;
+    typedef typename outerProduct1<Type>::type prodType;
+    prodType result = sumSqr(tf());
+    tf.clear();
+    return result;
 }
 
-TMP_UNARY_FUNCTION(scalar, sumMag)
+
+template<class Type>
+typename typeOfMag<Type>::type
+sumMag(const UList<Type>& f)
+{
+    typedef typename typeOfMag<Type>::type magType;
+    magType result = Zero;
+    if (f.size())
+    {
+        TFOR_ALL_S_OP_FUNC_F(magType, result, +=, mag, Type, f)
+    }
+    return result;
+}
+
+TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
 
 
 template<class Type>
 Type sumCmptMag(const UList<Type>& f)
 {
-    Type SumMag = Zero;
+    Type result = Zero;
     if (f.size())
     {
-        TFOR_ALL_S_OP_FUNC_F(scalar, SumMag, +=, cmptMag, Type, f)
+        TFOR_ALL_S_OP_FUNC_F(Type, result, +=, cmptMag, Type, f)
     }
-    return SumMag;
+    return result;
 }
 
 TMP_UNARY_FUNCTION(Type, sumCmptMag)
@@ -530,24 +567,24 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
 G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
 G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
 
-G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
-G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
+G_UNARY_FUNCTION(typename outerProduct1<Type>::type, gSumSqr, sumSqr, sum)
+G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
 
 #undef G_UNARY_FUNCTION
 
 
 template<class Type>
-typename pTraits<Type>::cmptType gSumProd
+typename scalarProduct<Type, Type>::type gSumProd
 (
     const UList<Type>& f1,
     const UList<Type>& f2,
     const label comm
 )
 {
-    typedef typename pTraits<Type>::cmptType outType;
+    typedef typename scalarProduct<Type, Type>::type prodType;
 
-    outType result = sumProd(f1, f2);
-    reduce(result, sumOp<outType>(), Pstream::msgType(), comm);
+    prodType result = sumProd(f1, f2);
+    reduce(result, sumOp<prodType>(), Pstream::msgType(), comm);
     return result;
 }
 
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
index c64aa4798a..e80d01df2b 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
@@ -94,23 +94,35 @@ sqr(const tmp<Field<Type>>& tf);
 
 
 template<class Type>
-void magSqr(Field<scalar>& res, const UList<Type>& f);
+void magSqr
+(
+    Field<typename typeOfMag<Type>::type>& res,
+    const UList<Type>& f
+);
 
 template<class Type>
-tmp<Field<scalar>> magSqr(const UList<Type>& f);
+tmp<Field<typename typeOfMag<Type>::type>>
+magSqr(const UList<Type>& f);
 
 template<class Type>
-tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf);
+tmp<Field<typename typeOfMag<Type>::type>>
+magSqr(const tmp<Field<Type>>& tf);
 
 
 template<class Type>
-void mag(Field<scalar>& res, const UList<Type>& f);
+void mag
+(
+    Field<typename typeOfMag<Type>::type>& res,
+    const UList<Type>& f
+);
 
 template<class Type>
-tmp<Field<scalar>> mag(const UList<Type>& f);
+tmp<Field<typename typeOfMag<Type>::type>>
+mag(const UList<Type>& f);
 
 template<class Type>
-tmp<Field<scalar>> mag(const tmp<Field<Type>>& tf);
+tmp<Field<typename typeOfMag<Type>::type>>
+mag(const tmp<Field<Type>>& tf);
 
 
 template<class Type>
@@ -197,7 +209,8 @@ TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 
 template<class Type>
-typename pTraits<Type>::cmptType sumProd
+typename scalarProduct<Type, Type>::type
+sumProd
 (
     const UList<Type>& f1,
     const UList<Type>& f2
@@ -207,14 +220,15 @@ template<class Type>
 Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2);
 
 template<class Type>
-scalar sumSqr(const UList<Type>& f);
-
-TMP_UNARY_FUNCTION(scalar, sumSqr)
+typename outerProduct1<Type>::type sumSqr(const UList<Type>& f);
 
 template<class Type>
-scalar sumMag(const UList<Type>& f);
+typename outerProduct1<Type>::type sumSqr(const tmp<Field<Type>>& tf);
 
-TMP_UNARY_FUNCTION(scalar, sumMag)
+template<class Type>
+typename typeOfMag<Type>::type sumMag(const UList<Type>& f);
+
+TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
 
 template<class Type>
 Type sumCmptMag(const UList<Type>& f);
@@ -244,13 +258,14 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
 G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
 G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
 
-G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
-G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
+G_UNARY_FUNCTION(typename outerProduct1<Type>::type, gSumSqr, sumSqr, sum)
+G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
 
 #undef G_UNARY_FUNCTION
 
+
 template<class Type>
-typename pTraits<Type>::cmptType gSumProd
+typename scalarProduct<Type, Type>::type gSumProd
 (
     const UList<Type>& f1,
     const UList<Type>& f2,
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
index 6c2a4d3b0c..b750be483d 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
@@ -257,7 +257,7 @@ sqr(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf)
 template<class Type, template<class> class PatchField, class GeoMesh>
 void magSqr
 (
-    GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
@@ -268,13 +268,16 @@ void magSqr
 
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+magSqr
 (
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres =
-        tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
+        tmp<GeometricField<magType, PatchField, GeoMesh>>::New
         (
             IOobject
             (
@@ -294,30 +297,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+magSqr
 (
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 {
-    const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
-
-    auto tres =
-        tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
-        (
-            IOobject
-            (
-                "magSqr(" + gf.name() + ')',
-                gf.instance(),
-                gf.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            gf.mesh(),
-            sqr(gf.dimensions())
-        );
-
-    magSqr(tres.ref(), gf);
-
+    auto tres = magSqr(tgf.cref());
     tgf.clear();
 
     return tres;
@@ -327,7 +313,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
 template<class Type, template<class> class PatchField, class GeoMesh>
 void mag
 (
-    GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
@@ -338,13 +324,16 @@ void mag
 
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+mag
 (
     const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 {
+    typedef typename typeOfMag<Type>::type magType;
+
     auto tres =
-        tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
+        tmp<GeometricField<magType, PatchField, GeoMesh>>::New
         (
             IOobject
             (
@@ -364,30 +353,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
 }
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+mag
 (
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 {
-    const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
-
-    auto tres =
-        tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
-        (
-            IOobject
-            (
-                "mag(" + gf.name() + ')',
-                gf.instance(),
-                gf.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            gf.mesh(),
-            gf.dimensions()
-        );
-
-    mag(tres.ref(), gf);
-
+    auto tres = mag(tgf.cref());
     tgf.clear();
 
     return tres;
@@ -559,7 +531,7 @@ dimensioned<returnType> func                                                   \
 
 UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
 UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
-UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
+UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
 
 #undef UNARY_REDUCTION_FUNCTION
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
index 142ffbee1f..e617526c6d 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
@@ -142,18 +142,20 @@ sqr(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf);
 template<class Type, template<class> class PatchField, class GeoMesh>
 void magSqr
 (
-    GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
     const GeometricField<Type, PatchField, GeoMesh>& gf
 );
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+magSqr
 (
     const GeometricField<Type, PatchField, GeoMesh>& gf
 );
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+magSqr
 (
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 );
@@ -161,18 +163,20 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
 template<class Type, template<class> class PatchField, class GeoMesh>
 void mag
 (
-    GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
     const GeometricField<Type, PatchField, GeoMesh>& gf
 );
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+mag
 (
     const GeometricField<Type, PatchField, GeoMesh>& gf
 );
 
 template<class Type, template<class> class PatchField, class GeoMesh>
-tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
+tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
+mag
 (
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 );
@@ -252,7 +256,7 @@ dimensioned<returnType> func                                                   \
 
 UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
 UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
-UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
+UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
 
 #undef UNARY_REDUCTION_FUNCTION
 
diff --git a/src/OpenFOAM/primitives/VectorSpace/products.H b/src/OpenFOAM/primitives/VectorSpace/products.H
index 36a0ce0c3f..8e7de12508 100644
--- a/src/OpenFOAM/primitives/VectorSpace/products.H
+++ b/src/OpenFOAM/primitives/VectorSpace/products.H
@@ -103,6 +103,16 @@ public:
 };
 
 
+//- Outer-product of identical types
+template<class arg1>
+class outerProduct1
+:
+    public outerProduct<arg1, arg1>
+{
+public:
+};
+
+
 template<class arg1, class arg2>
 class crossProduct
 {