ENH: turbulentDFSEMInlet: improve Reynolds-stress realizability checks

ENH: turbulentDFSEMInlet: new realizability checking function for scalar fields

src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C

@@ -811,71 +811,76 @@ turbulentDFSEMInletFvPatchVectorField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
+void Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
 (
-    const symmTensorField& Rf
+    const symmTensorField& R
 )
 {
-    // Perform checks of the stress tensor based on Cholesky decomposition
+    constexpr label maxDiffs = 5;
-    // constraints
+    label nDiffs = 0;
 
-    forAll(Rf, facei)
+    // (S:Eq. 4a-4c)
+    forAll(R, i)
     {
-        const symmTensor& R = Rf[facei];
+        bool diff = false;
 
-        if (R.xx() <= 0)
+        if (maxDiffs < nDiffs)
         {
-            FatalErrorInFunction
+            Info<< "More than " << maxDiffs << " times"
-                << "Reynolds stress " << R << " at face " << facei
+                << " Reynolds-stress realizability checks failed."
-                << " does not obey the constraint: R_xx > 0"
+                << " Skipping further comparisons." << endl;
-                << exit(FatalError);
+            return;
         }
 
-        const scalar a_xx = sqrt(R.xx());
+        const symmTensor& r = R[i];
 
-        const scalar a_xy = R.xy()/a_xx;
+        if (r.xx() < 0)
-
-        const scalar a_yy_2 = R.yy() - sqr(a_xy);
-
-        if (a_yy_2 < 0)
         {
-            FatalErrorInFunction
+            WarningInFunction
-                << "Reynolds stress " << R << " at face " << facei
+                << "Reynolds stress " << r << " at index " << i
-                << " leads to an invalid Cholesky decomposition due to the "
+                << " does not obey the constraint: Rxx >= 0"
-                << "constraint R_yy - sqr(a_xy) >= 0"
-                << exit(FatalError);
-        }
-
-        const scalar a_yy = Foam::sqrt(a_yy_2);
-
-        const scalar a_xz = R.xz()/a_xx;
-
-        const scalar a_yz = (R.yz() - a_xy*a_xz)/a_yy;
-
-        const scalar a_zz_2 = R.zz() - sqr(a_xz) - sqr(a_yz);
-
-        if (a_zz_2 < 0)
-        {
-            FatalErrorInFunction
-                << "Reynolds stress " << R << " at face " << facei
-                << " leads to an invalid Cholesky decomposition due to the "
-                << "constraint R_zz - sqr(a_xz) - sqr(a_yz) >= 0"
-                << exit(FatalError);
-        }
-
-        const scalar a_zz = Foam::sqrt(a_zz_2);
-
-        if (debug)
-        {
-            Pout<< "R: " << R
-                << " a_xx:" << a_xx << " a_xy:" << a_xy << " a_xz:" << a_xy
-                <<                     " a_yy:" << a_yy << " a_yz:" << a_yz
-                <<                                         " a_zz:" << a_zz
                 << endl;
+            diff = true;
+        }
+
+        if ((r.xx()*r.yy() - sqr(r.xy())) < 0)
+        {
+            WarningInFunction
+                << "Reynolds stress " << r << " at index " << i
+                << " does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0"
+                << endl;
+            diff = true;
+        }
+
+        if (det(r) < 0)
+        {
+            WarningInFunction
+                << "Reynolds stress " << r << " at index " << i
+                << " does not obey the constraint: det(R) >= 0"
+                << endl;
+            diff = true;
+        }
+
+        if (diff)
+        {
+            ++nDiffs;
+        }
     }
 }
 
-    return true;
+
+void Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
+(
+    const scalarField& R
+)
+{
+    if (min(R) <= 0)
+    {
+        FatalErrorInFunction
+            << "Reynolds stresses contain at least one "
+            << "nonpositive element. min(R) = " << min(R)
+            << exit(FatalError);
+    }
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H

@@ -361,8 +361,16 @@ public:
 
     // Member Functions
 
-        //- Return true if input Reynold stresses are valid
+        //- Check if input Reynolds stresses are valid
-        static bool checkStresses(const symmTensorField& Rf);
+        //  Realizability conditions (tag:S):
+        //      Schumann, U. (1977).
+        //      Realizability of Reynolds‐stress turbulence models.
+        //      The Physics of Fluids, 20(5), 721-725.
+        //      DOI:10.1063/1.861942
+        static void checkStresses(const symmTensorField& R);
+
+        //- Check if input Reynolds stresses are valid
+        static void checkStresses(const scalarField& R);
 
 
         // Mapping

src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/turbulentDigitalFilterInletFvPatchField.C

@@ -28,6 +28,7 @@ License
 #include "turbulentDigitalFilterInletFvPatchField.H"
 #include "addToRunTimeSelectionTable.H"
 #include "faceAreaWeightAMI.H"
+#include "turbulentDFSEMInletFvPatchVectorField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -284,6 +285,11 @@ turbulentDigitalFilterInletFvPatchField
             << "supported by the digital filter method."
             << endl;
     }
+
+    const scalar t = this->db().time().timeOutputValue();
+    const Field<TypeR> R(Rptr_->value(t));
+
+    turbulentDFSEMInletFvPatchVectorField::checkStresses(R);
 }