ENH: added a variant of betaMax that computes it

based on the product of the Reynolds and Darcy numbers

src/optimisation/adjointOptimisation/adjoint/Make/files

@@ -193,6 +193,7 @@ $(topoVars)/dynamicTopODesignVariables/dynamicTopODesignVariables.C
 $(topoVars)/betaMax/betaMax/betaMax.C
 $(topoVars)/betaMax/value/betaMaxValue.C
 $(topoVars)/betaMax/Darcy/betaMaxDarcy.C
+$(topoVars)/betaMax/ReynoldsDarcy/betaMaxReynoldsDarcy.C
 $(topoVars)/betaMax/stepRamp/betaMaxStepRamp.C
 $(topoVars)/topOZones/topOZones.C
 $(topoVars)/regularisation/fieldRegularisation.C

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/Darcy/betaMaxDarcy.C

@@ -27,8 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "betaMaxDarcy.H"
-#include "EdgeMap.H"
-#include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -40,174 +38,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::betaMaxDarcy::computeLength(const dictionary& dict) const
-{
-    scalar length = Zero;
-    // If length is not provided explicitly, loop over the provided patches
-    // and compute the hydraulic diamater
-    const dictionary& DarcyDict = dict.subDict(type() + "Coeffs");
-    if (!DarcyDict.readIfPresent("length", length))
-    {
-        const labelHashSet inletPatches =
-            mesh_.boundaryMesh().patchSet
-            (
-                DarcyDict.get<wordRes>("inletPatches")
-            );
-
-        // If 2D, use the inlet area divided by the depth in the empty direction
-        if (mesh_.nGeometricD() != label(3))
-        {
-            // Accumulate area
-            for (const label pI : inletPatches)
-            {
-                const fvPatch& patch = mesh_.boundary()[pI];
-                length += gSum(patch.magSf());
-            }
-
-            // Divide with the span in the empty direction
-            const Vector<label>& geometricD = mesh_.geometricD();
-            const boundBox& bounds = mesh_.bounds();
-            forAll(geometricD, iDir)
-            {
-                if (geometricD[iDir] == -1)
-                {
-                    length /= bounds.span()[iDir];
-                }
-            }
-        }
-        // If 3D, use the inlet hydraulic diameter
-        else
-        {
-            scalar perimeter = Zero;
-            scalar area = Zero;
-            for (const label pI : inletPatches)
-            {
-                const polyPatch& patch = mesh_.boundaryMesh()[pI];
-                // Accumulate perimeter
-                const edgeList& edges = patch.edges();
-                const vectorField& points = patch.localPoints();
-                const label nInternalEdges = patch.nInternalEdges();
-                const label nEdges = patch.nEdges();
-                // Processor edges should not be accounted for
-                boolList isProcessorEdge = markProcessorEdges(patch);
-                label nActiveEdges(0);
-                forAll(isProcessorEdge, beI)
-                {
-                    if (!isProcessorEdge[beI])
-                    {
-                        perimeter += edges[nInternalEdges + beI].mag(points);
-                        nActiveEdges++;
-                    }
-                }
-
-                if (debug > 1)
-                {
-                    Info<< "patch:: " << patch.name() << nl
-                        << "Number of boundary edges "
-                        << returnReduce(nEdges - nInternalEdges, sumOp<label>())
-                        << ", Number of non-processor edges "
-                        << returnReduce(nActiveEdges, sumOp<label>()) << endl;
-                }
-
-                // Accumulate area
-                area += sum(patch.magFaceAreas());
-            }
-
-            reduce(perimeter, sumOp<scalar>());
-            reduce(area, sumOp<scalar>());
-
-            length = scalar(4)*area/perimeter;
-        }
-    }
-
-    return length;
-}
-
-
-Foam::boolList Foam::betaMaxDarcy::markProcessorEdges
-(
-    const polyPatch& patch
-) const
-{
-    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-    const label nNonProcessor = pbm.nNonProcessor();
-
-    // Processor edges will artificially increase the perimeter of the inlet.
-    // We need to purge them.
-
-    // Mark all edges connected to a processor patch
-    label nProcEdges(0);
-    for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
-    {
-        const polyPatch& procPatch = pbm[procI];
-        nProcEdges += procPatch.nEdges() - procPatch.nInternalEdges();
-    }
-
-    EdgeMap<bool> isInletEdge(nProcEdges);
-
-    for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
-    {
-        const polyPatch& procPatch = pbm[procI];
-        const labelList& procMp = procPatch.meshPoints();
-        const label procInternalEdges = procPatch.nInternalEdges();
-        const label procEdges = procPatch.nEdges();
-
-        for (label edgeI = procInternalEdges; edgeI < procEdges; ++edgeI)
-        {
-            const edge& e = procPatch.edges()[edgeI];
-            const edge meshE = edge(procMp[e[0]], procMp[e[1]]);
-            isInletEdge.insert(meshE, false);
-        }
-    }
-
-    // Loop over inlet edges to mark common edges with processor patches
-    const label nInternalEdges = patch.nInternalEdges();
-    const label nEdges = patch.nEdges();
-
-    const edgeList& edges = patch.edges();
-    const labelList& mp = patch.meshPoints();
-    for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
-    {
-        const edge& e = edges[edgeI];
-        const edge meshE = edge(mp[e[0]], mp[e[1]]);
-        auto iter = isInletEdge.find(meshE);
-
-        if (iter.found())
-        {
-            iter.val() = true;
-        }
-    }
-
-    // A special case is a processor patch intersecting the inlet patches on a
-    // (true) boundary edge of the latter. Use syncEdgeMap to make sure these
-    // edges don't make it to the final list
-    syncTools::syncEdgeMap
-    (
-        pbm.mesh(),
-        isInletEdge,
-        andEqOp<bool>()
-    );
-
-    // Now loop over the inlet faces once more and mark the common edges with
-    // processor boundaries
-    boolList isProcessorEdge(nEdges - nInternalEdges, false);
-    for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
-    {
-        const edge& e = edges[edgeI];
-        const edge meshE = edge(mp[e[0]], mp[e[1]]);
-        const auto iter = isInletEdge.cfind(meshE);
-
-        if (iter.found() && iter.val())
-        {
-            isProcessorEdge[edgeI - nInternalEdges] = true;
-        }
-    }
-
-    return isProcessorEdge;
-}
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::betaMaxDarcy::betaMaxDarcy
@@ -238,7 +68,7 @@ Foam::betaMaxDarcy::betaMaxDarcy
         ).get<dimensionedScalar>("nu").value();
 
     value_ = nu/DarcyNumber_/length_/length_;
