ENH: potentialFoam: add region functionality. Fixes #1223.

Also implements combination of -region and -dry-run

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -132,6 +132,7 @@ int main(int argc, char *argv[])
         "Execute functionObjects"
     );
 
+    #include "addRegionOption.H"
     #include "addCheckCaseOptions.H"
     #include "setRootCaseLists.H"
     #include "createTime.H"

applications/solvers/basic/potentialFoam/potentialFoam.C.orig

@@ -0,0 +1,254 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    potentialFoam
+
+Group
+    grpBasicSolvers
+
+Description
+    Potential flow solver which solves for the velocity potential, to
+    calculate the flux-field, from which the velocity field is obtained by
+    reconstructing the flux.
+
+    \heading Solver details
+    The potential flow solution is typically employed to generate initial fields
+    for full Navier-Stokes codes.  The flow is evolved using the equation:
+
+    \f[
+        \laplacian \Phi = \div(\vec{U})
+    \f]
+
+    Where:
+    \vartable
+        \Phi      | Velocity potential [m2/s]
+        \vec{U}   | Velocity [m/s]
+    \endvartable
+
+    The corresponding pressure field could be calculated from the divergence
+    of the Euler equation:
+
+    \f[
+        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+
+    but this generates excessive pressure variation in regions of large
+    velocity gradient normal to the flow direction.  A better option is to
+    calculate the pressure field corresponding to velocity variation along the
+    stream-lines:
+
+    \f[
+        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+    where the flow direction tensor \f$\vec{F}\f$ is obtained from
+    \f[
+        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
+    \f]
+
+    \heading Required fields
+    \plaintable
+        U         | Velocity [m/s]
+    \endplaintable
+
+    \heading Optional fields
+    \plaintable
+        p         | Kinematic pressure [m2/s2]
+        Phi       | Velocity potential [m2/s]
+                  | Generated from p (if present) or U if not present
+    \endplaintable
+
+    \heading Options
+    \plaintable
+        -writep   | write the Euler pressure
+        -writePhi | Write the final velocity potential
+        -initialiseUBCs | Update the velocity boundaries before solving for Phi
+    \endplaintable
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pisoControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Potential flow solver which solves for the velocity potential"
+    );
+
+    argList::addOption
+    (
+        "pName",
+        "pName",
+        "Name of the pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "initialiseUBCs",
+        "Initialise U boundary conditions"
+    );
+
+    argList::addBoolOption
+    (
+        "writePhi",
+        "Write the final velocity potential field"
+    );
+
+    argList::addBoolOption
+    (
+        "writep",
+        "Calculate and write the Euler pressure field"
+    );
+
+    argList::addBoolOption
+    (
+        "withFunctionObjects",
+        "Execute functionObjects"
+    );
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pisoControl potentialFlow(mesh, "potentialFlow");
+
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< nl << "Calculating potential flow" << endl;
+
+    // Since solver contains no time loop it would never execute
+    // function objects so do it ourselves
+    runTime.functionObjects().start();
+
+    MRF.makeRelative(phi);
+    adjustPhi(phi, U, p);
+
+    // Non-orthogonal velocity potential corrector loop
+    while (potentialFlow.correctNonOrthogonal())
+    {
+        fvScalarMatrix PhiEqn
+        (
+            fvm::laplacian(dimensionedScalar("1", dimless, 1), Phi)
+         ==
+            fvc::div(phi)
+        );
+
+        PhiEqn.setReference(PhiRefCell, PhiRefValue);
+        PhiEqn.solve();
+
+        if (potentialFlow.finalNonOrthogonalIter())
+        {
+            phi -= PhiEqn.flux();
+        }
+    }
+
+    MRF.makeAbsolute(phi);
+
+    Info<< "Continuity error = "
+        << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
+        << endl;
+
+    U = fvc::reconstruct(phi);
+    U.correctBoundaryConditions();
+
+    Info<< "Interpolated velocity error = "
+        << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
+        << endl;
+
+    // Write U and phi
+    U.write();
+    phi.write();
+
+    // Optionally write Phi
+    if (args.found("writePhi"))
+    {
+        Phi.write();
+    }
+
+    // Calculate the pressure field from the Euler equation
+    if (args.found("writep"))
+    {
+        Info<< nl << "Calculating approximate pressure field" << endl;
+
+        label pRefCell = 0;
+        scalar pRefValue = 0.0;
+        setRefCell
+        (
+            p,
+            potentialFlow.dict(),
+            pRefCell,
+            pRefValue
+        );
+
+        // Calculate the flow-direction filter tensor
+        volScalarField magSqrU(magSqr(U));
+        volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
+
+        // Calculate the divergence of the flow-direction filtered div(U*U)
+        // Filtering with the flow-direction generates a more reasonable
+        // pressure distribution in regions of high velocity gradient in the
+        // direction of the flow
+        volScalarField divDivUU
+        (
+            fvc::div
+            (
+                F & fvc::div(phi, U),
+                "div(div(phi,U))"
+            )
+        );
+
+        // Solve a Poisson equation for the approximate pressure
+        while (potentialFlow.correctNonOrthogonal())
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::laplacian(p) + divDivUU
+            );
+
+            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.solve();
+        }
+
+        p.write();
+    }
+
+    runTime.functionObjects().end();
+
+    runTime.printExecutionTime(Info);
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/include/createMesh.H

