From 5afc0db5a82c384821f3e51fbe997ef052490404 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 2 Feb 2009 12:23:41 +0000 Subject: [PATCH 1/3] Solaris port --- etc/bashrc | 11 +++++++++++ etc/cshrc | 13 ++++++++++++- etc/settings.csh | 9 +++++++++ etc/settings.sh | 10 ++++++++++ src/OSspecific/Unix/Make/files | 6 +++++- wmake/rules/SunOS64Gcc/X | 3 +++ wmake/rules/SunOS64Gcc/c | 16 ++++++++++++++++ wmake/rules/SunOS64Gcc/c++ | 21 +++++++++++++++++++++ wmake/rules/SunOS64Gcc/c++Debug | 2 ++ wmake/rules/SunOS64Gcc/c++Opt | 2 ++ wmake/rules/SunOS64Gcc/c++Prof | 2 ++ wmake/rules/SunOS64Gcc/cDebug | 2 ++ wmake/rules/SunOS64Gcc/cOpt | 2 ++ wmake/rules/SunOS64Gcc/cProf | 2 ++ wmake/rules/SunOS64Gcc/general | 11 +++++++++++ wmake/rules/SunOS64Gcc/mplib | 3 +++ wmake/rules/SunOS64Gcc/mplibFJMPI | 3 +++ wmake/rules/SunOS64Gcc/mplibOPENMPI | 3 +++ 18 files changed, 119 insertions(+), 2 deletions(-) create mode 100644 wmake/rules/SunOS64Gcc/X create mode 100644 wmake/rules/SunOS64Gcc/c create mode 100644 wmake/rules/SunOS64Gcc/c++ create mode 100644 wmake/rules/SunOS64Gcc/c++Debug create mode 100644 wmake/rules/SunOS64Gcc/c++Opt create mode 100644 wmake/rules/SunOS64Gcc/c++Prof create mode 100644 wmake/rules/SunOS64Gcc/cDebug create mode 100644 wmake/rules/SunOS64Gcc/cOpt create mode 100644 wmake/rules/SunOS64Gcc/cProf create mode 100644 wmake/rules/SunOS64Gcc/general create mode 100644 wmake/rules/SunOS64Gcc/mplib create mode 100644 wmake/rules/SunOS64Gcc/mplibFJMPI create mode 100644 wmake/rules/SunOS64Gcc/mplibOPENMPI diff --git a/etc/bashrc b/etc/bashrc index 7c9f9ecad3..e9f9c6f9cf 100644 --- a/etc/bashrc +++ b/etc/bashrc @@ -172,6 +172,17 @@ Linux) esac ;; +SunOS) + WM_ARCH=SunOS64 + export WM_COMPILER_LIB_ARCH=64 + export WM_CC='gcc' + export WM_CXX='g++' + export WM_CFLAGS='-mabi=64 -fPIC' + export WM_CXXFLAGS='-mabi=64 -fPIC' + export WM_LDFLAGS='-mabi=64 -G0' + export WM_MPLIB=FJMPI + ;; + *) # an unsupported operating system cat < Date: Mon, 2 Feb 2009 12:24:32 +0000 Subject: [PATCH 2/3] cutting plane --- .../postProcessing/sampling/sample/sampleDict | 19 + src/sampling/Make/files | 1 + .../sampledCuttingPlane/sampledCuttingPlane.C | 420 ++++++++++++++++++ .../sampledCuttingPlane/sampledCuttingPlane.H | 258 +++++++++++ .../sampledCuttingPlaneTemplates.C | 97 ++++ 5 files changed, 795 insertions(+) create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict index 1616bbb662..b605f29e57 100644 --- a/applications/utilities/postProcessing/sampling/sample/sampleDict +++ b/applications/utilities/postProcessing/sampling/sample/sampleDict @@ -177,6 +177,25 @@ surfaces interpolate false; regularise false; // do not simplify } + + triangleCut + { + // Cutingplaneusing iso surface + type cuttingPlane; + planeType pointAndNormal; + pointAndNormalDict + { + basePoint (0.4 0 0.4); + normalVector (1 0.2 0.2); + } + interpolate true; + + //zone ABC; // Optional: zone only + //exposedPatchName fixedWalls; // Optional: zone only + + // regularise false; // Optional: do not simplify + } + ); diff --git a/src/sampling/Make/files b/src/sampling/Make/files index 43bc910a3c..004f81d435 100644 --- a/src/sampling/Make/files +++ b/src/sampling/Make/files @@ -29,6 +29,7 @@ sampledSurface/isoSurface/sampledIsoSurface.C sampledSurface/isoSurface/isoSurfaceCell.C sampledSurface/isoSurface/sampledIsoSurfaceCell.C sampledSurface/distanceSurface/distanceSurface.C +sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C sampledSurface/sampledSurface/sampledSurface.C sampledSurface/sampledSurfaces/sampledSurfaces.C sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C new file mode 100644 index 0000000000..85dbeca595 --- /dev/null +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C @@ -0,0 +1,420 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 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 + +\*---------------------------------------------------------------------------*/ + +#include "sampledCuttingPlane.H" +#include "dictionary.H" +#include "volFields.H" +#include "volPointInterpolation.H" +#include "addToRunTimeSelectionTable.H" +#include "fvMesh.H" +#include "isoSurface.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(sampledCuttingPlane, 0); + addNamedToRunTimeSelectionTable(sampledSurface, sampledCuttingPlane, word, cuttingPlane); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sampledCuttingPlane::createGeometry() +{ + if (debug) + { + Pout<< "sampledCuttingPlane::createGeometry :updating geometry." + << endl; + } + + // Clear any stored topologies + facesPtr_.clear(); + isoSurfPtr_.ptr(); + pointDistance_.clear(); + cellDistancePtr_.clear(); + + + // Get any subMesh + if (zoneID_.index() != -1 && !subMeshPtr_.valid()) + { + const polyBoundaryMesh& patches = mesh().boundaryMesh(); + + // Patch to put exposed internal faces into + label exposedPatchI = patches.findPatchID(exposedPatchName_); + + if (debug) + { + Info<< "Allocating subset of size " + << mesh().cellZones()[zoneID_.index()].size() + << " with exposed faces into patch " + << patches[exposedPatchI].name() << endl; + } + + subMeshPtr_.reset + ( + new fvMeshSubset(static_cast(mesh())) + ); + subMeshPtr_().setLargeCellSubset + ( + labelHashSet(mesh().cellZones()[zoneID_.index()]), + exposedPatchI + ); + } + + + // Select either the submesh or the underlying mesh + const fvMesh& fvm = + ( + subMeshPtr_.valid() + ? subMeshPtr_().subMesh() + : static_cast(mesh()) + ); + + + // Distance to cell centres + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + cellDistancePtr_.reset + ( + new volScalarField + ( + IOobject + ( + "cellDistance", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + fvm, + dimensionedScalar("zero", dimLength, 0) + ) + ); + volScalarField& cellDistance = cellDistancePtr_(); + + // Internal field + { + const pointField& cc = fvm.C(); + scalarField& fld = cellDistance.internalField(); + + forAll(cc, i) + { + // Signed distance + fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal(); + } + } + + // Patch fields + { + forAll(fvm.C().boundaryField(), patchI) + { + const pointField& cc = fvm.C().boundaryField()[patchI]; + fvPatchScalarField& fld = cellDistance.boundaryField()[patchI]; + + forAll(fld, i) + { + fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal(); + } + } + } + + + // On processor patches the mesh.C() will already be the cell centre + // on the opposite side so no need to swap cellDistance. + + + // Distance to points + pointDistance_.setSize(fvm.nPoints()); + { + const pointField& pts = fvm.points(); + + forAll(pointDistance_, i) + { + pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal(); + } + } + + + if (debug) + { + Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl; + cellDistance.write(); + pointScalarField pDist + ( + IOobject + ( + "pointDistance", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + pointMesh::New(fvm), + dimensionedScalar("zero", dimLength, 0) + ); + pDist.internalField() = pointDistance_; + + Pout<< "Writing point distance:" << pDist.objectPath() << endl; + pDist.write(); + } + + + //- Direct from cell field and point field. + isoSurfPtr_.reset + ( + new isoSurface + ( + cellDistance, + pointDistance_, + 0.0, + regularise_ + ) + ); + + if (debug) + { + print(Pout); + Pout<< endl; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sampledCuttingPlane::sampledCuttingPlane +( + const word& name, + const polyMesh& mesh, + const dictionary& dict +) +: + sampledSurface(name, mesh, dict), + plane_(dict), + mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)), + regularise_(dict.lookupOrDefault("regularise", true)), + zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()), + exposedPatchName_(word::null), + needsUpdate_(true), + subMeshPtr_(NULL), + cellDistancePtr_(NULL), + isoSurfPtr_(NULL), + facesPtr_(NULL) +{ + if (zoneID_.index() != -1) + { + dict.lookup("exposedPatchName") >> exposedPatchName_; + + if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1) + { + FatalErrorIn + ( + "sampledCuttingPlane::sampledCuttingPlane" + "(const word&, const polyMesh&, const dictionary&)" + ) << "Cannot find patch " << exposedPatchName_ + << " in which to put exposed faces." << endl + << "Valid patches are " << mesh.boundaryMesh().names() + << exit(FatalError); + } + + if (debug && zoneID_.index() != -1) + { + Info<< "Restricting to cellZone " << zoneID_.name() + << " with exposed internal faces into patch " + << exposedPatchName_ << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sampledCuttingPlane::~sampledCuttingPlane() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::sampledCuttingPlane::needsUpdate() const +{ + return needsUpdate_; +} + + +bool Foam::sampledCuttingPlane::expire() +{ + if (debug) + { + Pout<< "sampledCuttingPlane::expire :" + << " have-facesPtr_:" << facesPtr_.valid() + << " needsUpdate_:" << needsUpdate_ << endl; + } + + // Clear any stored topologies + facesPtr_.clear(); + + // already marked as expired + if (needsUpdate_) + { + return false; + } + + needsUpdate_ = true; + return true; +} + + +bool Foam::sampledCuttingPlane::update() +{ + if (debug) + { + Pout<< "sampledCuttingPlane::update :" + << " have-facesPtr_:" << facesPtr_.valid() + << " needsUpdate_:" << needsUpdate_ << endl; + } + + if (!needsUpdate_) + { + return false; + } + + createGeometry(); + + needsUpdate_ = false; + return true; +} + + +Foam::tmp +Foam::sampledCuttingPlane::sample +( + const volScalarField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp +Foam::sampledCuttingPlane::sample +( + const volVectorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp +Foam::sampledCuttingPlane::sample +( + const volSphericalTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp +Foam::sampledCuttingPlane::sample +( + const volSymmTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp +Foam::sampledCuttingPlane::sample +( + const volTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp +Foam::sampledCuttingPlane::interpolate +( + const interpolation& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp +Foam::sampledCuttingPlane::interpolate +( + const interpolation& interpolator +) const +{ + return interpolateField(interpolator); +} + +Foam::tmp +Foam::sampledCuttingPlane::interpolate +( + const interpolation& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp +Foam::sampledCuttingPlane::interpolate +( + const interpolation& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp +Foam::sampledCuttingPlane::interpolate +( + const interpolation& interpolator +) const +{ + return interpolateField(interpolator); +} + + +void Foam::sampledCuttingPlane::print(Ostream& os) const +{ + os << "sampledCuttingPlane: " << name() << " :" + << " plane:" << plane_ + << " faces:" << faces().size() + << " points:" << points().size(); +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H new file mode 100644 index 0000000000..fcf14db8bd --- /dev/null +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H @@ -0,0 +1,258 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 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 + +Class + Foam::sampledCuttingPlane + +Description + A sampledSurface defined by a plane + +SourceFiles + sampledCuttingPlane.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sampledCuttingPlane_H +#define sampledCuttingPlane_H + +#include "sampledSurface.H" +#include "isoSurface.H" +#include "plane.H" +#include "ZoneIDs.H" +#include "fvMeshSubset.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class sampledCuttingPlane Declaration +\*---------------------------------------------------------------------------*/ + +class sampledCuttingPlane +: + public sampledSurface +{ + // Private data + + //- Plane + const plane plane_; + + //- Merge tolerance + const scalar mergeTol_; + + //- Whether to coarsen + const Switch regularise_; + + //- zone name/index (if restricted to zones) + mutable cellZoneID zoneID_; + + //- for zones: patch to put exposed faces into + mutable word exposedPatchName_; + + //- Track if the surface needs an update + mutable bool needsUpdate_; + + + //- Optional subsetted mesh + autoPtr subMeshPtr_; + + //- Distance to cell centres + autoPtr cellDistancePtr_; + + //- Distance to points + scalarField pointDistance_; + + //- Constructed iso surface + autoPtr isoSurfPtr_; + + //- triangles converted to faceList + mutable autoPtr facesPtr_; + + + // Private Member Functions + + //- Create iso surface + void createGeometry(); + + //- sample field on faces + template + tmp > sampleField + ( + const GeometricField& vField + ) const; + + + template + tmp > + interpolateField(const interpolation&) const; + + +public: + + //- Runtime type information + TypeName("sampledCuttingPlane"); + + + // Constructors + + //- Construct from dictionary + sampledCuttingPlane + ( + const word& name, + const polyMesh& mesh, + const dictionary& dict + ); + + + // Destructor + + virtual ~sampledCuttingPlane(); + + + // Member Functions + + //- Does the surface need an update? + virtual bool needsUpdate() const; + + //- Mark the surface as needing an update. + // May also free up unneeded data. + // Return false if surface was already marked as expired. + virtual bool expire(); + + //- Update the surface as required. + // Do nothing (and return false) if no update was needed + virtual bool update(); + + //- Points of surface + virtual const pointField& points() const + { + return surface().points(); + } + + //- Faces of surface + virtual const faceList& faces() const + { + if (facesPtr_.empty()) + { + const triSurface& s = surface(); + + facesPtr_.reset(new faceList(s.size())); + + forAll(s, i) + { + facesPtr_()[i] = s[i].triFaceFace(); + } + } + return facesPtr_; + } + + + const isoSurface& surface() const + { + return isoSurfPtr_(); + } + + //- sample field on surface + virtual tmp sample + ( + const volScalarField& + ) const; + + //- sample field on surface + virtual tmp sample + ( + const volVectorField& + ) const; + + //- sample field on surface + virtual tmp sample + ( + const volSphericalTensorField& + ) const; + + //- sample field on surface + virtual tmp sample + ( + const volSymmTensorField& + ) const; + + //- sample field on surface + virtual tmp sample + ( + const volTensorField& + ) const; + + + //- interpolate field on surface + virtual tmp interpolate + ( + const interpolation& + ) const; + + //- interpolate field on surface + virtual tmp interpolate + ( + const interpolation& + ) const; + + //- interpolate field on surface + virtual tmp interpolate + ( + const interpolation& + ) const; + + //- interpolate field on surface + virtual tmp interpolate + ( + const interpolation& + ) const; + + //- interpolate field on surface + virtual tmp interpolate + ( + const interpolation& + ) const; + + //- Write + virtual void print(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "sampledCuttingPlaneTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C new file mode 100644 index 0000000000..f83aeca4e5 --- /dev/null +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C @@ -0,0 +1,97 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 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 + +\*---------------------------------------------------------------------------*/ + +#include "volPointInterpolation.H" +#include "sampledCuttingPlane.H" +#include "isoSurface.H" +#include "volFieldsFwd.H" +#include "pointFields.H" +#include "volPointInterpolation.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +Foam::tmp > +Foam::sampledCuttingPlane::sampleField +( + const GeometricField& vField +) const +{ + return tmp >(new Field(vField, surface().meshCells())); +} + + +template +Foam::tmp > +Foam::sampledCuttingPlane::interpolateField +( + const interpolation& interpolator +) const +{ + // Get fields to sample. Assume volPointInterpolation! + const GeometricField& volFld = + interpolator.psi(); + + if (subMeshPtr_.valid()) + { + tmp > tvolSubFld = + subMeshPtr_().interpolate(volFld); + + const GeometricField& volSubFld = + tvolSubFld(); + + tmp > tpointSubFld = + volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld); + + // Sample. + return surface().interpolate + ( + cellDistancePtr_(), + pointDistance_, + volSubFld, + tpointSubFld() + ); + } + else + { + tmp > tpointFld + ( + volPointInterpolation::New(volFld.mesh()).interpolate(volFld) + ); + + // Sample. + return surface().interpolate + ( + cellDistancePtr_(), + pointDistance_, + volFld, + tpointFld() + ); + } +} + + +// ************************************************************************* // From 338b72b1ebbe9b83f81c3d021f9d71c0205c5600 Mon Sep 17 00:00:00 2001 From: henry Date: Tue, 3 Feb 2009 16:51:07 +0000 Subject: [PATCH 3/3] Improved the flux-velocity correspondence for cases where hydrostatic balance is important e.g. in atmospheric flows. --- .../solvers/heatTransfer/buoyantFoam/UEqn.H | 16 ++- .../heatTransfer/buoyantFoam/buoyantFoam.C | 28 ++--- .../heatTransfer/buoyantFoam/createFields.H | 1 + .../solvers/heatTransfer/buoyantFoam/pEqn.H | 103 ++++++++++-------- .../heatTransfer/buoyantSimpleFoam/UEqn.H | 10 +- .../buoyantSimpleFoam/createFields.H | 1 + .../heatTransfer/buoyantSimpleFoam/pEqn.H | 98 +++++++++-------- .../buoyantSimpleRadiationFoam/Make/options | 1 + .../buoyantSimpleRadiationFoam/UEqn.H | 18 --- .../buoyantSimpleRadiationFoam/createFields.H | 1 + .../buoyantSimpleRadiationFoam/pEqn.H | 54 --------- 11 files changed, 151 insertions(+), 180 deletions(-) delete mode 100644 applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H delete mode 100644 applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H diff --git a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H index 988268c4ef..24cfbf4f68 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H @@ -9,4 +9,18 @@ UEqn.relax(); - solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh); + if (momentumPredictor) + { + solve + ( + UEqn + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C index cafb677900..8b3cca8d21 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C +++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C @@ -41,13 +41,12 @@ Description int main(int argc, char *argv[]) { - -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "readEnvironmentalProperties.H" -# include "createFields.H" -# include "initContinuityErrs.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readEnvironmentalProperties.H" + #include "createFields.H" + #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,23 +54,24 @@ int main(int argc, char *argv[]) while (runTime.run()) { -# include "readPISOControls.H" -# include "compressibleCourantNo.H" -//# include "setDeltaT.H" + #include "readTimeControls.H" + #include "readPISOControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; -# include "rhoEqn.H" + #include "rhoEqn.H" -# include "UEqn.H" + #include "UEqn.H" // --- PISO loop for (int corr=0; corrcorrect(); diff --git a/applications/solvers/heatTransfer/buoyantFoam/createFields.H b/applications/solvers/heatTransfer/buoyantFoam/createFields.H index 3d604d4daf..aa32325c34 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H @@ -57,6 +57,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H index 36378f7f6f..9feededcd3 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H @@ -1,60 +1,67 @@ -bool closedVolume = pd.needReference(); - -rho = thermo->rho(); - -volScalarField rUA = 1.0/UEqn.A(); -U = rUA*UEqn.H(); - -phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) - ) - - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn + bool closedVolume = pd.needReference(); + + rho = thermo->rho(); + + volScalarField rUA = 1.0/UEqn.A(); + surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); + + U = rUA*UEqn.H(); + + surfaceScalarField phiU ( - fvm::ddt(psi, pd) - + fvc::ddt(psi)*pRef - + fvc::ddt(psi, rho)*gh - + fvc::div(phi) - - fvm::laplacian(rho*rUA, pd) + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ) ); - if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - pdEqn.solve(mesh.solver(pd.name() + "Final")); - } - else - { - pdEqn.solve(mesh.solver(pd.name())); + fvScalarMatrix pdEqn + ( + fvm::ddt(psi, pd) + + fvc::ddt(psi)*pRef + + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + - fvm::laplacian(rhorUAf, pd) + ); + + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pdEqn.solve(mesh.solver(pd.name() + "Final")); + } + else + { + pdEqn.solve(mesh.solver(pd.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi += pdEqn.flux(); + } } - if (nonOrth == nNonOrthCorr) + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + + p == pd + rho*gh + pRef; + dpdt = fvc::ddt(p); + + #include "rhoEqn.H" + #include "compressibleContinuityErrs.H" + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi += pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); + rho = thermo->rho(); } -} -p == pd + rho*gh + pRef; -dpdt = fvc::ddt(p); - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); - rho = thermo->rho(); } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H index 5bc4c4738c..71a57cd2a9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H @@ -11,7 +11,15 @@ eqnResidual = solve ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index 184242be81..d26da5cb32 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -53,6 +53,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 8c3621205b..3c257d249d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,55 +1,65 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); - -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); + volScalarField rUA = 1.0/UEqn().A(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) + U = rUA*UEqn().H(); + UEqn.clear(); + + phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); + bool closedVolume = adjustPhi(phi, U, p); + surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + phi -= buoyancyPhi; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); + fvScalarMatrix pdEqn + ( + fvm::laplacian(rhorUAf, pd) == fvc::div(phi) + ); + + pdEqn.setReference(pdRefCell, pdRefValue); + + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pdEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pdEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + // Calculate the conservative fluxes + phi -= pdEqn.flux(); + + // Explicitly relax pressure for momentum corrector + pd.relax(); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf); + U.correctBoundaryConditions(); + } } - if (nonOrth == nNonOrthCorr) + #include "continuityErrs.H" + + p == pd + rho*gh + pRef; + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi -= pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); } -} -#include "continuityErrs.H" + rho = thermo->rho(); + rho.relax(); + Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); } - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options index 22d41fdb80..89d7787af6 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I../buoyantSimpleFoam \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H deleted file mode 100644 index 5bc4c4738c..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H +++ /dev/null @@ -1,18 +0,0 @@ - // Solve the Momentum equation - - tmp UEqn - ( - fvm::div(phi, U) - - fvm::Sp(fvc::div(phi), U) - + turbulence->divDevRhoReff(U) - ); - - UEqn().relax(); - - eqnResidual = solve - ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh - ).initialResidual(); - - maxResidual = max(eqnResidual, maxResidual); - diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H index 76e3805bca..62c06ec38d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H @@ -54,6 +54,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H deleted file mode 100644 index 2a713bac62..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H +++ /dev/null @@ -1,54 +0,0 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) -{ - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); - - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) - { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); - } - - if (nonOrth == nNonOrthCorr) - { - phi -= pdEqn.flux(); - } -} - -#include "continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); - pd == p - (rho*gh + pRef); -} - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;