From ceed53775d5fcf96e207252f1fa39ff90ff73027 Mon Sep 17 00:00:00 2001 From: Andrew Heather <> Date: Fri, 5 Jun 2020 14:55:45 +0100 Subject: [PATCH] ENH: Updated Curle function object The function object now computes the acoustic pressure at a list of user specified locations, or from the face centres of a user-supplied surface. When operating on an input surface, the output can be written back to the surface or as a list of point values. Example of function object specification: Curle1 { type Curle; libs ("libfieldFunctionObjects.so"); ... patches (surface1 surface2); c0 330; // Input - either points or surface input points; observerPositions ((0 0 0)(1 0 0)); //input surface; //surface "inputSurface.obj" // Output - either points or surface output points; //output surface; //surfaceType ensight; } Where the entries comprise: Property | Description | Required | Default value type | Type name: Curle | yes | p | Pressure field name | no | p patches | Sound generation patch names | yes | c0 | Reference speed of sound | yes | input | Input type | yes | observerPositions | List of observer positions (x y z) | no | surface | Input surface file name | no | output | Output type | yes | surfaceType | Output surface type | no | --- src/functionObjects/field/Curle/Curle.C | 301 ++++++++++++++++-------- src/functionObjects/field/Curle/Curle.H | 103 +++++--- 2 files changed, 270 insertions(+), 134 deletions(-) diff --git a/src/functionObjects/field/Curle/Curle.C b/src/functionObjects/field/Curle/Curle.C index d3a16edbb5..b249470668 100644 --- a/src/functionObjects/field/Curle/Curle.C +++ b/src/functionObjects/field/Curle/Curle.C @@ -43,57 +43,12 @@ namespace functionObjects } } - -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // - -bool Foam::functionObjects::Curle::calc() -{ - if (foundObject(fieldName_)) - { - // Evaluate pressure force time derivative - - const volScalarField& p = lookupObject(fieldName_); - const volScalarField dpdt(scopedName("dpdt"), fvc::ddt(p)); - const volScalarField::Boundary& dpdtBf = dpdt.boundaryField(); - const surfaceVectorField::Boundary& SfBf = mesh_.Sf().boundaryField(); - - dimensionedVector dfdt("dfdt", p.dimensions()*dimArea/dimTime, Zero); - - for (const label patchi : patchSet_) - { - dfdt.value() += sum(dpdtBf[patchi]*SfBf[patchi]); - } - - reduce(dfdt.value(), sumOp()); - - - // Construct and store Curle acoustic pressure - - const volVectorField& C = mesh_.C(); - - auto tpDash = tmp::New - ( - IOobject - ( - resultName_, - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh_, - dimensionedScalar(p.dimensions(), Zero) - ); - auto& pDash = tpDash.ref(); - - const volVectorField d(scopedName("d"), C - x0_); - pDash = (d/magSqr(d) & dfdt)/(4.0*mathematical::pi*c0_); - - return store(resultName_, tpDash); - } - - return false; -} +const Foam::Enum +Foam::functionObjects::Curle::modeTypeNames_ +({ + { modeType::point, "point" }, + { modeType::surface, "surface" }, +}); // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -105,14 +60,17 @@ Foam::functionObjects::Curle::Curle const dictionary& dict ) : - fieldExpression(name, runTime, dict, "p"), + fvMeshFunctionObject(name, runTime, dict), + writeFile(mesh_, name), + pName_("p"), patchSet_(), - x0_("x0", dimLength, Zero), - c0_("c0", dimVelocity, Zero) + observerPositions_(), + c0_(0), + rawFilePtrs_(), + inputSurface_(), + surfaceWriterPtr_(nullptr) { read(dict); - - setResultName(typeName, fieldName_); } @@ -120,57 +78,192 @@ Foam::functionObjects::Curle::Curle bool Foam::functionObjects::Curle::read(const dictionary& dict) { - if (fieldExpression::read(dict)) + if (!(fvMeshFunctionObject::read(dict) && writeFile::read(dict))) { - patchSet_ = - mesh_.boundaryMesh().patchSet - ( - dict.get("patches") - ); - - if (patchSet_.empty()) - { - WarningInFunction - << "No patches defined" - << endl; - - return false; - } - - // Read the reference speed of sound - dict.readEntry("c0", c0_); - - if (c0_.value() < VSMALL) - { - FatalErrorInFunction - << "Reference speed of sound = " << c0_ - << " cannot be negative or zero." - << abort(FatalError); - } - - // Set the location of the effective point source to the area-average - // of the patch face centres - const volVectorField::Boundary& Cbf = mesh_.C().boundaryField(); - const surfaceScalarField::Boundary& magSfBf = - mesh_.magSf().boundaryField(); - - x0_.value() = vector::zero; - scalar sumMagSf = 0; - for (auto patchi : patchSet_) - { - x0_.value() += sum(Cbf[patchi]*magSfBf[patchi]); - sumMagSf += sum(magSfBf[patchi]); - } - - reduce(x0_.value(), sumOp()); - reduce(sumMagSf, sumOp()); - - x0_.value() /= sumMagSf + ROOTVSMALL; - - return true; + return false; } - return false; + dict.readIfPresent("p", pName_); + + patchSet_ = mesh_.boundaryMesh().patchSet(dict.get("patches")); + + if (patchSet_.empty()) + { + WarningInFunction + << "No patches defined" + << endl; + + return false; + } + + const modeType inputMode = modeTypeNames_.get("input", dict); + + switch (inputMode) + { + case modeType::point: + { + observerPositions_ = dict.get>("observerPositions"); + break; + } + case modeType::surface: + { + const fileName fName = (dict.get("surface")).expand(); + inputSurface_ = MeshedSurface(fName); + + observerPositions_ = inputSurface_.Cf(); + } + } + + if (observerPositions_.empty()) + { + WarningInFunction + << "No observer positions defined" + << endl; + + return false; + } + + const auto outputMode = modeTypeNames_.get("output", dict); + + switch (outputMode) + { + case modeType::point: + { + rawFilePtrs_.setSize(observerPositions_.size()); + forAll(rawFilePtrs_, filei) + { + rawFilePtrs_.set + ( + filei, + createFile("observer" + Foam::name(filei)) + ); + + if (rawFilePtrs_.set(filei)) + { + OFstream& os = rawFilePtrs_[filei]; + + writeHeaderValue(os, "Observer", filei); + writeHeaderValue(os, "Position", observerPositions_[filei]); + writeCommented(os, "Time"); + writeTabbed(os, "p(Curle)"); + os << endl; + } + } + break; + } + case modeType::surface: + { + if (inputMode != modeType::surface) + { + FatalIOErrorInFunction(dict) + << "Surface output is only available when input mode is " + << "type '" << modeTypeNames_[modeType::surface] << "'" + << abort(FatalIOError); + } + + const word surfaceType(dict.get("surfaceType")); + surfaceWriterPtr_ = surfaceWriter::New + ( + surfaceType, + dict.subOrEmptyDict("formatOptions").subOrEmptyDict(surfaceType) + ); + + // Use outputDir/TIME/surface-name + surfaceWriterPtr_->useTimeDir() = true; + } + } + + // Read the reference speed of sound + dict.readEntry("c0", c0_); + + if (c0_ < VSMALL) + { + FatalErrorInFunction + << "Reference speed of sound = " << c0_ + << " cannot be negative or zero." + << abort(FatalError); + } + + return true; +} + + +bool Foam::functionObjects::Curle::execute() +{ + if (!foundObject(pName_)) + { + return false; + } + + const volScalarField& p = lookupObject(pName_); + const volScalarField::Boundary& pBf = p.boundaryField(); + const volScalarField dpdt(scopedName("dpdt"), fvc::ddt(p)); + const volScalarField::Boundary& dpdtBf = dpdt.boundaryField(); + const surfaceVectorField::Boundary& SfBf = mesh_.Sf().boundaryField(); + const surfaceVectorField::Boundary& CfBf = mesh_.Cf().boundaryField(); + + scalarField pDash(observerPositions_.size(), 0); + + for (auto patchi : patchSet_) + { + const scalarField& pp = pBf[patchi]; + const scalarField& dpdtp = dpdtBf[patchi]; + const vectorField& Sfp = SfBf[patchi]; + const vectorField& Cfp = CfBf[patchi]; + + forAll(observerPositions_, pointi) + { + const vectorField r(Cfp - observerPositions_[pointi]); + const scalarField invMagR(1/(mag(r) + ROOTVSMALL)); + + pDash[pointi] += + sum((pp*sqr(invMagR) + invMagR/c0_*dpdtp)*(Sfp & r)); + } + } + Pstream::listCombineGather(pDash, plusEqOp()); + Pstream::listCombineScatter(pDash); + + if (surfaceWriterPtr_) + { + if (Pstream::master()) + { + // Time-aware, with time spliced into the output path + surfaceWriterPtr_->beginTime(time_); + + surfaceWriterPtr_->open + ( + inputSurface_.points(), + inputSurface_.surfFaces(), + (baseFileDir()/name()/"surface"), + false // serial - already merged + ); + + surfaceWriterPtr_->write("p(Curle)", pDash); + + surfaceWriterPtr_->endTime(); + surfaceWriterPtr_->clear(); + } + } + else + { + forAll(observerPositions_, pointi) + { + if (rawFilePtrs_.set(pointi)) + { + OFstream& os = rawFilePtrs_[pointi]; + writeCurrentTime(os); + os << pDash[pointi] << endl; + } + } + } + + return true; +} + + +bool Foam::functionObjects::Curle::write() +{ + return true; } diff --git a/src/functionObjects/field/Curle/Curle.H b/src/functionObjects/field/Curle/Curle.H index 2ac98ec05d..1050440769 100644 --- a/src/functionObjects/field/Curle/Curle.H +++ b/src/functionObjects/field/Curle/Curle.H @@ -35,15 +35,21 @@ Description Curle's analogy is implemented as: \f[ - p' = \frac{1}{4 \pi c_0}\frac{\vec d}{|\vec d|^2}\cdot\frac{d\vec F}{dt} + p' = \frac{1}{4 \pi} + \frac{\vec{r}}{| \vec{r} | ^2} \cdot + \left( + \frac{\vec{F}}{| \vec{r} |} + + \frac{1}{c_0}\frac{d(\vec{F})}{d(t)} + \right) \f] where \vartable - p' | Curle's acoustic pressure [Pa] or [Pa \f$(m^3/\rho)\f$] - c_0 | Reference speed of sound [m/s] - \vec d | Distance vector to observer locations [m] - \vec F | Force [N] or [N (\f$m^3/\rho\f$)] + p' | Curle's acoustic pressure [Pa] or [Pa (m3/kg)] + c_0 | Reference speed of sound [m/s] + \vec{r} | Distance vector to observer locations [m] + \vec{F} | Force [N] or [N (m3/kg)] + t | time [s] \endvartable Operands: @@ -62,16 +68,27 @@ Usage \verbatim Curle1 { - // Mandatory entries (unmodifiable) type Curle; - libs (fieldFunctionObjects); - - // Mandatory entries (runtime modifiable) - patches ( ... ) - c0 343; - - // Optional (inherited) entries + libs ("libfieldFunctionObjects.so"); ... + patches (surface1 surface2); + c0 330; + + + // Input - either points or surface + + input points; + observerPositions ((0 0 0)(1 0 0)); + + //input surface; + //surface "inputSurface.obj" + + + // Output - either points or surface + output points; + + //output surface; + //surfaceType ensight; } \endverbatim @@ -80,8 +97,14 @@ Usage Property | Description | Type | Req'd | Dflt type | Type name: Curle | word | yes | - libs | Library name: fieldFunctionObjects | word | yes | - - patches | Names of the operand patches | wordList | yes | - - c0 | Reference speed of sound [m/s] | scalar | yes | - + p | Pressure field name | word | no | p + patches | Sound generation patch names | wordList | yes | - + c0 | Reference speed of sound | scalar | yes | - + input | Input type | word | yes | - + observerPositions | List of observer positions (x y z) | pointList | no |- + surface | Input surface file name | word | no | - + output | Output type | word | yes | - + surfaceType | Output surface type | word | no | - \endtable The inherited entries are elaborated in: @@ -93,7 +116,7 @@ Usage See also - Foam::functionObject - Foam::functionObjects::fvMeshFunctionObject - - Foam::functionObjects::fieldExpression + - Foam::functionObjects::writeFile - ExtendedCodeGuide::functionObjects::field::Curle SourceFiles @@ -104,9 +127,11 @@ SourceFiles #ifndef functionObjects_Curle_H #define functionObjects_Curle_H -#include "fieldExpression.H" -#include "dimensionedScalar.H" -#include "dimensionedVector.H" +#include "fvMeshFunctionObject.H" +#include "writeFile.H" +#include "Enum.H" +#include "MeshedSurface.H" +#include "surfaceWriter.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -121,26 +146,40 @@ namespace functionObjects class Curle : - public fieldExpression + public fvMeshFunctionObject, + public writeFile { + enum class modeType + { + point, + surface + }; + + static const Enum modeTypeNames_; + + // Private Data - // Read from dictionary + //- Name of pressure field; default = p + word pName_; - //- Patches to integrate forces over - labelHashSet patchSet_; + //- Patches to integrate forces over + labelHashSet patchSet_; - //- Area-averaged centre of patch faces - dimensionedVector x0_; + //- Observer positions + List observerPositions_; - //- Reference speed of sound - dimensionedScalar c0_; + //- Reference speed of sound + scalar c0_; + //- Output files when sampling points + PtrList rawFilePtrs_; -protected: + //- Input surface when sampling a surface + MeshedSurface inputSurface_; - //- Calculate the acoustic pressure field and return true if successful - virtual bool calc(); + //- Ouput surface when sampling a surface + autoPtr surfaceWriterPtr_; public: @@ -174,6 +213,10 @@ public: //- Read the Curle data virtual bool read(const dictionary&); + + virtual bool execute(); + + virtual bool write(); };