-    DebugInfo
+    Info
         << "Computed a betaMax value of " << value_
         << " based on a length of " << length_ << nl << endl;
 }

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/Darcy/betaMaxDarcy.H

@@ -95,15 +95,6 @@ protected:
         scalar length_;
 
 
-    // Protected Member Functions
-
-        //- Compute the characteristic length
-        scalar computeLength(const dictionary& dict) const;
-
-        //- Mark all common inlet - processor edges
-        boolList markProcessorEdges(const polyPatch& patch) const;
-
-
 public:
 
     //- Runtime type information

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/ReynoldsDarcy/betaMaxReynoldsDarcy.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2023 PCOpt/NTUA
+    Copyright (C) 2020-2023 FOSS GP
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "betaMaxReynoldsDarcy.H"
+#include "EdgeMap.H"
+#include "syncTools.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(betaMaxReynoldsDarcy, 0);
+    addToRunTimeSelectionTable(betaMax, betaMaxReynoldsDarcy, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::betaMaxReynoldsDarcy::betaMaxReynoldsDarcy
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    betaMax(mesh, dict),
+    ReynoldsDarcyNumber_
+    (
+        dict.subDict(type() + "Coeffs").getOrDefault<scalar>("ReDa", 1.e-5)
+    ),
+    length_(computeLength(dict)),
+    Uref_(dict.subDict(type() + "Coeffs").get<scalar>("Uref"))
+{
+    value_ = Uref_/ReynoldsDarcyNumber_/length_;
+    Info
+        << "Computed a betaMax value of " << value_
+        << " based on a length of " << length_ << nl << endl;
+}
+
+
+// ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/ReynoldsDarcy/betaMaxReynoldsDarcy.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020-2023 PCOpt/NTUA
+    Copyright (C) 2020-2023 FOSS GP
+-------------------------------------------------------------------------------
+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::betaMaxReynoldsDarcy
+
+Description
+    Compute betaMax through the products of the Reynolds and Darcy numbers,
+    quantifying the momentum-to-porous forces ratio
+        ReDa = Uref/betaMax/L
+    where Uref is a reference velocity and L is a characteristic length.
+    The latter is either supplied directly or computed as the
+    - 2D: area of the inlet patches divided by the span in the empty direction
+    - 3D: the hydraulic diameter of the inlet patches
+
+SourceFiles
+    betaMaxReynoldsDarcy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef betaMaxReynoldsDarcy_H
+#define betaMaxReynoldsDarcy_H
+
+#include "betaMax.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class betaMaxReynoldsDarcy Declaration
+\*---------------------------------------------------------------------------*/
+
+class betaMaxReynoldsDarcy
+:
+    public betaMax
+{
+
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        betaMaxReynoldsDarcy(const betaMaxReynoldsDarcy&) = delete;
+
+        //- No copy assignment
+        void operator=(const betaMaxReynoldsDarcy&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- The Darcy number expressing the ratio of viscous to porous forces
+        scalar ReynoldsDarcyNumber_;
+
+        //- Characteristic length of the case
+        //  Either supplied directly or computed as the hydraulic diameter of
+        //  the inlet
+        scalar length_;
+
+        //- Reference velocity
+        scalar Uref_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ReynoldsDarcy");
+
+
+    // Constructors
+
+        //- Construct from components
+        betaMaxReynoldsDarcy
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~betaMaxReynoldsDarcy() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.C

@@ -27,6 +27,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "betaMax.H"
+#include "EdgeMap.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -37,6 +39,176 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::betaMax::computeLength(const dictionary& dict) const
+{
+    scalar length = Zero;
+    // If length is not provided explicitly, loop over the provided patches
+    // and compute the hydraulic diamater
+    const dictionary& DarcyDict = dict.subDict(type() + "Coeffs");
+    if (!DarcyDict.readIfPresent("length", length))
+    {
+        const labelHashSet inletPatches =
+            mesh_.boundaryMesh().patchSet
+            (
+                DarcyDict.get<wordRes>("inletPatches")
+            );
+
+        // If 2D, use the inlet area divided by the depth in the empty direction
+        if (mesh_.nGeometricD() != label(3))
+        {
+            // Accumulate area
+            for (const label pI : inletPatches)
+            {
+                const fvPatch& patch = mesh_.boundary()[pI];
+                length += gSum(patch.magSf());
+            }
+
+            // Divide with the span in the empty direction
+            const Vector<label>& geometricD = mesh_.geometricD();
+            const boundBox& bounds = mesh_.bounds();
+            forAll(geometricD, iDir)
+            {
+                if (geometricD[iDir] == -1)
+                {
+                    length /= bounds.span()[iDir];
+                }
+            }
+        }
+        // If 3D, use the inlet hydraulic diameter
+        else
+        {
+            scalar perimeter = Zero;
+            scalar area = Zero;
+            for (const label pI : inletPatches)
+            {
+                const polyPatch& patch = mesh_.boundaryMesh()[pI];
+                // Accumulate perimeter
+                const edgeList& edges = patch.edges();
+                const vectorField& points = patch.localPoints();
+                const label nInternalEdges = patch.nInternalEdges();
+                const label nEdges = patch.nEdges();
+                // Processor edges should not be accounted for
+                boolList isProcessorEdge = markProcessorEdges(patch);
+                label nActiveEdges(0);
+                forAll(isProcessorEdge, beI)
+                {
+                    if (!isProcessorEdge[beI])
+                    {
+                        perimeter += edges[nInternalEdges + beI].mag(points);
+                        nActiveEdges++;
+                    }
+                }
+
+                if (debug > 1)
+                {
+                    Info<< "patch:: " << patch.name() << nl
+                        << "Number of boundary edges "
+                        << returnReduce(nEdges - nInternalEdges, sumOp<label>())
+                        << ", Number of non-processor edges "
+                        << returnReduce(nActiveEdges, sumOp<label>()) << endl;
+                }
+
+                // Accumulate area
+                area += sum(patch.magFaceAreas());
+            }
+
+            reduce(perimeter, sumOp<scalar>());
+            reduce(area, sumOp<scalar>());
+
+            length = scalar(4)*area/perimeter;
+        }
+    }
+
+    return length;
+}
+
+
+Foam::boolList Foam::betaMax::markProcessorEdges
+(
+    const polyPatch& patch
+) const
+{
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+    const label nNonProcessor = pbm.nNonProcessor();
+
+    // Processor edges will artificially increase the perimeter of the inlet.
+    // We need to purge them.
+
+    // Mark all edges connected to a processor patch
+    label nProcEdges(0);
+    for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
+    {
+        const polyPatch& procPatch = pbm[procI];
+        nProcEdges += procPatch.nEdges() - procPatch.nInternalEdges();
+    }
+
+    EdgeMap<bool> isInletEdge(nProcEdges);
+
+    for (label procI = nNonProcessor; procI < pbm.size() ; ++procI)
+    {
+        const polyPatch& procPatch = pbm[procI];
+        const labelList& procMp = procPatch.meshPoints();
+        const label procInternalEdges = procPatch.nInternalEdges();
+        const label procEdges = procPatch.nEdges();
+
+        for (label edgeI = procInternalEdges; edgeI < procEdges; ++edgeI)
+        {
+            const edge& e = procPatch.edges()[edgeI];
+            const edge meshE = edge(procMp[e[0]], procMp[e[1]]);
+            isInletEdge.insert(meshE, false);
+        }
+    }
+
+    // Loop over inlet edges to mark common edges with processor patches
+    const label nInternalEdges = patch.nInternalEdges();
+    const label nEdges = patch.nEdges();
+
+    const edgeList& edges = patch.edges();
+    const labelList& mp = patch.meshPoints();
+    for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
+    {
+        const edge& e = edges[edgeI];
+        const edge meshE = edge(mp[e[0]], mp[e[1]]);
+        auto iter = isInletEdge.find(meshE);
+
+        if (iter.found())
+        {
+            iter.val() = true;
+        }
+    }
+
+    // A special case is a processor patch intersecting the inlet patches on a
+    // (true) boundary edge of the latter. Use syncEdgeMap to make sure these
+    // edges don't make it to the final list
+    syncTools::syncEdgeMap
+    (
+        pbm.mesh(),
+        isInletEdge,
+        andEqOp<bool>()
+    );
+
+    // Now loop over the inlet faces once more and mark the common edges with
+    // processor boundaries
+    boolList isProcessorEdge(nEdges - nInternalEdges, false);
+    for (label edgeI = nInternalEdges; edgeI < nEdges; ++edgeI)
+    {
+        const edge& e = edges[edgeI];
+        const edge meshE = edge(mp[e[0]], mp[e[1]]);
+        const auto iter = isInletEdge.cfind(meshE);
+
+        if (iter.found() && iter.val())
+        {
+            isProcessorEdge[edgeI - nInternalEdges] = true;
+        }
+    }
+
+    return isProcessorEdge;
+}
+
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::betaMax::betaMax

src/optimisation/adjointOptimisation/adjoint/optimisation/designVariables/topODesignVariables/betaMax/betaMax/betaMax.H

@@ -75,6 +75,15 @@ protected:
         scalar value_;
 
 
+    // Protected Member Functions
+
+        //- Compute the characteristic length
+        scalar computeLength(const dictionary& dict) const;
+
+        //- Mark all common inlet - processor edges
+        boolList markProcessorEdges(const polyPatch& patch) const;
+
+
 public:
 
     //- Runtime type information