@@ -1,4 +1,5 @@
 Foam::autoPtr<Foam::fvMesh> meshPtr(nullptr);
+Foam::word regionName = Foam::fvMesh::defaultRegion;
 
 if (args.found("dry-run") || args.found("dry-run-write"))
 {
@@ -6,10 +7,16 @@ if (args.found("dry-run") || args.found("dry-run-write"))
         << "Operating in 'dry-run' mode: case will run for 1 time step.  "
         << "All checks assumed OK on a clean exit" << Foam::endl;
 
+    // Allow region in combination with dry-run
+    args.readIfPresent("region", regionName);
+
     Foam::FieldBase::allowConstructFromLargerSize = true;
 
     // Create a simplified 1D mesh and attempt to re-create boundary conditions
-    meshPtr.reset(new Foam::simplifiedMeshes::columnFvMesh(runTime));
+    meshPtr.reset
+    (
+        new Foam::simplifiedMeshes::columnFvMesh(runTime, regionName)
+    );
 
     // Stop after 1 iteration of the simplified mesh
 
@@ -29,10 +36,19 @@ if (args.found("dry-run") || args.found("dry-run-write"))
     Foam::functionObject::outputPrefix = "postProcessing-dry-run";
 }
 else
+{
+    if (args.readIfPresent("region", regionName))
+    {
+        Foam::Info
+            << "Create mesh " << regionName << " for time = "
+            << runTime.timeName() << Foam::nl << Foam::endl;
+    }
+    else
     {
         Foam::Info
             << "Create mesh for time = "
             << runTime.timeName() << Foam::nl << Foam::endl;
+    }
 
     meshPtr.reset
     (
@@ -40,7 +56,7 @@ else
         (
             Foam::IOobject
             (
-                Foam::fvMesh::defaultRegion,
+                regionName,
                 runTime.timeName(),
                 runTime,
                 Foam::IOobject::MUST_READ

src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C

@@ -67,14 +67,14 @@ Foam::simplifiedMeshes::simplifiedDynamicFvMeshBase::New
         if (cstrIter.found())
         {
             Info<< "Selecting simplified mesh model " << modelType << endl;
-            return autoPtr<dynamicFvMesh>(cstrIter()(io.time()));
+            return autoPtr<dynamicFvMesh>(cstrIter()(io.time(), io.name()));
         }
     }
 
     Info<< "Selecting simplified mesh model " << staticFvMesh::typeName << endl;
     return autoPtr<dynamicFvMesh>
     (
-        new SimplifiedDynamicFvMesh<staticFvMesh>(io.time())
+        new SimplifiedDynamicFvMesh<staticFvMesh>(io.time(), io.name())
     );
 }
 

src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.H

@@ -67,9 +67,10 @@ public:
             dynamicFvMesh,
             time,
             (
-                const Time& runTime
+                const Time& runTime,
+                const word& regionName
             ),
-            (runTime)
+            (runTime, regionName)
         );
 
 
@@ -101,7 +102,7 @@ public:
     ClassNameNoDebug(DynamicMeshType::typeName_.c_str());
 
     //- Constructor
-    SimplifiedDynamicFvMesh(const Time& runTime);
+    SimplifiedDynamicFvMesh(const Time& runTime, const word& regionName);
 
     //- Update the mesh for both mesh motion and topology change
     virtual bool update()

src/dynamicFvMesh/simplifiedDynamicFvMesh/simplifiedDynamicFvMeshTemplates.C

@@ -29,16 +29,17 @@ template<class DynamicMeshType>
 Foam::simplifiedMeshes::SimplifiedDynamicFvMesh<DynamicMeshType>::
 SimplifiedDynamicFvMesh
 (
-    const Time& runTime
+    const Time& runTime,
+    const word& regionName
 )
 :
     simplifiedDynamicFvMeshBase(),
-    columnFvMeshInfo(runTime),
+    columnFvMeshInfo(runTime, regionName),
     DynamicMeshType
     (
         IOobject
         (
-            fvMesh::defaultRegion,
+            regionName,
             runTime.constant(),
             runTime,
             IOobject::NO_READ, // Do not read any existing mesh

src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.C

@@ -62,7 +62,7 @@ bool Foam::simplifiedMeshes::columnFvMeshInfo::setPatchEntries
     (
         "boundary",
         localInstance_,
-        polyMesh::meshSubDir,
+        regionPrefix_ + polyMesh::meshSubDir,
         runTime,
         IOobject::MUST_READ,
         IOobject::NO_WRITE,
@@ -97,7 +97,12 @@ bool Foam::simplifiedMeshes::columnFvMeshInfo::setPatchEntries
     else
     {
         // No boundary file - try reading from a field
-        IOobjectList objects(runTime, runTime.timeName());
+        IOobjectList objects
+        (
+            runTime,
+            runTime.timeName(),
+            (regionName_ == fvMesh::defaultRegion ? "" : regionName_)
+        );
 
         if (objects.empty())
         {
@@ -202,7 +207,7 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::initialise(const Time& runTime)
         (
             "points",
             localInstance_,
-            polyMesh::meshSubDir,
+            regionPrefix_ + polyMesh::meshSubDir,
             runTime,
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
@@ -371,12 +376,7 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::addLocalPatches
 
     if (debug)
     {
-        Pout<< "patches:" << nl << endl;
+        Pout<< "patches:" << nl << mesh.boundaryMesh() << endl;
-        forAll(patches, patchi)
-        {
-            Pout<< "patch: " << patches[patchi]->name() << nl
-                << *patches[patchi] << endl;
-        }
     }
 }
 
@@ -400,13 +400,24 @@ void Foam::simplifiedMeshes::columnFvMeshInfo::initialiseZones(fvMesh& mesh)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo(const Time& runTime)
+Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo
+(
+    const Time& runTime,
+    const word& regionName
+)
 :
+    regionName_(regionName),
+    regionPrefix_
+    (
+        regionName_ == fvMesh::defaultRegion
+      ? ""
+      : regionName_ + '/'
+    ),
     localInstance_
     (
         runTime.findInstance
         (
-            polyMesh::meshSubDir,
+            regionPrefix_ + polyMesh::meshSubDir,
             "boundary",
             IOobject::READ_IF_PRESENT
         )
@@ -428,14 +439,18 @@ Foam::simplifiedMeshes::columnFvMeshInfo::columnFvMeshInfo(const Time& runTime)
 }
 
 
-Foam::simplifiedMeshes::columnFvMesh::columnFvMesh(const Time& runTime)
+Foam::simplifiedMeshes::columnFvMesh::columnFvMesh
+(
+    const Time& runTime,
+    const word& regionName
+)
 :
-    columnFvMeshInfo(runTime),
+    columnFvMeshInfo(runTime, regionName),
     simplifiedFvMesh
     (
         IOobject
         (
-            fvMesh::defaultRegion,
+            regionName,
             runTime.constant(),
             runTime,
             IOobject::NO_READ, // Do not read any existing mesh

src/finiteVolume/fvMesh/simplifiedFvMesh/columnFvMesh/columnFvMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -73,6 +73,12 @@ protected:
 
     // Protected data
 
+        //- Region of existing mesh
+        const word regionName_;
+
+        //- Additional prefix for region. Empty if default region
+        const fileName regionPrefix_;
+
         //- Location of existing mesh (if present)
         const word localInstance_;
 
@@ -118,7 +124,11 @@ public:
     ClassName("columnFvMeshInfo");
 
     // Constructor
-    columnFvMeshInfo(const Time& runTime);
+    columnFvMeshInfo
+    (
+        const Time& runTime,
+        const word& regionName
+    );
 };
 
 
@@ -134,7 +144,11 @@ public:
     TypeName("columnFvMesh");
 
     //- Constructor
-    columnFvMesh(const Time& runTime);
+    columnFvMesh
+    (
+        const Time& runTime,
+        const word& regionName = fvMesh::defaultRegion
+    );
 };