diff --git a/applications/test/vector/Test-vector.C b/applications/test/vector/Test-vector.C
index 799dfd83f1fcb63b0c51991f6fbeec5a1d9686a9..f5b29425ca5de240d4afbcf518b6d8d6ee67e02d 100644
--- a/applications/test/vector/Test-vector.C
+++ b/applications/test/vector/Test-vector.C
@@ -45,6 +45,7 @@ void printInfo(const vector& vec)
         << " sum:" << cmptSum(vec)
         << " prod:" << cmptProduct(vec)
         << " mag:" << cmptMag(vec)
+        << " magSqr:" << cmptMagSqr(vec)
         << nl << nl;
 }
 
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
index b6a498640780f04fa0421fdf17caf3edbf061c53..acae8b4f37857d5b0feb69fade174c6f1a501402 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
@@ -322,6 +322,30 @@ tmp<Field<Type>> cmptMag(const tmp<Field<Type>>& tf)
 }
 
 
+template<class Type>
+void cmptMagSqr(Field<Type>& res, const UList<Type>& f)
+{
+    TFOR_ALL_F_OP_FUNC_F(Type, res, =, cmptMagSqr, Type, f)
+}
+
+template<class Type>
+tmp<Field<Type>> cmptMagSqr(const UList<Type>& f)
+{
+    auto tres = tmp<Field<Type>>::New(f.size());
+    cmptMagSqr(tres.ref(), f);
+    return tres;
+}
+
+template<class Type>
+tmp<Field<Type>> cmptMagSqr(const tmp<Field<Type>>& tf)
+{
+    auto tres = New(tf);
+    cmptMagSqr(tres.ref(), tf());
+    tf.clear();
+    return tres;
+}
+
+
 #define TMP_UNARY_FUNCTION(ReturnType, Func)                                   \
                                                                                \
 template<class Type>                                                           \
diff --git a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
index e80d01df2bfd0967b4cc0718b53c1a9f491d04f0..aef3f53c89a1ae22eae20fea75c102e1e2801416 100644
--- a/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
+++ b/src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
@@ -167,6 +167,16 @@ template<class Type>
 tmp<Field<Type>> cmptMag(const tmp<Field<Type>>& tf);
 
 
+template<class Type>
+void cmptMagSqr(Field<Type>& res, const UList<Type>& f);
+
+template<class Type>
+tmp<Field<Type>> cmptMagSqr(const UList<Type>& f);
+
+template<class Type>
+tmp<Field<Type>> cmptMagSqr(const tmp<Field<Type>>& tf);
+
+
 #define TMP_UNARY_FUNCTION(ReturnType, Func)                                   \
                                                                                \
 /** \brief  Apply the \c Func() function on the tmp field */                   \
diff --git a/src/OpenFOAM/primitives/Scalar/Scalar.H b/src/OpenFOAM/primitives/Scalar/Scalar.H
index 7caea028a0697118af94c48dc45c087a753ac54b..2d96ab3054f4f41d52dad6a75debe2042a59eaae 100644
--- a/src/OpenFOAM/primitives/Scalar/Scalar.H
+++ b/src/OpenFOAM/primitives/Scalar/Scalar.H
@@ -382,6 +382,12 @@ inline Scalar cmptMag(const Scalar s)
 }
 
 
+inline Scalar cmptMagSqr(const Scalar s)
+{
+    return magSqr(s);
+}
+
+
 inline Scalar sqrtSumSqr(const Scalar a, const Scalar b)
 {
     const Scalar maga = mag(a);
diff --git a/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H b/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
index d81a43cdfc55fcb5f7453d6fa1327dec41ae16de..641aafaad42dae2c9fd810544fd1ccc254551691 100644
--- a/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
+++ b/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
@@ -635,6 +635,18 @@ inline Form cmptMag
 }
 
 
+template<class Form, class Cmpt, direction Ncmpts>
+inline Form cmptMagSqr
+(
+    const VectorSpace<Form, Cmpt, Ncmpts>& vs
+)
+{
+    Form v;
+    VectorSpaceOps<Ncmpts,0>::eqOp(v, vs, eqMagSqrOp<Cmpt>());
+    return v;
+}
+
+
 template<class Form, class Cmpt, direction Ncmpts>
 inline Form max
 (
diff --git a/src/OpenFOAM/primitives/ops/ops.H b/src/OpenFOAM/primitives/ops/ops.H
index 7161d90cd5b91b61b8804d86439bd2c12a289ed1..9a007e4f40dbf6998be3bbad0c0f83f889e429fc 100644
--- a/src/OpenFOAM/primitives/ops/ops.H
+++ b/src/OpenFOAM/primitives/ops/ops.H
@@ -74,6 +74,7 @@ EqOp(multiplyEq, x *= y)
 EqOp(divideEq, x /= y)
 EqOp(eqSqr, x = sqr(y))
 EqOp(eqMag, x = mag(y))
+EqOp(eqMagSqr, x = magSqr(y))
 EqOp(plusEqMagSqr, x += magSqr(y))
 EqOp(maxEq, x = max(x, y))
 EqOp(minEq, x = min(x, y))