diff --git a/etc/caseDicts/postProcessing/fields/energySpectrum b/etc/caseDicts/postProcessing/fields/energySpectrum new file mode 100644 index 0000000000000000000000000000000000000000..4df60270faf7a9a2f2de7469183057feaa63d5b8 --- /dev/null +++ b/etc/caseDicts/postProcessing/fields/energySpectrum @@ -0,0 +1,19 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Web: www.OpenFOAM.com + \\/ M anipulation | +------------------------------------------------------------------------------- +Description + Calculates the energy spectrum for a box of turbulence + +\*---------------------------------------------------------------------------*/ + +type energySpectrum; +libs ("librandomProcessesFunctionObjects.so"); + +executeControl writeTime; +writeControl writeTime; + +// ************************************************************************* // diff --git a/src/functionObjects/randomProcesses/Allwmake b/src/functionObjects/randomProcesses/Allwmake new file mode 100755 index 0000000000000000000000000000000000000000..c81ffedf7817dc95192fe7649ba032bf6e7488c9 --- /dev/null +++ b/src/functionObjects/randomProcesses/Allwmake @@ -0,0 +1,16 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory +targetType=libso # Preferred library type +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments +. $WM_PROJECT_DIR/wmake/scripts/have_fftw + +#------------------------------------------------------------------------------ + +if have_fftw +then + wmake $targetType +else + echo "==> skip randomProcesses library (no FFTW)" +fi + +#------------------------------------------------------------------------------ diff --git a/src/functionObjects/randomProcesses/Make/files b/src/functionObjects/randomProcesses/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..e13e8afc998f172001da30bccfe8e7c54013061a --- /dev/null +++ b/src/functionObjects/randomProcesses/Make/files @@ -0,0 +1,3 @@ +energySpectrum/energySpectrum.C + +LIB = $(FOAM_LIBBIN)/librandomProcessesFunctionObjects diff --git a/src/functionObjects/randomProcesses/Make/options b/src/functionObjects/randomProcesses/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..4c713c968a763f22cc0bbc61b22ea179f8b56466 --- /dev/null +++ b/src/functionObjects/randomProcesses/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/randomProcesses/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +LIB_LIBS = \ + -lrandomProcesses \ + -lfiniteVolume \ + -lmeshTools diff --git a/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C new file mode 100644 index 0000000000000000000000000000000000000000..7c19426a6e998906e4f62f40f96cc99f05bce937 --- /dev/null +++ b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.C @@ -0,0 +1,263 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "energySpectrum.H" +#include "fft.H" +#include "fvMesh.H" +#include "boundBox.H" +#include "OFstream.H" +#include "mathematicalConstants.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + defineTypeNameAndDebug(energySpectrum, 0); + addToRunTimeSelectionTable(functionObject, energySpectrum, dictionary); +} +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::functionObjects::energySpectrum::writeFileHeader(Ostream& os) +{ + writeHeader(os, "Turbulence energy spectra"); + + writeCommented(os, "kappa E(kappa)"); + + os << endl; +} + + +void Foam::functionObjects::energySpectrum::calcAndWriteSpectrum +( + const vectorField& U, + const vectorField& C, + const vector& c0, + const vector& deltaC, + const Vector<label>& N, + const scalar kappaNorm +) +{ + Log << "Computing FFT" << endl; + const complexVectorField Uf + ( + fft::forwardTransform + ( + ReComplexField(U), + List<label>({N.x(), N.y(), N.z()}) + ) + /scalar(cmptProduct(N)) + ); + + + Log << "Computing wave numbers" << endl; + const label nMax = cmptMax(N); + scalarField kappa(nMax); + forAll(kappa, kappai) + { + kappa[kappai] = kappai*kappaNorm; + } + + Log << "Computing energy spectrum" << endl; + scalarField E(nMax, 0); + const scalarField Ec(0.5*magSqr(Uf)); + forAll(C, celli) + { + point Cc(C[celli] - c0); + label i = round((Cc.x() - 0.5*deltaC.x())/(deltaC.x())*(N.x() - 1)); + label j = round((Cc.y() - 0.5*deltaC.y())/(deltaC.y())*(N.y() - 1)); + label k = round((Cc.z() - 0.5*deltaC.z())/(deltaC.z())*(N.z() - 1)); + label kappai = round(Foam::sqrt(scalar(i*i + j*j + k*k))); + + E[kappai] += Ec[celli]; + } + + E /= kappaNorm; + + Log << "Writing spectrum" << endl; + autoPtr<OFstream> osPtr = createFile(name(), time_.value()); + OFstream& os = osPtr.ref(); + writeFileHeader(os); + + forAll(kappa, kappai) + { + os << kappa[kappai] << tab << E[kappai] << nl; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::functionObjects::energySpectrum::energySpectrum +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + fvMeshFunctionObject(name, runTime, dict), + writeFile(mesh_, name), + cellAddr_(mesh_.nCells()), + UName_("U"), + N_(Zero), + c0_(Zero), + deltaC_(Zero), + kappaNorm_(0) +{ + read(dict); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::functionObjects::energySpectrum::read(const dictionary& dict) +{ + fvMeshFunctionObject::read(dict); + writeFile::read(dict); + + dict.readIfPresent("U", UName_); + + const boundBox meshBb(mesh_.bounds()); + + // Assume all cells are the same size... + const cell& c = mesh_.cells()[0]; + boundBox cellBb(boundBox::invertedBox); + forAll(c, facei) + { + const face& f = mesh_.faces()[c[facei]]; + cellBb.add(mesh_.points(), f); + } + + const vector L(meshBb.max() - meshBb.min()); + const vector nCellXYZ(cmptDivide(L, cellBb.max() - cellBb.min())); + + N_ = Vector<label> + ( + round(nCellXYZ.x()), + round(nCellXYZ.z()), + round(nCellXYZ.z()) + ); + + // Check that the mesh is a structured box + vector cellDx(cellBb.max() - cellBb.min()); + vector expectedMax(N_.x()*cellDx.x(), N_.y()*cellDx.y(), N_.z()*cellDx.z()); + vector relativeSize(cmptDivide(L, expectedMax)); + for (direction i = 0; i < 3; ++i) + { + if (mag(relativeSize[i] - 1) > 1e-3) + { + FatalErrorInFunction + << name() << " function object is only appropriate for " + << "isotropic structured IJK meshes. Mesh extents: " << L + << ", computed IJK mesh extents: " << expectedMax + << exit(FatalError); + } + } + Log << "Mesh extents (deltax,deltay,deltaz): " << L << endl; + Log << "Number of cells (Nx,Ny,Nz): " << N_ << endl; + + // Map into i-j-k co-ordinates + const vectorField& C = mesh_.C(); + c0_ = returnReduce(min(C), minOp<vector>()); + const vector cMax = returnReduce(max(C), maxOp<vector>()); + deltaC_ = cMax - c0_; + + forAll(C, celli) + { + label i = round((C[celli].x() - c0_.x())/(deltaC_.x())*(N_.x() - 1)); + label j = round((C[celli].y() - c0_.y())/(deltaC_.y())*(N_.y() - 1)); + label k = round((C[celli].z() - c0_.z())/(deltaC_.z())*(N_.z() - 1)); + + cellAddr_[celli] = k + j*N_.y() + i*N_.y()*N_.z(); + } + + kappaNorm_ = constant::mathematical::twoPi/cmptMax(L); + + return true; +} + + +bool Foam::functionObjects::energySpectrum::execute() +{ + return true; +} + + +bool Foam::functionObjects::energySpectrum::write() +{ + const auto& U = mesh_.lookupObject<volVectorField>(UName_); + const vectorField& Uc = U.primitiveField(); + const vectorField& C = mesh_.C(); + + if (Pstream::parRun()) + { + PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking); + + UOPstream toProc(0, pBufs); + + toProc << Uc << C << cellAddr_; + + pBufs.finishedSends(); + + if (Pstream::master()) + { + const label nGlobalCells(cmptProduct(N_)); + vectorField Uijk(nGlobalCells); + vectorField Cijk(nGlobalCells); + + for (label proci = 0; proci < Pstream::nProcs(); ++proci) + { + UIPstream fromProc(proci, pBufs); + vectorField Up; + vectorField Cp; + labelList cellAddrp; + fromProc >> Up >> Cp >> cellAddrp; + UIndirectList<vector>(Uijk, cellAddrp) = Up; + UIndirectList<vector>(Cijk, cellAddrp) = Cp; + } + + calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_); + } + + } + else + { + vectorField Uijk(Uc, cellAddr_); + vectorField Cijk(C, cellAddr_); + + calcAndWriteSpectrum(Uijk, Cijk, c0_, deltaC_, N_, kappaNorm_); + } + + return true; +} + + +// ************************************************************************* // diff --git a/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.H b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.H new file mode 100644 index 0000000000000000000000000000000000000000..9958d6edbe7f412a74defc63adb2a9bcba0a8f7d --- /dev/null +++ b/src/functionObjects/randomProcesses/energySpectrum/energySpectrum.H @@ -0,0 +1,178 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +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::functionObjects::energySpectrum + +Group + grpFieldFunctionObjects + +Description + Calculates the energy spectrum for a structured IJK mesh + +Usage + Example of function object specification: + \verbatim + energySpectrum1 + { + type energySpectrum; + libs ("libfieldFunctionObjects.so"); + ... + log yes; + } + \endverbatim + + Where the entries comprise: + \table + Property | Description | Required | Default value + type | type name: energySpectrum | yes | + log | write info to standard output | no | yes + \endtable + + Output data is written to the file \<timeDir\>/energySpectrum.dat + +See also + Foam::functionObjects::fvMeshFunctionObject + Foam::functionObjects::writeFile + +SourceFiles + energySpectrum.C + energySpectrumTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef functionObjects_energySpectrum_H +#define functionObjects_energySpectrum_H + +#include "fvMeshFunctionObject.H" +#include "writeFile.H" +#include "Vector.H" +#include "vectorField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + +/*---------------------------------------------------------------------------*\ + Class energySpectrum Declaration +\*---------------------------------------------------------------------------*/ + +class energySpectrum +: + public fvMeshFunctionObject, + public writeFile +{ +protected: + + // Protected data + + //- I-J-K mesh addressing + labelList cellAddr_; + + //- Name of velocity field, default = U + word UName_; + + //- Number of cells in I-J-K directions + Vector<label> N_; + + //- Reference point + vector c0_; + + //- Cell length scale + vector deltaC_; + + //- Wave number + scalar kappaNorm_; + + + // Protected Member Functions + + + //- Output file header information + virtual void writeFileHeader(Ostream& os); + + //- Calculate and write the spectrum + void calcAndWriteSpectrum + ( + const vectorField& U, + const vectorField& C, + const vector& c0, + const vector& deltaC, + const Vector<label>& N, + const scalar kappaNorm + ); + + //- No copy construct + energySpectrum(const energySpectrum&) = delete; + + //- No copy assignment + void operator=(const energySpectrum&) = delete; + + +public: + + //- Runtime type information + TypeName("energySpectrum"); + + + // Constructors + + //- Construct from Time and dictionary + energySpectrum + ( + const word& name, + const Time& runTime, + const dictionary& dict + ); + + + //- Destructor + virtual ~energySpectrum() = default; + + + // Member Functions + + //- Read the field min/max data + virtual bool read(const dictionary&); + + //- Execute, currently does nothing + virtual bool execute(); + + //- Write the energySpectrum + virtual bool write(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace functionObjects +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //