From 55f287cd83e5df1573c64e04e0a12fcf304d8f4a Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 16 Mar 2022 15:57:31 +0100
Subject: [PATCH] ENH: limitVelocity, report count and percent of cells/faces
 (#2414)

- percent of cells is taken relative to selection size.
- percent of faces is taken relative to the number of boundary faces
  that do not fix velocity themselves.

ENH: avoid correctBoundaryConditions() if values were not limited
---
 .../corrections/limitVelocity/limitVelocity.C | 47 ++++++++++---
 .../velocityDampingConstraint.C               | 17 +++--
 .../limitTemperature/limitTemperature.C       | 66 ++++++++++++-------
 .../corrections/limitVelocity/limitVelocity.C | 63 +++++++++++++++---
 4 files changed, 149 insertions(+), 44 deletions(-)

diff --git a/src/faOptions/corrections/limitVelocity/limitVelocity.C b/src/faOptions/corrections/limitVelocity/limitVelocity.C
index 8fbb5756c0..186205ce64 100644
--- a/src/faOptions/corrections/limitVelocity/limitVelocity.C
+++ b/src/faOptions/corrections/limitVelocity/limitVelocity.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,41 +83,70 @@ void Foam::fa::limitVelocity::correct(areaVectorField& U)
 {
     const scalar maxSqrU = sqr(max_);
 
+    // Count nTotFaces ourselves
+    // (maybe only applying on a subset)
+    label nFacesAbove(0);
+    const label nTotFaces(returnReduce(faces_.size(), sumOp<label>()));
+
     vectorField& Uif = U.primitiveFieldRef();
 
     for (const label facei : faces_)
     {
-        const scalar magSqrUi = magSqr(Uif[facei]);
+        auto& Uval = Uif[facei];
+
+        const scalar magSqrUi = magSqr(Uval);
 
         if (magSqrUi > maxSqrU)
         {
-            Uif[facei] *= sqrt(maxSqrU/max(magSqrUi, SMALL));
+            Uval *= sqrt(maxSqrU/max(magSqrUi, SMALL));
+            ++nFacesAbove;
         }
     }
 
     // Handle boundaries in the case of 'all'
+    label nEdgesAbove(0);
+
     if (!faceSetOption::useSubMesh())
     {
         for (faPatchVectorField& Up : U.boundaryFieldRef())
         {
             if (!Up.fixesValue())
             {
-                forAll(Up, facei)
+                for (auto& Uval : Up)
                 {
-                    const scalar magSqrUi = magSqr(Up[facei]);
+                    const scalar magSqrUi = magSqr(Uval);
 
                     if (magSqrUi > maxSqrU)
                     {
-                        Up[facei] *= sqrt(maxSqrU/max(magSqrUi, SMALL));
+                        Uval *= sqrt(maxSqrU/max(magSqrUi, SMALL));
+                        ++nEdgesAbove;
                     }
                 }
             }
         }
     }
 
-    // We've changed internal values so give boundary conditions opportunity
-    // to correct.
-    U.correctBoundaryConditions();
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
+
+
+    reduce(nFacesAbove, sumOp<label>());
+    reduce(nEdgesAbove, sumOp<label>());
+
+    Info<< type() << ' ' << name_ << " Limited "
+        << nFacesAbove << " ("
+        << percent(nFacesAbove, nTotFaces)
+        << "%) of faces, with max limit " << max_ << endl;
+
+    if (nFacesAbove || nEdgesAbove)
+    {
+        // We've changed internal values so give
+        // boundary conditions opportunity to correct
+        U.correctBoundaryConditions();
+    }
 }
 
 
diff --git a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
index 372b5dedde..3794742ecc 100644
--- a/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
+++ b/src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C
@@ -62,7 +62,10 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
     const volVectorField& U = eqn.psi();
     scalarField& diag = eqn.diag();
 
-    label nDamped = 0;
+    // Count nTotCells ourselves
+    // (maybe only applying on a subset)
+    label nDamped(0);
+    const label nTotCells(returnReduce(cells_.size(), sumOp<label>()));
 
     for (label celli : cells_)
     {
@@ -79,10 +82,16 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
 
     reduce(nDamped, sumOp<label>());
 
-    Info<< type() << " " << name_ << " damped "
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
+
+    Info<< type() << ' ' << name_ << " damped "
         << nDamped << " ("
-        << 100*scalar(nDamped)/mesh_.globalData().nTotalCells()
-        << "%) of cells" << endl;
+        << percent(nDamped, nTotCells)
+        << "%) of cells, with max limit " << UMax_ << endl;
 }
 
 
diff --git a/src/fvOptions/corrections/limitTemperature/limitTemperature.C b/src/fvOptions/corrections/limitTemperature/limitTemperature.C
index a96d5d697d..b6dc4c333e 100644
--- a/src/fvOptions/corrections/limitTemperature/limitTemperature.C
+++ b/src/fvOptions/corrections/limitTemperature/limitTemperature.C
@@ -126,43 +126,53 @@ void Foam::fv::limitTemperature::correct(volScalarField& he)
     scalar Tmin0 = min(T);
     scalar Tmax0 = max(T);
 
-    label nOverTmax = 0;
-    label nLowerTmin = 0;
+    // Count nTotCells ourselves
+    // (maybe only applying on a subset)
+    label nBelowMin(0);
+    label nAboveMax(0);
+    const label nTotCells(returnReduce(cells_.size(), sumOp<label>()));
 
     forAll(cells_, i)
     {
         const label celli = cells_[i];
         if (hec[celli] < heMin[i])
         {
-            nLowerTmin++;
+            hec[celli] = heMin[i];
+            ++nBelowMin;
         }
         else if (hec[celli] > heMax[i])
         {
-            nOverTmax++;
+            hec[celli] = heMax[i];
+            ++nAboveMax;
         }
-        hec[celli]= max(min(hec[celli], heMax[i]), heMin[i]);
     }
 
-    reduce(nOverTmax, sumOp<label>());
-    reduce(nLowerTmin, sumOp<label>());
+    reduce(nBelowMin, sumOp<label>());
+    reduce(nAboveMax, sumOp<label>());
 
     reduce(Tmin0, minOp<scalar>());
     reduce(Tmax0, maxOp<scalar>());
 
-    Info<< type() << " " << name_ << " Lower limited "
-        << nLowerTmin << " ("
-        << 100*scalar(nLowerTmin)/mesh_.globalData().nTotalCells()
-        << "%) of cells" << endl;
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
 
-    Info<< type() << " " << name_ << " Upper limited "
-        << nOverTmax << " ("
-        << 100*scalar(nOverTmax)/mesh_.globalData().nTotalCells()
-        << "%) of cells" << endl;
+    Info<< type() << ' ' << name_ << " Lower limited " << nBelowMin << " ("
+        << percent(nBelowMin, nTotCells)
+        << "%) of cells, with min limit " << Tmin_ << endl;
+
+    Info<< type() << ' ' << name_ << " Upper limited " << nAboveMax << " ("
+        << percent(nAboveMax, nTotCells)
+        << "%) of cells, with max limit " << Tmax_ << endl;
+
+    Info<< type() << ' ' << name_ << " Unlimited Tmin " << Tmin0 << endl;
+    Info<< type() << ' ' << name_ << " Unlimited Tmax " << Tmax0 << endl;
 
-    Info<< type() << " " << name_ << " Unlimited Tmax " << Tmax0 << nl
-        <<  "Unlimited Tmin " << Tmin0 << endl;
 
     // Handle boundaries in the case of 'all'
+    bool changedValues = (nBelowMin || nAboveMax);
     if (!cellSetOption::useSubMesh())
     {
         volScalarField::Boundary& bf = he.boundaryFieldRef();
@@ -183,16 +193,28 @@ void Foam::fv::limitTemperature::correct(volScalarField& he)
 
                 forAll(hep, facei)
                 {
-                    hep[facei] =
-                        max(min(hep[facei], heMaxp[facei]), heMinp[facei]);
+                    if (hep[facei] < heMinp[facei])
+                    {
+                        hep[facei] = heMinp[facei];
+                        changedValues = true;
+                    }
+                    else if (hep[facei] > heMaxp[facei])
+                    {
+                        hep[facei] = heMaxp[facei];
+                        changedValues = true;
+                    }
                 }
             }
         }
     }
 
