From 3ca3f962dc3590080276dd4f7f70b29cf0423044 Mon Sep 17 00:00:00 2001 From: Bernardo Pacini Date: Mon, 1 Jul 2024 21:21:14 -0400 Subject: [PATCH 1/2] Adding new post-processing routines to compute surface forces --- .../calcForceCompressible/Make/files | 3 + .../calcForceCompressible/Make/options | 21 +++ .../calcForceCompressible/calcForce.C | 147 ++++++++++++++++ .../calcForceCompressible/createFields.H | 55 ++++++ .../calcForceIncompressible/Make/files | 3 + .../calcForceIncompressible/Make/options | 20 +++ .../calcForceIncompressible/calcForce.C | 162 ++++++++++++++++++ .../calcForceIncompressible/createFields.H | 48 ++++++ 8 files changed, 459 insertions(+) create mode 100755 src/utilities/postProcessing/calcForceCompressible/Make/files create mode 100755 src/utilities/postProcessing/calcForceCompressible/Make/options create mode 100644 src/utilities/postProcessing/calcForceCompressible/calcForce.C create mode 100755 src/utilities/postProcessing/calcForceCompressible/createFields.H create mode 100755 src/utilities/postProcessing/calcForceIncompressible/Make/files create mode 100755 src/utilities/postProcessing/calcForceIncompressible/Make/options create mode 100755 src/utilities/postProcessing/calcForceIncompressible/calcForce.C create mode 100755 src/utilities/postProcessing/calcForceIncompressible/createFields.H diff --git a/src/utilities/postProcessing/calcForceCompressible/Make/files b/src/utilities/postProcessing/calcForceCompressible/Make/files new file mode 100755 index 00000000..4fb73b3b --- /dev/null +++ b/src/utilities/postProcessing/calcForceCompressible/Make/files @@ -0,0 +1,3 @@ +calcForce.C + +EXE = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedBins/calcForceCompressible diff --git a/src/utilities/postProcessing/calcForceCompressible/Make/options b/src/utilities/postProcessing/calcForceCompressible/Make/options new file mode 100755 index 00000000..e837721a --- /dev/null +++ b/src/utilities/postProcessing/calcForceCompressible/Make/options @@ -0,0 +1,21 @@ +EXE_INC = \ + -std=c++11 \ + -Wno-old-style-cast \ + -Wno-conversion-null \ + -Wno-deprecated-copy \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lfiniteVolume \ + -lsampling \ + -lmeshTools \ + -lfvOptions + diff --git a/src/utilities/postProcessing/calcForceCompressible/calcForce.C b/src/utilities/postProcessing/calcForceCompressible/calcForce.C new file mode 100644 index 00000000..d6a17f10 --- /dev/null +++ b/src/utilities/postProcessing/calcForceCompressible/calcForce.C @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + + DAFoam : Discrete Adjoint with OpenFOAM + Version : v3 + + Description: + Calculating force per surface area. This will be used to plot the + spanwise force distribution + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "argList.H" +#include "autoPtr.H" +#include "Time.H" +#include "timeSelector.H" +#include "TimePaths.H" +#include "ListOps.H" +#include "fvMesh.H" +#include "OFstream.H" +#include "simpleControl.H" +#include "fvOptions.H" +#include "fluidThermo.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char* argv[]) +{ + Info << "Computing forces...." << endl; + + argList::addOption( + "patchNames", + "'(wing)'", + "List of patch names to compute"); + + argList::addOption( + "time", + "1000", + "Time instance to compute, if not provided runs all times"); + +#include "setRootCase.H" +#include "createTime.H" + + if (!args.optionFound("time")) + { + Info << "Time not set, running all times." << endl; + } + + // Create the processor databases + autoPtr timePaths; + timePaths = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); + + List patchNames; + if (args.optionFound("patchNames")) + { + patchNames = wordReList(args.optionLookup("patchNames")()); + } + else + { + Info << "patchNames not set! Exit." << endl; + return 1; + } + + forAll(timeDirs, iTime) + { + runTime.setTime(timeDirs[iTime].value(), 0); + +#include "createMesh.H" +#include "createFields.H" + + // Initialize surface field for face-centered forces + volVectorField volumeForceField( + IOobject( + "force", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + mesh, + dimensionedVector("surfaceForce", dimensionSet(1, 1, -2, 0, 0, 0, 0), vector::zero), + fixedValueFvPatchScalarField::typeName); + + // this code is pulled from: + // src/functionObjects/forcces/forces.C + // modified slightly + vector forces(vector::zero); + + const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); + const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); + + volSymmTensorField devRhoReff( + IOobject( + IOobject::groupName("devRhoReff", U.group()), + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + (-rho * nuEff) * dev(twoSymm(fvc::grad(U)))); + + const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); + + vector totalForce(vector::zero); + forAll(patchNames, cI) + { + // get the patch id label + label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); + if (patchI < 0) + { + Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; + return 1; + } + // create a shorter handle for the boundary patch + const fvPatch& patch = mesh.boundary()[patchI]; + // normal force + vectorField fN(Sfb[patchI] * p.boundaryField()[patchI] * rho); + // tangential force + vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); + // sum them up + forAll(patch, faceI) + { + // compute forces + forces.x() = fN[faceI].x() + fT[faceI].x(); + forces.y() = fN[faceI].y() + fT[faceI].y(); + forces.z() = fN[faceI].z() + fT[faceI].z(); + + // project force direction + volumeForceField.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; + + totalForce.x() += forces.x(); + totalForce.y() += forces.y(); + totalForce.z() += forces.z(); + } + } + volumeForceField.write(); + + Info << "Total force: " << totalForce << endl; + + Info << "Computing forces.... Completed!" << endl; + } + return 0; +} + +// ************************************************************************* // diff --git a/src/utilities/postProcessing/calcForceCompressible/createFields.H b/src/utilities/postProcessing/calcForceCompressible/createFields.H new file mode 100755 index 00000000..0ec75421 --- /dev/null +++ b/src/utilities/postProcessing/calcForceCompressible/createFields.H @@ -0,0 +1,55 @@ +autoPtr pThermo +( + fluidThermo::New(mesh) +); +fluidThermo& thermo = pThermo(); +thermo.validate(args.executable(), "h", "e"); + +volScalarField& p = thermo.p(); + +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info<< "Reading field nut\n" << endl; +volScalarField nut +( + IOobject + ( + "nut", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +volScalarField nuEff = nut + thermo.nu(); + +#include "compressibleCreatePhi.H" + + diff --git a/src/utilities/postProcessing/calcForceIncompressible/Make/files b/src/utilities/postProcessing/calcForceIncompressible/Make/files new file mode 100755 index 00000000..5ed44cb8 --- /dev/null +++ b/src/utilities/postProcessing/calcForceIncompressible/Make/files @@ -0,0 +1,3 @@ +calcForce.C + +EXE = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedBins/calcForceIncompressible diff --git a/src/utilities/postProcessing/calcForceIncompressible/Make/options b/src/utilities/postProcessing/calcForceIncompressible/Make/options new file mode 100755 index 00000000..86de1a4a --- /dev/null +++ b/src/utilities/postProcessing/calcForceIncompressible/Make/options @@ -0,0 +1,20 @@ +EXE_INC = \ + -std=c++11 \ + -Wno-old-style-cast \ + -Wno-conversion-null \ + -Wno-deprecated-copy \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + + +EXE_LIBS = \ + -lturbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools \ + -lfvOptions \ + -lsampling + diff --git a/src/utilities/postProcessing/calcForceIncompressible/calcForce.C b/src/utilities/postProcessing/calcForceIncompressible/calcForce.C new file mode 100755 index 00000000..7b58b940 --- /dev/null +++ b/src/utilities/postProcessing/calcForceIncompressible/calcForce.C @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + + DAFoam : Discrete Adjoint with OpenFOAM + Version : v3 + + Description: + Calculating force per surface area. This will be used to plot the + spanwise force distribution + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "argList.H" +#include "autoPtr.H" +#include "Time.H" +#include "timeSelector.H" +#include "TimePaths.H" +#include "ListOps.H" +#include "fvMesh.H" +#include "OFstream.H" +#include "simpleControl.H" +#include "fvOptions.H" +#include "singlePhaseTransportModel.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char* argv[]) +{ + Info << "Computing forces...." << endl; + + argList::addOption( + "patchNames", + "'(wing)'", + "List of patch names to compute"); + + argList::addOption( + "time", + "1000", + "Time instance to compute, if not provided runs all times"); + + argList::addOption( + "rho", + "1.0", + "Freestream density"); + +#include "setRootCase.H" +#include "createTime.H" + + if (!args.optionFound("time")) + { + Info << "Time not set, running all times." << endl; + } + + // Create the processor databases + autoPtr timePaths; + timePaths = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); + + List patchNames; + if (args.optionFound("patchNames")) + { + patchNames = wordReList(args.optionLookup("patchNames")()); + } + else + { + Info << "patchNames not set! Exit." << endl; + return 1; + } + + scalar rho = 1.0; + if (args.optionFound("rho")) + { + rho = readScalar(args.optionLookup("rho")()); + } + else + { + Info << "Density not set! Defaulting to 1.0." << endl; + } + + forAll(timeDirs, iTime) + { + runTime.setTime(timeDirs[iTime].value(), 0); + +#include "createMesh.H" +#include "createFields.H" + + // Initialize surface field for face-centered forces + volVectorField volumeForceField( + IOobject( + "force", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + mesh, + dimensionedVector("surfaceForce", dimensionSet(1, 1, -2, 0, 0, 0, 0), vector::zero), + fixedValueFvPatchScalarField::typeName); + + // this code is pulled from: + // src/functionObjects/forcces/forces.C + // modified slightly + vector forces(vector::zero); + + const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); + const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); + + volSymmTensorField devRhoReff( + IOobject( + IOobject::groupName("devRhoReff", U.group()), + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + -nuEff * dev(twoSymm(fvc::grad(U)))); + + const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); + + vector totalForce(vector::zero); + forAll(patchNames, cI) + { + // get the patch id label + label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); + if (patchI < 0) + { + Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; + return 1; + } + // create a shorter handle for the boundary patch + const fvPatch& patch = mesh.boundary()[patchI]; + // normal force + vectorField fN(Sfb[patchI] * p.boundaryField()[patchI] * rho); + // tangential force + vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); + // sum them up + forAll(patch, faceI) + { + // compute forces + forces.x() = fN[faceI].x() + fT[faceI].x(); + forces.y() = fN[faceI].y() + fT[faceI].y(); + forces.z() = fN[faceI].z() + fT[faceI].z(); + + // project force direction + volumeForceField.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; + + totalForce.x() += forces.x(); + totalForce.y() += forces.y(); + totalForce.z() += forces.z(); + } + } + volumeForceField.write(); + + Info << "Total force: " << totalForce << endl; + + Info << "Computing forces.... Completed!" << endl; + } + return 0; +} + +// ************************************************************************* // diff --git a/src/utilities/postProcessing/calcForceIncompressible/createFields.H b/src/utilities/postProcessing/calcForceIncompressible/createFields.H new file mode 100755 index 00000000..13f6d410 --- /dev/null +++ b/src/utilities/postProcessing/calcForceIncompressible/createFields.H @@ -0,0 +1,48 @@ +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 +); + +Info<< "Reading field nut\n" << endl; +volScalarField nut +( + IOobject + ( + "nut", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "createPhi.H" + +singlePhaseTransportModel laminarTransport(U, phi); + +volScalarField nuEff = nut + laminarTransport.nu(); + From 32066373ae0f2699a28d485ec9e04931f26a7ed2 Mon Sep 17 00:00:00 2001 From: Bernardo Pacini Date: Tue, 2 Jul 2024 08:37:51 -0400 Subject: [PATCH 2/2] Replacing previous version of forcePerS --- .../calcForceCompressible/Make/files | 3 - .../calcForceCompressible/Make/options | 21 -- .../calcForceCompressible/calcForce.C | 147 -------------- .../calcForceCompressible/createFields.H | 55 ------ .../calcForceIncompressible/Make/files | 3 - .../calcForceIncompressible/Make/options | 20 -- .../calcForceIncompressible/calcForce.C | 162 ---------------- .../calcForceIncompressible/createFields.H | 48 ----- .../calcForcePerSCompressible/calcForcePerS.C | 172 ++++++++--------- .../calcForcePerS.C | 179 +++++++++--------- 10 files changed, 177 insertions(+), 633 deletions(-) delete mode 100755 src/utilities/postProcessing/calcForceCompressible/Make/files delete mode 100755 src/utilities/postProcessing/calcForceCompressible/Make/options delete mode 100644 src/utilities/postProcessing/calcForceCompressible/calcForce.C delete mode 100755 src/utilities/postProcessing/calcForceCompressible/createFields.H delete mode 100755 src/utilities/postProcessing/calcForceIncompressible/Make/files delete mode 100755 src/utilities/postProcessing/calcForceIncompressible/Make/options delete mode 100755 src/utilities/postProcessing/calcForceIncompressible/calcForce.C delete mode 100755 src/utilities/postProcessing/calcForceIncompressible/createFields.H diff --git a/src/utilities/postProcessing/calcForceCompressible/Make/files b/src/utilities/postProcessing/calcForceCompressible/Make/files deleted file mode 100755 index 4fb73b3b..00000000 --- a/src/utilities/postProcessing/calcForceCompressible/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -calcForce.C - -EXE = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedBins/calcForceCompressible diff --git a/src/utilities/postProcessing/calcForceCompressible/Make/options b/src/utilities/postProcessing/calcForceCompressible/Make/options deleted file mode 100755 index e837721a..00000000 --- a/src/utilities/postProcessing/calcForceCompressible/Make/options +++ /dev/null @@ -1,21 +0,0 @@ -EXE_INC = \ - -std=c++11 \ - -Wno-old-style-cast \ - -Wno-conversion-null \ - -Wno-deprecated-copy \ - -I$(LIB_SRC)/transportModels/compressible/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/finiteVolume/cfdTools \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude - -EXE_LIBS = \ - -lcompressibleTransportModels \ - -lfluidThermophysicalModels \ - -lspecie \ - -lfiniteVolume \ - -lsampling \ - -lmeshTools \ - -lfvOptions - diff --git a/src/utilities/postProcessing/calcForceCompressible/calcForce.C b/src/utilities/postProcessing/calcForceCompressible/calcForce.C deleted file mode 100644 index d6a17f10..00000000 --- a/src/utilities/postProcessing/calcForceCompressible/calcForce.C +++ /dev/null @@ -1,147 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - DAFoam : Discrete Adjoint with OpenFOAM - Version : v3 - - Description: - Calculating force per surface area. This will be used to plot the - spanwise force distribution - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "argList.H" -#include "autoPtr.H" -#include "Time.H" -#include "timeSelector.H" -#include "TimePaths.H" -#include "ListOps.H" -#include "fvMesh.H" -#include "OFstream.H" -#include "simpleControl.H" -#include "fvOptions.H" -#include "fluidThermo.H" - -using namespace Foam; - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char* argv[]) -{ - Info << "Computing forces...." << endl; - - argList::addOption( - "patchNames", - "'(wing)'", - "List of patch names to compute"); - - argList::addOption( - "time", - "1000", - "Time instance to compute, if not provided runs all times"); - -#include "setRootCase.H" -#include "createTime.H" - - if (!args.optionFound("time")) - { - Info << "Time not set, running all times." << endl; - } - - // Create the processor databases - autoPtr timePaths; - timePaths = autoPtr::New(args.rootPath(), args.caseName()); - - const instantList timeDirs(timeSelector::select(timePaths->times(), args)); - - List patchNames; - if (args.optionFound("patchNames")) - { - patchNames = wordReList(args.optionLookup("patchNames")()); - } - else - { - Info << "patchNames not set! Exit." << endl; - return 1; - } - - forAll(timeDirs, iTime) - { - runTime.setTime(timeDirs[iTime].value(), 0); - -#include "createMesh.H" -#include "createFields.H" - - // Initialize surface field for face-centered forces - volVectorField volumeForceField( - IOobject( - "force", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - mesh, - dimensionedVector("surfaceForce", dimensionSet(1, 1, -2, 0, 0, 0, 0), vector::zero), - fixedValueFvPatchScalarField::typeName); - - // this code is pulled from: - // src/functionObjects/forcces/forces.C - // modified slightly - vector forces(vector::zero); - - const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); - const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); - - volSymmTensorField devRhoReff( - IOobject( - IOobject::groupName("devRhoReff", U.group()), - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - (-rho * nuEff) * dev(twoSymm(fvc::grad(U)))); - - const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); - - vector totalForce(vector::zero); - forAll(patchNames, cI) - { - // get the patch id label - label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); - if (patchI < 0) - { - Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; - return 1; - } - // create a shorter handle for the boundary patch - const fvPatch& patch = mesh.boundary()[patchI]; - // normal force - vectorField fN(Sfb[patchI] * p.boundaryField()[patchI] * rho); - // tangential force - vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); - // sum them up - forAll(patch, faceI) - { - // compute forces - forces.x() = fN[faceI].x() + fT[faceI].x(); - forces.y() = fN[faceI].y() + fT[faceI].y(); - forces.z() = fN[faceI].z() + fT[faceI].z(); - - // project force direction - volumeForceField.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; - - totalForce.x() += forces.x(); - totalForce.y() += forces.y(); - totalForce.z() += forces.z(); - } - } - volumeForceField.write(); - - Info << "Total force: " << totalForce << endl; - - Info << "Computing forces.... Completed!" << endl; - } - return 0; -} - -// ************************************************************************* // diff --git a/src/utilities/postProcessing/calcForceCompressible/createFields.H b/src/utilities/postProcessing/calcForceCompressible/createFields.H deleted file mode 100755 index 0ec75421..00000000 --- a/src/utilities/postProcessing/calcForceCompressible/createFields.H +++ /dev/null @@ -1,55 +0,0 @@ -autoPtr pThermo -( - fluidThermo::New(mesh) -); -fluidThermo& thermo = pThermo(); -thermo.validate(args.executable(), "h", "e"); - -volScalarField& p = thermo.p(); - -volScalarField rho -( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - thermo.rho() -); - -Info<< "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -Info<< "Reading field nut\n" << endl; -volScalarField nut -( - IOobject - ( - "nut", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -volScalarField nuEff = nut + thermo.nu(); - -#include "compressibleCreatePhi.H" - - diff --git a/src/utilities/postProcessing/calcForceIncompressible/Make/files b/src/utilities/postProcessing/calcForceIncompressible/Make/files deleted file mode 100755 index 5ed44cb8..00000000 --- a/src/utilities/postProcessing/calcForceIncompressible/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -calcForce.C - -EXE = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedBins/calcForceIncompressible diff --git a/src/utilities/postProcessing/calcForceIncompressible/Make/options b/src/utilities/postProcessing/calcForceIncompressible/Make/options deleted file mode 100755 index 86de1a4a..00000000 --- a/src/utilities/postProcessing/calcForceIncompressible/Make/options +++ /dev/null @@ -1,20 +0,0 @@ -EXE_INC = \ - -std=c++11 \ - -Wno-old-style-cast \ - -Wno-conversion-null \ - -Wno-deprecated-copy \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude - - -EXE_LIBS = \ - -lturbulenceModels \ - -lincompressibleTransportModels \ - -lfiniteVolume \ - -lmeshTools \ - -lfvOptions \ - -lsampling - diff --git a/src/utilities/postProcessing/calcForceIncompressible/calcForce.C b/src/utilities/postProcessing/calcForceIncompressible/calcForce.C deleted file mode 100755 index 7b58b940..00000000 --- a/src/utilities/postProcessing/calcForceIncompressible/calcForce.C +++ /dev/null @@ -1,162 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - DAFoam : Discrete Adjoint with OpenFOAM - Version : v3 - - Description: - Calculating force per surface area. This will be used to plot the - spanwise force distribution - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "argList.H" -#include "autoPtr.H" -#include "Time.H" -#include "timeSelector.H" -#include "TimePaths.H" -#include "ListOps.H" -#include "fvMesh.H" -#include "OFstream.H" -#include "simpleControl.H" -#include "fvOptions.H" -#include "singlePhaseTransportModel.H" - -using namespace Foam; - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char* argv[]) -{ - Info << "Computing forces...." << endl; - - argList::addOption( - "patchNames", - "'(wing)'", - "List of patch names to compute"); - - argList::addOption( - "time", - "1000", - "Time instance to compute, if not provided runs all times"); - - argList::addOption( - "rho", - "1.0", - "Freestream density"); - -#include "setRootCase.H" -#include "createTime.H" - - if (!args.optionFound("time")) - { - Info << "Time not set, running all times." << endl; - } - - // Create the processor databases - autoPtr timePaths; - timePaths = autoPtr::New(args.rootPath(), args.caseName()); - - const instantList timeDirs(timeSelector::select(timePaths->times(), args)); - - List patchNames; - if (args.optionFound("patchNames")) - { - patchNames = wordReList(args.optionLookup("patchNames")()); - } - else - { - Info << "patchNames not set! Exit." << endl; - return 1; - } - - scalar rho = 1.0; - if (args.optionFound("rho")) - { - rho = readScalar(args.optionLookup("rho")()); - } - else - { - Info << "Density not set! Defaulting to 1.0." << endl; - } - - forAll(timeDirs, iTime) - { - runTime.setTime(timeDirs[iTime].value(), 0); - -#include "createMesh.H" -#include "createFields.H" - - // Initialize surface field for face-centered forces - volVectorField volumeForceField( - IOobject( - "force", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - mesh, - dimensionedVector("surfaceForce", dimensionSet(1, 1, -2, 0, 0, 0, 0), vector::zero), - fixedValueFvPatchScalarField::typeName); - - // this code is pulled from: - // src/functionObjects/forcces/forces.C - // modified slightly - vector forces(vector::zero); - - const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); - const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); - - volSymmTensorField devRhoReff( - IOobject( - IOobject::groupName("devRhoReff", U.group()), - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - -nuEff * dev(twoSymm(fvc::grad(U)))); - - const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); - - vector totalForce(vector::zero); - forAll(patchNames, cI) - { - // get the patch id label - label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); - if (patchI < 0) - { - Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; - return 1; - } - // create a shorter handle for the boundary patch - const fvPatch& patch = mesh.boundary()[patchI]; - // normal force - vectorField fN(Sfb[patchI] * p.boundaryField()[patchI] * rho); - // tangential force - vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); - // sum them up - forAll(patch, faceI) - { - // compute forces - forces.x() = fN[faceI].x() + fT[faceI].x(); - forces.y() = fN[faceI].y() + fT[faceI].y(); - forces.z() = fN[faceI].z() + fT[faceI].z(); - - // project force direction - volumeForceField.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; - - totalForce.x() += forces.x(); - totalForce.y() += forces.y(); - totalForce.z() += forces.z(); - } - } - volumeForceField.write(); - - Info << "Total force: " << totalForce << endl; - - Info << "Computing forces.... Completed!" << endl; - } - return 0; -} - -// ************************************************************************* // diff --git a/src/utilities/postProcessing/calcForceIncompressible/createFields.H b/src/utilities/postProcessing/calcForceIncompressible/createFields.H deleted file mode 100755 index 13f6d410..00000000 --- a/src/utilities/postProcessing/calcForceIncompressible/createFields.H +++ /dev/null @@ -1,48 +0,0 @@ -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 -); - -Info<< "Reading field nut\n" << endl; -volScalarField nut -( - IOobject - ( - "nut", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "createPhi.H" - -singlePhaseTransportModel laminarTransport(U, phi); - -volScalarField nuEff = nut + laminarTransport.nu(); - diff --git a/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C b/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C index d2021d0d..750283ff 100755 --- a/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C +++ b/src/utilities/postProcessing/calcForcePerSCompressible/calcForcePerS.C @@ -11,7 +11,11 @@ #include "fvCFD.H" #include "argList.H" +#include "autoPtr.H" #include "Time.H" +#include "timeSelector.H" +#include "TimePaths.H" +#include "ListOps.H" #include "fvMesh.H" #include "OFstream.H" #include "simpleControl.H" @@ -24,40 +28,31 @@ using namespace Foam; int main(int argc, char* argv[]) { - Info << "Computing forcePerS...." << endl; + Info << "Computing forces...." << endl; argList::addOption( "patchNames", - "'(inlet)'", + "'(wing)'", "List of patch names to compute"); - argList::addOption( - "forceDir", - "'(0 0 1)'", - "Force direction"); - argList::addOption( "time", "1000", - "Tme instance to compute"); + "Time instance to compute, if not provided runs all times"); #include "setRootCase.H" #include "createTime.H" - scalar time; - if (args.optionFound("time")) + if (!args.optionFound("time")) { - time = readScalar(args.optionLookup("time")()); + Info << "Time not set, running all times." << endl; } - else - { - Info << "time not set! Exit." << endl; - return 1; - } - runTime.setTime(time, 0); -#include "createMesh.H" -#include "createFields.H" + // Create the processor databases + autoPtr timePaths; + timePaths = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); List patchNames; if (args.optionFound("patchNames")) @@ -70,83 +65,82 @@ int main(int argc, char* argv[]) return 1; } - List forceDir1; - if (args.optionFound("forceDir")) - { - forceDir1 = scalarList(args.optionLookup("forceDir")()); - } - else + forAll(timeDirs, iTime) { - Info << "forceDir not set! Exit." << endl; - return 1; - } - vector forceDir(vector::zero); - forceDir.x() = forceDir1[0]; - forceDir.y() = forceDir1[1]; - forceDir.z() = forceDir1[2]; - - volScalarField forcePerS( - IOobject( - "forcePerS", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - mesh, - dimensionedScalar("forcePerS", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0), - fixedValueFvPatchScalarField::typeName); - - // this code is pulled from: - // src/functionObjects/forcces/forces.C - // modified slightly - vector forces(vector::zero); - - const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); - const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); - - volSymmTensorField devRhoReff( - IOobject( - IOobject::groupName("devRhoReff", U.group()), - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - (-rho * nuEff) * dev(twoSymm(fvc::grad(U)))); + runTime.setTime(timeDirs[iTime].value(), 0); - const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); +#include "createMesh.H" +#include "createFields.H" - forAll(patchNames, cI) - { - // get the patch id label - label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); - if (patchI < 0) - { - Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; - return 1; - } - // create a shorter handle for the boundary patch - const fvPatch& patch = mesh.boundary()[patchI]; - // normal force - vectorField fN( - Sfb[patchI] * p.boundaryField()[patchI]); - // tangential force - vectorField fT(Sfb[patchI] & devRhoReffb[patchI]); - // sum them up - forAll(patch, faceI) + // Initialize surface field for face-centered forces + volVectorField forcePerS( + IOobject( + "forcePerS", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + mesh, + dimensionedVector("forcePerS", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector::zero), + fixedValueFvPatchScalarField::typeName); + + // this code is pulled from: + // src/functionObjects/forcces/forces.C + // modified slightly + vector forces(vector::zero); + + const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); + const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); + + volSymmTensorField devRhoReff( + IOobject( + IOobject::groupName("devRhoReff", U.group()), + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + (-rho * nuEff) * dev(twoSymm(fvc::grad(U)))); + + const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); + + vector totalForce(vector::zero); + forAll(patchNames, cI) { - forces.x() = fN[faceI].x() + fT[faceI].x(); - forces.y() = fN[faceI].y() + fT[faceI].y(); - forces.z() = fN[faceI].z() + fT[faceI].z(); - scalar force = forces & forceDir; - forcePerS.boundaryFieldRef()[patchI][faceI] = force / magSfb[patchI][faceI]; + // get the patch id label + label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); + if (patchI < 0) + { + Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; + return 1; + } + // create a shorter handle for the boundary patch + const fvPatch& patch = mesh.boundary()[patchI]; + // normal force + vectorField fN(Sfb[patchI] * p.boundaryField()[patchI]); + // tangential force + vectorField fT(Sfb[patchI] & devRhoReffb[patchI]); + // sum them up + forAll(patch, faceI) + { + // compute forces + forces.x() = fN[faceI].x() + fT[faceI].x(); + forces.y() = fN[faceI].y() + fT[faceI].y(); + forces.z() = fN[faceI].z() + fT[faceI].z(); + + // project force direction + forcePerS.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; + + totalForce.x() += forces.x(); + totalForce.y() += forces.y(); + totalForce.z() += forces.z(); + } } - } - forcePerS.write(); - - Info << "Force: " << forces << endl; + forcePerS.write(); - Info << "Computing forcePerS.... Completed!" << endl; + Info << "Total force: " << totalForce << endl; + Info << "Computing forces.... Completed!" << endl; + } return 0; } diff --git a/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C b/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C index f004e293..781891b2 100755 --- a/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C +++ b/src/utilities/postProcessing/calcForcePerSIncompressible/calcForcePerS.C @@ -11,7 +11,11 @@ #include "fvCFD.H" #include "argList.H" +#include "autoPtr.H" #include "Time.H" +#include "timeSelector.H" +#include "TimePaths.H" +#include "ListOps.H" #include "fvMesh.H" #include "OFstream.H" #include "simpleControl.H" @@ -24,40 +28,36 @@ using namespace Foam; int main(int argc, char* argv[]) { - Info << "Computing forcePerS...." << endl; + Info << "Computing forces...." << endl; argList::addOption( "patchNames", - "'(inlet)'", + "'(wing)'", "List of patch names to compute"); - argList::addOption( - "forceDir", - "'(0 0 1)'", - "Force direction"); - argList::addOption( "time", "1000", - "Tme instance to compute"); + "Time instance to compute, if not provided runs all times"); + + argList::addOption( + "rho", + "1.0", + "Freestream density"); #include "setRootCase.H" #include "createTime.H" - scalar time; - if (args.optionFound("time")) + if (!args.optionFound("time")) { - time = readScalar(args.optionLookup("time")()); + Info << "Time not set, running all times." << endl; } - else - { - Info << "time not set! Exit." << endl; - return 1; - } - runTime.setTime(time, 0); -#include "createMesh.H" -#include "createFields.H" + // Create the processor databases + autoPtr timePaths; + timePaths = autoPtr::New(args.rootPath(), args.caseName()); + + const instantList timeDirs(timeSelector::select(timePaths->times(), args)); List patchNames; if (args.optionFound("patchNames")) @@ -70,83 +70,92 @@ int main(int argc, char* argv[]) return 1; } - List forceDir1; - if (args.optionFound("forceDir")) + scalar rho = 1.0; + if (args.optionFound("rho")) { - forceDir1 = scalarList(args.optionLookup("forceDir")()); + rho = readScalar(args.optionLookup("rho")()); } else { - Info << "forceDir not set! Exit." << endl; - return 1; + Info << "Density not set! Defaulting to 1.0." << endl; } - vector forceDir(vector::zero); - forceDir.x() = forceDir1[0]; - forceDir.y() = forceDir1[1]; - forceDir.z() = forceDir1[2]; - - volScalarField forcePerS( - IOobject( - "forcePerS", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - mesh, - dimensionedScalar("forcePerS", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0), - fixedValueFvPatchScalarField::typeName); - - // this code is pulled from: - // src/functionObjects/forcces/forces.C - // modified slightly - vector forces(vector::zero); - - const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); - const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); - - volSymmTensorField devRhoReff( - IOobject( - IOobject::groupName("devRhoReff", U.group()), - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE), - -nuEff * dev(twoSymm(fvc::grad(U)))); - - const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); - forAll(patchNames, cI) + forAll(timeDirs, iTime) { - // get the patch id label - label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); - if (patchI < 0) - { - Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; - return 1; - } - // create a shorter handle for the boundary patch - const fvPatch& patch = mesh.boundary()[patchI]; - // normal force - vectorField fN( - Sfb[patchI] * p.boundaryField()[patchI]); - // tangential force - vectorField fT(Sfb[patchI] & devRhoReffb[patchI]); - // sum them up - forAll(patch, faceI) + runTime.setTime(timeDirs[iTime].value(), 0); + +#include "createMesh.H" +#include "createFields.H" + + // Initialize surface field for face-centered forces + volVectorField forcePerS( + IOobject( + "forcePerS", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + mesh, + dimensionedVector("forcePerS", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector::zero), + fixedValueFvPatchScalarField::typeName); + + // this code is pulled from: + // src/functionObjects/forcces/forces.C + // modified slightly + vector forces(vector::zero); + + const surfaceVectorField::Boundary& Sfb = mesh.Sf().boundaryField(); + const surfaceScalarField::Boundary& magSfb = mesh.magSf().boundaryField(); + + volSymmTensorField devRhoReff( + IOobject( + IOobject::groupName("devRhoReff", U.group()), + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE), + -nuEff * dev(twoSymm(fvc::grad(U)))); + + const volSymmTensorField::Boundary& devRhoReffb = devRhoReff.boundaryField(); + + vector totalForce(vector::zero); + forAll(patchNames, cI) { - forces.x() = fN[faceI].x() + fT[faceI].x(); - forces.y() = fN[faceI].y() + fT[faceI].y(); - forces.z() = fN[faceI].z() + fT[faceI].z(); - scalar force = forces & forceDir; - forcePerS.boundaryFieldRef()[patchI][faceI] = force / magSfb[patchI][faceI]; + // get the patch id label + label patchI = mesh.boundaryMesh().findPatchID(patchNames[cI]); + if (patchI < 0) + { + Info << "ERROR: Patch name " << patchNames[cI] << " not found in constant/polyMesh/boundary! Exit!" << endl; + return 1; + } + // create a shorter handle for the boundary patch + const fvPatch& patch = mesh.boundary()[patchI]; + // normal force + vectorField fN(Sfb[patchI] * p.boundaryField()[patchI] * rho); + // tangential force + vectorField fT(Sfb[patchI] & devRhoReffb[patchI] * rho); + // sum them up + forAll(patch, faceI) + { + // compute forces + forces.x() = fN[faceI].x() + fT[faceI].x(); + forces.y() = fN[faceI].y() + fT[faceI].y(); + forces.z() = fN[faceI].z() + fT[faceI].z(); + + // project force direction + forcePerS.boundaryFieldRef()[patchI][faceI] = forces / magSfb[patchI][faceI]; + + totalForce.x() += forces.x(); + totalForce.y() += forces.y(); + totalForce.z() += forces.z(); + } } - } - forcePerS.write(); - - Info << "Force: " << forces << endl; + forcePerS.write(); - Info << "Computing forcePerS.... Completed!" << endl; + Info << "Total force: " << totalForce << endl; + Info << "Computing forces.... Completed!" << endl; + } return 0; }