Data (e.g. slopes) extraction code for 2D flat plate case

utilities/extractData_2D/Make/files

@@ -0,0 +1,3 @@
+extractData_2D.C
+
+EXE = $(FOAM_USER_APPBIN)/extractData_2D

utilities/extractData_2D/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -DFULLDEBUG -g -O0
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools
+

utilities/extractData_2D/createFields.H

@@ -0,0 +1,49 @@
+Info<< "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Read transportProperties
+Info << "Reading transportProperties\n" << endl;
+IOdictionary transportProperties
+(
+	IOobject
+	(
+		"transportProperties",
+		runTime.constant(),
+		mesh,
+		IOobject::MUST_READ_IF_MODIFIED,
+		IOobject::NO_WRITE
+	)
+);
+
+// Read the constant nu
+dimensionedScalar nu
+(
+	"nu",
+	dimViscosity,
+	transportProperties.lookup("nu")
+);

utilities/extractData_2D/extractData_2D.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 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
+    extractData_2D
+
+Description
+    Extract data in 2D simulation
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Extract data in 2D flat plate."
+    );
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"    
+
+    #include "createFields.H"
+    #include "findCellFaceLabels.H"
+    #include "saveDataInit.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+    for (int timei = 1; timei < timeDirs.size(); timei++)
+    {
+        runTime.setTime(timeDirs[timei], timei);
+        Info << "Time = " << runTime.timeName() << endl;
+
+        #include "createFields.H"
+        #include "saveDataFields.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

utilities/extractData_2D/findCellFaceLabels.H

@@ -0,0 +1,88 @@
+// Required geometrical information
+label surfaceID(0);
+forAll (mesh.boundary(), patchI)
+{
+    if (mesh.boundary()[patchI].name() == "bottomWall")
+    {
+        surfaceID = patchI;
+    }
+}
+
+// Coordinates of cell centers and faces
+const vectorField Cf(mesh.Cf().internalField());
+const vectorField Cc(mesh.C().internalField());
+
+// Steps for getting IDs in y-coordinates
+int steps = 175;
+
+const polyPatch& bouPatch = mesh.Cf().boundaryField()[surfaceID].patch().patch();
+
+// Get global IDs of faces in boundary patch
+labelList patchFaceIDs(bouPatch.size());
+forAll (bouPatch,faceI)
+{
+    patchFaceIDs[faceI] = bouPatch.start()+faceI;
+}
+
+// Get IDs of primary adjacent cells
+const labelList& adjacentCellIDs = bouPatch.faceCells();
+
+// Vectors for face IDs and cell IDs
+std::vector<label> faceIDs;
+std::vector<label> cellIDs;
+
+// Get IDs of first cell face
+forAll (patchFaceIDs, faceI)
+{
+    faceIDs.push_back(
+        mesh.cells()[adjacentCellIDs[faceI]].opposingFaceLabel
+        (
+            patchFaceIDs[faceI],mesh.faces()
+        )
+    );
+}
+
+// Get IDs of second cell center
+label globFace = -1;
+forAll (patchFaceIDs, faceI)
+{
+    globFace = faceIDs[faceI];
+
+    if (mesh.owner()[globFace] == adjacentCellIDs[faceI])
+    {
+        cellIDs.push_back(mesh.neighbour()[globFace]);
+    }
+    else
+    {
+        cellIDs.push_back(mesh.owner()[globFace]);
+    }
+}
+
+// Get IDs of cell centers and faces iteratively
+for (int i = 0; i < steps; i++)
+{
+    forAll (patchFaceIDs, faceI)
+    {
+        faceIDs.push_back(
+            mesh.cells()[cellIDs[i*patchFaceIDs.size() + faceI]].opposingFaceLabel
+            (
+                faceIDs[i*patchFaceIDs.size() + faceI],mesh.faces()
+            )
+        );
+    }
+
+    globFace = -1;
+    forAll (patchFaceIDs, faceI)
+    {
+        globFace = faceIDs[(i+1)*patchFaceIDs.size() + faceI];
+
+        if (mesh.owner()[globFace] == cellIDs[i*patchFaceIDs.size() + faceI])
+        {
+            cellIDs.push_back(mesh.neighbour()[globFace]);
+        }
+        else
+        {
+            cellIDs.push_back(mesh.owner()[globFace]);
+        }
+    }
+}

utilities/extractData_2D/saveDataFields.H

@@ -0,0 +1,55 @@
+// Interpolate face values from cell center velocities
+surfaceVectorField U_face = fvc::interpolate(U);
+surfaceVectorField U_sngrad = fvc::snGrad(U);
+
+// Vector for average velocity in x-direction
+std::vector<scalar> averageU;
+
+// Calculate average velocity
+for (int i = 0; i < steps; i++)
+{
+    forAll (patchFaceIDs, faceI) // For U, it is actually a cell center value. (cellI)
+    {
+        scalar sum = 0;
+
+        if (i == 0)
+        {
+            averageU.push_back(U[adjacentCellIDs[faceI]].x());
+        }
+
+        else
+        {
+            sum = U[adjacentCellIDs[faceI]].x()*Cf[faceIDs[faceI]].y();
+
+            for (int index = 0; index < i; index++)
+            {
+                sum += U[cellIDs[index*patchFaceIDs.size() + faceI]].x()
+                *(Cf[faceIDs[(index + 1)*patchFaceIDs.size() + faceI]].y() - 
+                Cf[faceIDs[index*patchFaceIDs.size() + faceI]].y());
+            }            
+            averageU.push_back(sum/Cf[faceIDs[(i + 1)*patchFaceIDs.size() + faceI]].y());
+        }
+    }
+}
+
+forAll (patchFaceIDs, faceI)
+{
+    outputFile_wall << Cc[adjacentCellIDs[faceI]].x() << ", ";
+    outputFile_wall << Cc[adjacentCellIDs[faceI]].y() << ", ";
+    outputFile_wall << U[adjacentCellIDs[faceI]].x() << ", ";
+    outputFile_wall << mag(U_sngrad.boundaryField()[surfaceID][faceI]) << ", ";
+    outputFile_wall << "\n";
+}
+
+for (int i = 0; i < steps; i++)
+{
+    forAll (patchFaceIDs, faceI)
+    {
+        outputFile_face << Cf[faceIDs[i*patchFaceIDs.size() + faceI]].x() << ", ";
+        outputFile_face << Cf[faceIDs[i*patchFaceIDs.size() + faceI]].y() << ", ";
+        outputFile_face << averageU[i*patchFaceIDs.size() + faceI] << ", ";
+        outputFile_face << mag(U_sngrad[faceIDs[i*patchFaceIDs.size() + faceI]]) << ", ";
+        outputFile_face << U_face[faceIDs[i*patchFaceIDs.size() + faceI]].x();
+        outputFile_face << "\n";
+    }
+}

utilities/extractData_2D/saveDataInit.H

@@ -0,0 +1,11 @@
+// Save files
+string filename_wall = "extractData_wall_2D.csv";
+OFstream outputFile_wall(filename_wall);
+outputFile_wall.precision(12);
+
+string filename_face = "extractData_face_2D.csv";
+OFstream outputFile_face(filename_face);
+outputFile_face.precision(12);
+
+outputFile_wall << "x_center,y_1stcell,Ux_1stcell,xslope_wall" << "\n";
+outputFile_face << "x_center,y_face,avgU,xslope_face,Ux_face" << "\n";