-    // We've changed internal values so give
-    // boundary conditions opportunity to correct
-    he.correctBoundaryConditions();
+
+    if (returnReduce(changedValues, orOp<bool>()))
+    {
+        // We've changed internal values so give
+        // boundary conditions opportunity to correct
+        he.correctBoundaryConditions();
+    }
 }
 
 
diff --git a/src/fvOptions/corrections/limitVelocity/limitVelocity.C b/src/fvOptions/corrections/limitVelocity/limitVelocity.C
index 2f8d1003ce..7f0d4a6bf2 100644
--- a/src/fvOptions/corrections/limitVelocity/limitVelocity.C
+++ b/src/fvOptions/corrections/limitVelocity/limitVelocity.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,41 +80,86 @@ void Foam::fv::limitVelocity::correct(volVectorField& U)
 {
     const scalar maxSqrU = sqr(max_);
 
+    // Count nTotCells ourselves
+    // (maybe only applying on a subset)
+    label nCellsAbove(0);
+    const label nTotCells(returnReduce(cells_.size(), sumOp<label>()));
+
     vectorField& Uif = U.primitiveFieldRef();
 
     for (const label celli : cells_)
     {
-        const scalar magSqrUi = magSqr(Uif[celli]);
+        auto& Uval = Uif[celli];
+
+        const scalar magSqrUi = magSqr(Uval);
 
         if (magSqrUi > maxSqrU)
         {
-            Uif[celli] *= sqrt(maxSqrU/magSqrUi);
+            Uval *= sqrt(maxSqrU/magSqrUi);
+            ++nCellsAbove;
         }
     }
 
     // Handle boundaries in the case of 'all'
+
+    label nFacesAbove(0);
+    label nTotFaces(0);
+
     if (!cellSetOption::useSubMesh())
     {
         for (fvPatchVectorField& Up : U.boundaryFieldRef())
         {
             if (!Up.fixesValue())
             {
-                forAll(Up, facei)
+                // Do not count patches that fix velocity themselves
+                nTotFaces += Up.size();
+
+                for (auto& Uval : Up)
                 {
-                    const scalar magSqrUi = magSqr(Up[facei]);
+                    const scalar magSqrUi = magSqr(Uval);
 
                     if (magSqrUi > maxSqrU)
                     {
-                        Up[facei] *= sqrt(maxSqrU/magSqrUi);
+                        Uval *= sqrt(maxSqrU/magSqrUi);
+                        ++nFacesAbove;
                     }
                 }
             }
         }
     }
 
-    // We've changed internal values so give
-    // boundary conditions opportunity to correct
-    U.correctBoundaryConditions();
+    // Percent, max 2 decimal places
+    const auto percent = [](scalar num, label denom) -> scalar
+    {
+        return (denom ? 1e-2*round(1e4*num/denom) : 0);
+    };
+
+
+    reduce(nCellsAbove, sumOp<label>());
+
+    // Report total numbers and percent
+    Info<< type() << ' ' << name_ << " Limited ";
+
+    Info<< nCellsAbove << " ("
+        << percent(nCellsAbove, nTotCells)
+        << "%) of cells";
+
+    reduce(nTotFaces, sumOp<label>());
+    reduce(nFacesAbove, sumOp<label>());
+    if (nTotFaces)
+    {
+        Info<< ", " << nFacesAbove << " ("
+            << percent(nFacesAbove, nTotFaces)
+            << "%) of faces";
+    }
+    Info<< ", with max limit " << max_ << endl;
+
+    if (nCellsAbove || nFacesAbove)
+    {
+        // We've changed internal values so give
+        // boundary conditions opportunity to correct
+        U.correctBoundaryConditions();
+    }
 }