|
[Sponsors] |
July 31, 2013, 12:46 |
Mass weighted average
|
#1 |
Member
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 16 |
Hiho,
I´m modifying the tool patchAverage with the intention to calculate the mass weighted average, but I get an error while compiling the code. The reason is the "&" operator, which should make a inner product of the patch normal vector and the velocity vector at the boundary. Here is the code (the line with the vector operation is marked in red): Code:
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright held by original author \\/ 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 2 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, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Application massAverage Description Calculates the average of the specified field over the specified patch. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: int main(int argc, char *argv[]) { timeSelector::addOptions(); argList::validArgs.append("fieldName"); argList::validArgs.append("patchName"); # include "setRootCase.H" # include "createTime.H" instantList timeDirs = timeSelector::select0(runTime, args); # include "createMesh.H" word fieldName(args.additionalArgs()[0]); word patchName(args.additionalArgs()[1]); //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),args.additionalArgs()[2]); //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),1.226); forAll(timeDirs, timeI) { runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; IOobject fieldHeader ( fieldName, runTime.timeName(), mesh, IOobject::MUST_READ ); IOobject U_header ( "U", runTime.timeName(), mesh, IOobject::MUST_READ ); // Check field exists if (fieldHeader.headerOk()) { mesh.readUpdate(); label patchi = mesh.boundaryMesh().findPatchID(patchName); if (patchi < 0) { FatalError << "Unable to find patch " << patchName << nl << exit(FatalError); } if (fieldHeader.headerClassName() == "volScalarField") { Info<< " Reading volScalarField " << fieldName << endl; volScalarField field(fieldHeader, mesh); volVectorField U(U_header, mesh); volScalarField UA = U.boundaryField()[patchi].patchInternalField() & mesh.Sf().boundaryField()[patchi]; //volScalarField UA = U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi]; scalar total_mass = gSum(UA.boundaryField()[patchi] * def_rho); scalar sumField = 0; if (total_mass > 0) { sumField = gSum ( UA.boundaryField()[patchi] * def_rho * field.boundaryField()[patchi] ) / total_mass; } Info<< " Average of " << fieldName << " over patch " << patchName << '[' << patchi << ']' << " = " << sumField << endl; Info<< " mesh.Sf().boundaryField()[patchi] " << mesh.Sf().boundaryField()[patchi] << "\n" << " U.boundaryField()[patchi] " << U.boundaryField()[patchi] << endl; } else { FatalError << "Only possible to average volScalarFields " << nl << exit(FatalError); } } else { Info<< " No field " << fieldName << endl; } Info<< endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // Any ideas? |
|
July 31, 2013, 17:51 |
|
#2 |
Member
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 16 |
The code is now running. It was a silly mistake.
This is the code: Code:
Application massAverage Description Calculates the average of the specified field over the specified density. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: int main(int argc, char *argv[]) { timeSelector::addOptions(); argList::validArgs.append("fieldName"); argList::validArgs.append("patchName"); # include "setRootCase.H" # include "createTime.H" instantList timeDirs = timeSelector::select0(runTime, args); # include "createMesh.H" word fieldName(args.additionalArgs()[0]); word patchName(args.additionalArgs()[1]); //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),args.additionalArgs()[2]); //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),scalar(1.226)); volScalarField def_rho ( IOobject ( "def_rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedScalar ( "def_rho", dimensionSet(1,-3,0,0,0,0,0), scalar(1.226) ) ); forAll(timeDirs, timeI) { runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; IOobject fieldHeader ( fieldName, runTime.timeName(), mesh, IOobject::MUST_READ ); IOobject U_header ( "U", runTime.timeName(), mesh, IOobject::MUST_READ ); // Check field exists if (fieldHeader.headerOk()) { mesh.readUpdate(); label patchi = mesh.boundaryMesh().findPatchID(patchName); if (patchi < 0) { FatalError << "Unable to find patch " << patchName << nl << exit(FatalError); } if (fieldHeader.headerClassName() == "volScalarField") { Info<< " Reading volScalarField " << fieldName << endl; volScalarField field(fieldHeader, mesh); volVectorField U(U_header, mesh); volScalarField UA ( IOobject ( "UA", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedScalar ( "UA", dimensionSet(0,2,-1,0,0,0,0), 0 ) ); UA.boundaryField()[patchi] = (U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi]); scalar total_mass = gSum(UA.boundaryField()[patchi] * def_rho.boundaryField()[patchi]); scalar sumField = 0; if (total_mass > 0) { sumField = gSum ( UA.boundaryField()[patchi] * def_rho * field.boundaryField()[patchi] ) / total_mass; } Info<< " Average of " << fieldName << " over patch " << patchName << '[' << patchi << ']' << " = " << sumField << endl; } else { FatalError << "Only possible to average volScalarFields " << nl << exit(FatalError); } } else { Info<< " No field " << fieldName << endl; } Info<< endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // Best regards, Pat |
|
August 1, 2013, 06:46 |
|
#3 |
Member
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 16 |
Since sum(scalar*rho*U*A)/sum(rho*U*A) = sum(scalar*U*A)/sum(U*A) for rho = const
and the inner product of U and A is saved in phi, I have changed my code: Code:
int main(int argc, char *argv[]) { timeSelector::addOptions(); argList::validArgs.append("fieldName"); argList::validArgs.append("patchName"); # include "setRootCase.H" # include "createTime.H" instantList timeDirs = timeSelector::select0(runTime, args); # include "createMesh.H" word fieldName(args.additionalArgs()[0]); word patchName(args.additionalArgs()[1]); forAll(timeDirs, timeI) { if(timeI == 0) continue; runTime.setTime(timeDirs[timeI], timeI); Info<< "Time = " << runTime.timeName() << endl; IOobject fieldHeader ( fieldName, runTime.timeName(), mesh, IOobject::MUST_READ ); IOobject phi_header ( "phi", runTime.timeName(), mesh, IOobject::MUST_READ ); // Check field exists if (fieldHeader.headerOk()) { mesh.readUpdate(); label patchi = mesh.boundaryMesh().findPatchID(patchName); if (patchi < 0) { FatalError << "Unable to find patch " << patchName << nl << exit(FatalError); } if (fieldHeader.headerClassName() == "volScalarField") { Info<< " Reading volScalarField " << fieldName << endl; volScalarField field(fieldHeader, mesh); surfaceScalarField phi(phi_header, mesh); scalar total_mass = gSum(phi.boundaryField()[patchi]); scalar sumField = 0; if (total_mass > 0) { sumField = gSum ( phi.boundaryField()[patchi] * field.boundaryField()[patchi] ) / total_mass; } Info<< " Average of " << fieldName << " over patch " << patchName << '[' << patchi << ']' << " = " << sumField << endl; } else { FatalError << "Only possible to average volScalarFields " << nl << exit(FatalError); } } else { Info<< " No field " << fieldName << endl; } Info<< endl; } Info<< "End\n" << endl; return 0; } You can now use it for openFoam1.6 and openFoam1.6-ext. Last edited by Pat84; August 5, 2013 at 10:36. |
|
April 16, 2014, 04:20 |
|
#4 |
Senior Member
Roman Thiele
Join Date: Aug 2009
Location: Eindhoven, NL
Posts: 374
Rep Power: 20 |
Hej,
I think the tool is ready needed. However, I was wondering if you can update it to be compiled in OF 2.3? I just checked to compile it. I get only one error message when compiling and I commented out one line in order to get it to compile. Line 53, what is it for? Code:
// if(timeDirs[timeI] == 0) continue; It only "knows" the standard patches, and does not know about the turbulence patch types. I think you can change that by actually including them in the compiling process. For example in my case it does not recognize the patch type compressible::turbulentHeatFluxTemperature.
__________________
~roman |
|
July 22, 2014, 15:13 |
|
#5 |
New Member
Jens
Join Date: Oct 2013
Posts: 2
Rep Power: 0 |
Hi,
good Tool. I´ve searched exactly for this. I am very lucky to find it In my opinion Line 53 should only skip the Time "Zero" and jump directly to the next time step. Because in Time "Zero" there is no "phi" file and nobody need the information about this time step. |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
y+ and u+ values with low-Re RANS turbulence models: utility + testcase | florian_krause | OpenFOAM | 114 | August 23, 2023 05:37 |
Area weighted or mass weighted average | SAM | Main CFD Forum | 31 | April 30, 2018 04:12 |
Diffference between mass weighted average and Area | giogio | FLUENT | 9 | March 6, 2018 08:24 |
mass weighted average temperature along a curve | preetam69 | Visualization & Post-Processing | 1 | January 14, 2014 11:33 |
mass weighted average temperature along a curve | preetam69 | FLUENT | 1 | July 5, 2013 02:48 |