From 37de1d4badb02b6fbf47cd9c4f524810afcb243c Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 12 Apr 2012 17:40:05 +0100 Subject: [PATCH 01/82] ENH: timeVaryingMappedFixedValuePoint: pointField version of timeVaryingMapped --- .../timeVaryingMappedFixedValueFvPatchField.C | 441 +++------------- .../timeVaryingMappedFixedValueFvPatchField.H | 47 +- src/fvMotionSolver/Make/files | 20 +- ...meVaryingMappedFixedValuePointPatchField.C | 482 ++++++++++++++++++ ...meVaryingMappedFixedValuePointPatchField.H | 195 +++++++ ...eVaryingMappedFixedValuePointPatchFields.C | 43 ++ ...eVaryingMappedFixedValuePointPatchFields.H | 57 +++ src/meshTools/Make/files | 1 + .../pointToPointPlanarInterpolation.C | 322 ++++++++++++ .../pointToPointPlanarInterpolation.H | 164 ++++++ ...pointToPointPlanarInterpolationTemplates.C | 71 +++ .../triSurfaceTools/triSurfaceTools.H | 4 +- 12 files changed, 1438 insertions(+), 409 deletions(-) create mode 100644 src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C create mode 100644 src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H create mode 100644 src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C create mode 100644 src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H create mode 100644 src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C create mode 100644 src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H create mode 100644 src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C index 61fb24ba77..1831bd903d 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,12 +25,7 @@ License #include "timeVaryingMappedFixedValueFvPatchField.H" #include "Time.H" -#include "triSurfaceTools.H" -#include "triSurface.H" -#include "vector2D.H" -#include "OFstream.H" #include "AverageIOField.H" -#include "Random.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -51,9 +46,6 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(iF.name()), setAverage_(false), perturb_(0), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -78,9 +70,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -104,9 +94,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(iF.name()), setAverage_(readBool(dict.lookup("setAverage"))), perturb_(dict.lookupOrDefault("perturb", 1E-5)), - referenceCS_(NULL), - nearestVertex_(0), - nearestVertexWeight_(0), + mapperPtr_(NULL), sampleTimes_(0), startSampleTime_(-1), startSampledValues_(0), @@ -139,9 +127,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(ptf.referenceCS_), - nearestVertex_(ptf.nearestVertex_), - nearestVertexWeight_(ptf.nearestVertexWeight_), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), startSampledValues_(ptf.startSampledValues_), @@ -165,9 +151,7 @@ timeVaryingMappedFixedValueFvPatchField fieldTableName_(ptf.fieldTableName_), setAverage_(ptf.setAverage_), perturb_(ptf.perturb_), - referenceCS_(ptf.referenceCS_), - nearestVertex_(ptf.nearestVertex_), - nearestVertexWeight_(ptf.nearestVertexWeight_), + mapperPtr_(ptf.mapperPtr_), sampleTimes_(ptf.sampleTimes_), startSampleTime_(ptf.startSampleTime_), startSampledValues_(ptf.startSampledValues_), @@ -192,6 +176,10 @@ void timeVaryingMappedFixedValueFvPatchField::autoMap startSampledValues_.autoMap(m); endSampledValues_.autoMap(m); } + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; } @@ -209,295 +197,11 @@ void timeVaryingMappedFixedValueFvPatchField::rmap startSampledValues_.rmap(tiptf.startSampledValues_, addr); endSampledValues_.rmap(tiptf.endSampledValues_, addr); -} - -template -void timeVaryingMappedFixedValueFvPatchField::readSamplePoints() -{ - // Read the sample points - - pointIOField samplePoints - ( - IOobject - ( - "points", - this->db().time().constant(), - "boundaryData"/this->patch().name(), - this->db(), - IOobject::MUST_READ, - IOobject::AUTO_WRITE, - false - ) - ); - - const fileName samplePointsFile = samplePoints.filePath(); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Read " << samplePoints.size() << " sample points from " - << samplePointsFile << endl; - } - - // Determine coordinate system from samplePoints - - if (samplePoints.size() < 3) - { - FatalErrorIn - ( - "timeVaryingMappedFixedValueFvPatchField::readSamplePoints()" - ) << "Only " << samplePoints.size() << " points read from file " - << samplePoints.objectPath() << nl - << "Need at least three non-colinear samplePoints" - << " to be able to interpolate." - << "\n on patch " << this->patch().name() - << " of points " << samplePoints.name() - << " in file " << samplePoints.objectPath() - << exit(FatalError); - } - - const point& p0 = samplePoints[0]; - - // Find furthest away point - vector e1; - label index1 = -1; - scalar maxDist = -GREAT; - - for (label i = 1; i < samplePoints.size(); i++) - { - const vector d = samplePoints[i] - p0; - scalar magD = mag(d); - - if (magD > maxDist) - { - e1 = d/magD; - index1 = i; - maxDist = magD; - } - } - // Find point that is furthest away from line p0-p1 - const point& p1 = samplePoints[index1]; - - label index2 = -1; - maxDist = -GREAT; - for (label i = 1; i < samplePoints.size(); i++) - { - if (i != index1) - { - const point& p2 = samplePoints[i]; - vector e2(p2 - p0); - e2 -= (e2&e1)*e1; - scalar magE2 = mag(e2); - - if (magE2 > maxDist) - { - index2 = i; - maxDist = magE2; - } - } - } - if (index2 == -1) - { - FatalErrorIn - ( - "timeVaryingMappedFixedValueFvPatchField::readSamplePoints()" - ) << "Cannot find points that make valid normal." << nl - << "Have so far points " << p0 << " and " << p1 - << "Need at least three sample points which are not in a line." - << "\n on patch " << this->patch().name() - << " of points " << samplePoints.name() - << " in file " << samplePoints.objectPath() - << exit(FatalError); - } - - vector n = e1^(samplePoints[index2]-p0); - n /= mag(n); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Used points " << p0 << ' ' << samplePoints[index1] - << ' ' << samplePoints[index2] - << " to define coordinate system with normal " << n << endl; - } - - referenceCS_.reset - ( - new coordinateSystem - ( - "reference", - p0, // origin - n, // normal - e1 // 0-axis - ) - ); - - tmp tlocalVertices - ( - referenceCS().localPosition(samplePoints) - ); - vectorField& localVertices = tlocalVertices(); - - const boundBox bb(localVertices, true); - const point bbMid(bb.midpoint()); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField :" - << " Perturbing points with " << perturb_ - << " fraction of a random position inside " << bb - << " to break any ties on regular meshes." - << nl << endl; - } - - Random rndGen(123456); - forAll(localVertices, i) - { - localVertices[i] += - perturb_ - *(rndGen.position(bb.min(), bb.max())-bbMid); - } - - // Determine triangulation - List localVertices2D(localVertices.size()); - forAll(localVertices, i) - { - localVertices2D[i][0] = localVertices[i][0]; - localVertices2D[i][1] = localVertices[i][1]; - } - - triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); - - tmp tlocalFaceCentres - ( - referenceCS().localPosition - ( - this->patch().patch().faceCentres() - ) - ); - const pointField& localFaceCentres = tlocalFaceCentres(); - - if (debug) - { - Pout<< "readSamplePoints :" - <<" Dumping triangulated surface to triangulation.stl" << endl; - s.write(this->db().time().path()/"triangulation.stl"); - - OFstream str(this->db().time().path()/"localFaceCentres.obj"); - Pout<< "readSamplePoints :" - << " Dumping face centres to " << str.name() << endl; - - forAll(localFaceCentres, i) - { - const point& p = localFaceCentres[i]; - str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; - } - } - - // Determine interpolation onto face centres. - triSurfaceTools::calcInterpolationWeights - ( - s, - localFaceCentres, // points to interpolate to - nearestVertex_, - nearestVertexWeight_ - ); - - - - // Read the times for which data is available - - const fileName samplePointsDir = samplePointsFile.path(); - - sampleTimes_ = Time::findTimes(samplePointsDir); - - if (debug) - { - Info<< "timeVaryingMappedFixedValueFvPatchField : In directory " - << samplePointsDir << " found times " << timeNames(sampleTimes_) - << endl; - } -} - - -template -wordList timeVaryingMappedFixedValueFvPatchField::timeNames -( - const instantList& times -) -{ - wordList names(times.size()); - - forAll(times, i) - { - names[i] = times[i].name(); - } - return names; -} - - -template -void timeVaryingMappedFixedValueFvPatchField::findTime -( - const fileName& instance, - const fileName& local, - const scalar timeVal, - label& lo, - label& hi -) const -{ - lo = startSampleTime_; - hi = -1; - - for (label i = startSampleTime_+1; i < sampleTimes_.size(); i++) - { - if (sampleTimes_[i].value() > timeVal) - { - break; - } - else - { - lo = i; - } - } - - if (lo == -1) - { - FatalErrorIn("findTime") - << "Cannot find starting sampling values for current time " - << timeVal << nl - << "Have sampling values for times " - << timeNames(sampleTimes_) << nl - << "In directory " - << this->db().time().constant()/"boundaryData"/this->patch().name() - << "\n on patch " << this->patch().name() - << " of field " << fieldTableName_ - << exit(FatalError); - } - - if (lo < sampleTimes_.size()-1) - { - hi = lo+1; - } - - - if (debug) - { - if (hi == -1) - { - Pout<< "findTime : Found time " << timeVal << " after" - << " index:" << lo << " time:" << sampleTimes_[lo].value() - << endl; - } - else - { - Pout<< "findTime : Found time " << timeVal << " inbetween" - << " index:" << lo << " time:" << sampleTimes_[lo].value() - << " and index:" << hi << " time:" << sampleTimes_[hi].value() - << endl; - } - } + // Clear interpolator + mapperPtr_.clear(); + startSampleTime_ = -1; + endSampleTime_ = -1; } @@ -507,22 +211,84 @@ void timeVaryingMappedFixedValueFvPatchField::checkTable() // Initialise if (startSampleTime_ == -1 && endSampleTime_ == -1) { - readSamplePoints(); + pointIOField samplePoints + ( + IOobject + ( + "points", + this->db().time().constant(), + "boundaryData"/this->patch().name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + const fileName samplePointsFile = samplePoints.filePath(); + + if (debug) + { + Info<< "timeVaryingMappedFixedValueFvPatchField :" + << " Read " << samplePoints.size() << " sample points from " + << samplePointsFile << endl; + } + + // Allocate the interpolator + mapperPtr_.reset + ( + new pointToPointPlanarInterpolation + ( + samplePoints, + this->patch().patch().faceCentres(), + perturb_ + ) + ); + + // Read the times for which data is available + const fileName samplePointsDir = samplePointsFile.path(); + sampleTimes_ = Time::findTimes(samplePointsDir); + + if (debug) + { + Info<< "timeVaryingMappedFixedValueFvPatchField : In directory " + << samplePointsDir << " found times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) + << endl; + } } + // Find current time in sampleTimes label lo = -1; label hi = -1; - findTime + bool foundTime = mapperPtr_().findTime ( - this->db().time().constant(), - "boundaryData"/this->patch().name(), + sampleTimes_, + startSampleTime_, this->db().time().value(), lo, hi ); + if (!foundTime) + { + FatalErrorIn + ( + "timeVaryingMappedFixedValueFvPatchField::checkTable" + ) << "Cannot find starting sampling values for current time " + << this->db().time().value() << nl + << "Have sampling values for times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl + << "In directory " + << this->db().time().constant()/"boundaryData"/this->patch().name() + << "\n on patch " << this->patch().name() + << " of field " << fieldTableName_ + << exit(FatalError); + } + + // Update sampled data fields. if (lo != startSampleTime_) @@ -573,7 +339,7 @@ void timeVaryingMappedFixedValueFvPatchField::checkTable() ); startAverage_ = vals.average(); - startSampledValues_ = interpolate(vals); + startSampledValues_ = mapperPtr_().interpolate(vals); } } @@ -617,61 +383,22 @@ void timeVaryingMappedFixedValueFvPatchField::checkTable() ) ); endAverage_ = vals.average(); - endSampledValues_ = interpolate(vals); + endSampledValues_ = mapperPtr_().interpolate(vals); } } } -template -tmp > timeVaryingMappedFixedValueFvPatchField::interpolate -( - const Field& sourceFld -) const -{ - tmp > tfld(new Field(nearestVertex_.size())); - Field& fld = tfld(); - - forAll(fld, i) - { - const FixedList& verts = nearestVertex_[i]; - const FixedList& w = nearestVertexWeight_[i]; - - if (verts[2] == -1) - { - if (verts[1] == -1) - { - // Use vertex0 only - fld[i] = sourceFld[verts[0]]; - } - else - { - // Use vertex 0,1 - fld[i] = - w[0]*sourceFld[verts[0]] - + w[1]*sourceFld[verts[1]]; - } - } - else - { - fld[i] = - w[0]*sourceFld[verts[0]] - + w[1]*sourceFld[verts[1]] - + w[2]*sourceFld[verts[2]]; - } - } - return tfld; -} - - template void timeVaryingMappedFixedValueFvPatchField::updateCoeffs() { + if (this->updated()) { return; } + checkTable(); // Interpolate between the sampled data @@ -768,7 +495,7 @@ void timeVaryingMappedFixedValueFvPatchField::write(Ostream& os) const { fvPatchField::write(os); os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl; - os.writeKeyword("peturb") << perturb_ << token::END_STATEMENT << nl; + os.writeKeyword("perturb") << perturb_ << token::END_STATEMENT << nl; if (fieldTableName_ != this->dimensionedInternalField().name()) { diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H index e3ac112154..1494efec35 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,9 +66,9 @@ SourceFiles #define timeVaryingMappedFixedValueFvPatchField_H #include "fixedValueFvPatchFields.H" -#include "coordinateSystem.H" #include "FixedList.H" #include "instantList.H" +#include "pointToPointPlanarInterpolation.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -95,16 +95,8 @@ class timeVaryingMappedFixedValueFvPatchField //- Fraction of perturbation (fraction of bounding box) to add scalar perturb_; - //- Coordinate system - autoPtr referenceCS_; - - //- Current interpolation addressing to face centres of underlying - // patch - List > nearestVertex_; - - //- Current interpolation factors to face centres of underlying - // patch - List > nearestVertexWeight_; + //- 2D interpolation + autoPtr mapperPtr_; //- List of boundaryData time directories instantList sampleTimes_; @@ -127,31 +119,6 @@ class timeVaryingMappedFixedValueFvPatchField //- If setAverage: end average value Type endAverage_; - - // Private Member Functions - - //- Get names of times - static wordList timeNames(const instantList&); - - //- Find times around current time - void findTime - ( - const fileName& instance, - const fileName& local, - const scalar timeVal, - label& lo, - label& hi - ) const; - - - //- Read boundary points and determine interpolation weights to patch - // faceCentres - void readSamplePoints(); - - //- Do actual interpolation using current weights - tmp > interpolate(const Field&) const; - - public: //- Runtime type information @@ -224,12 +191,6 @@ public: // Access - //- Return the coordinateSystem - const coordinateSystem& referenceCS() const - { - return referenceCS_; - } - //- Return startSampledValues const Field startSampledValues() { diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files index bea2f10656..43dc7513e7 100644 --- a/src/fvMotionSolver/Make/files +++ b/src/fvMotionSolver/Make/files @@ -26,13 +26,19 @@ motionDiffusivity/manipulators/exponential/exponentialDiffusivity.C fvPatchFields/derived/cellMotion/cellMotionFvPatchFields.C fvPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementFvPatchFields.C -pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C -pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C +derivedPoint = pointPatchFields/derived + +$(derivedPoint)/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C +$(derivedPoint)/angularOscillatingVelocity/angularOscillatingVelocityPointPatchVectorField.C + +$(derivedPoint)/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C +$(derivedPoint)/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C +$(derivedPoint)/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C +$(derivedPoint)/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C +$(derivedPoint)/waveDisplacement/waveDisplacementPointPatchVectorField.C + +$(derivedPoint)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C + -pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C -pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C -pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C -pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C -pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C LIB = $(FOAM_LIBBIN)/libfvMotionSolvers diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C new file mode 100644 index 0000000000..e1bdfc0910 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C @@ -0,0 +1,482 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "timeVaryingMappedFixedValuePointPatchField.H" +#include "Time.H" +#include "AverageIOField.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam:: +timeVaryingMappedFixedValuePointPatchField:: +timeVaryingMappedFixedValuePointPatchField +( + const pointPatch& p, + const DimensionedField& iF +) +: + fixedValuePointPatchField(p, iF), + fieldTableName_(iF.name()), + setAverage_(false), + perturb_(0), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits::zero) + +{} + + +template +Foam:: +timeVaryingMappedFixedValuePointPatchField:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField& ptf, + const pointPatch& p, + const DimensionedField& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField(ptf, p, iF, mapper), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits::zero) +{} + + +template +Foam:: +timeVaryingMappedFixedValuePointPatchField:: +timeVaryingMappedFixedValuePointPatchField +( + const pointPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValuePointPatchField(p, iF), + fieldTableName_(iF.name()), + setAverage_(readBool(dict.lookup("setAverage"))), + perturb_(dict.lookupOrDefault("perturb", 1E-5)), + mapperPtr_(NULL), + sampleTimes_(0), + startSampleTime_(-1), + startSampledValues_(0), + startAverage_(pTraits::zero), + endSampleTime_(-1), + endSampledValues_(0), + endAverage_(pTraits::zero) +{ + dict.readIfPresent("fieldTableName", fieldTableName_); + updateCoeffs(); +} + + +template +Foam:: +timeVaryingMappedFixedValuePointPatchField:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField& ptf +) +: + fixedValuePointPatchField(ptf), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(ptf.sampleTimes_), + startSampleTime_(ptf.startSampleTime_), + startSampledValues_(ptf.startSampledValues_), + startAverage_(ptf.startAverage_), + endSampleTime_(ptf.endSampleTime_), + endSampledValues_(ptf.endSampledValues_), + endAverage_(ptf.endAverage_) +{} + + +template +Foam:: +timeVaryingMappedFixedValuePointPatchField:: +timeVaryingMappedFixedValuePointPatchField +( + const timeVaryingMappedFixedValuePointPatchField& ptf, + const DimensionedField& iF +) +: + fixedValuePointPatchField(ptf, iF), + fieldTableName_(ptf.fieldTableName_), + setAverage_(ptf.setAverage_), + perturb_(ptf.perturb_), + mapperPtr_(ptf.mapperPtr_), + sampleTimes_(ptf.sampleTimes_), + startSampleTime_(ptf.startSampleTime_), + startSampledValues_(ptf.startSampledValues_), + startAverage_(ptf.startAverage_), + endSampleTime_(ptf.endSampleTime_), + endSampledValues_(ptf.endSampledValues_), + endAverage_(ptf.endAverage_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::timeVaryingMappedFixedValuePointPatchField::checkTable() +{ + // Initialise + if (startSampleTime_ == -1 && endSampleTime_ == -1) + { + const polyMesh& pMesh = this->patch().boundaryMesh().mesh()(); + + // Read the initial point position + pointField meshPts; + + if (pMesh.pointsInstance() == pMesh.facesInstance()) + { + meshPts = pointField(pMesh.points(), this->patch().meshPoints()); + } + else + { + // Load points from facesInstance + if (debug) + { + Info<< "Reloading points0 from " << pMesh.facesInstance() + << endl; + } + + pointIOField points0 + ( + IOobject + ( + "points", + pMesh.facesInstance(), + polyMesh::meshSubDir, + pMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) + ); + meshPts = pointField(points0, this->patch().meshPoints()); + } + + pointIOField samplePoints + ( + IOobject + ( + "points", + this->db().time().constant(), + "boundaryData"/this->patch().name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + mapperPtr_.reset + ( + new pointToPointPlanarInterpolation + ( + samplePoints, + meshPts, + perturb_ + ) + ); + + // Read the times for which data is available + + const fileName samplePointsFile = samplePoints.filePath(); + const fileName samplePointsDir = samplePointsFile.path(); + sampleTimes_ = Time::findTimes(samplePointsDir); + + if (debug) + { + Info<< "timeVaryingMappedFixedValuePointPatchField : In directory " + << samplePointsDir << " found times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) + << endl; + } + } + + // Find current time in sampleTimes + label lo = -1; + label hi = -1; + + bool foundTime = mapperPtr_().findTime + ( + sampleTimes_, + startSampleTime_, + this->db().time().value(), + lo, + hi + ); + + if (!foundTime) + { + FatalErrorIn + ( + "timeVaryingMappedFixedValuePointPatchField::checkTable" + ) << "Cannot find starting sampling values for current time " + << this->db().time().value() << nl + << "Have sampling values for times " + << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl + << "In directory " + << this->db().time().constant()/"boundaryData"/this->patch().name() + << "\n on patch " << this->patch().name() + << " of field " << fieldTableName_ + << exit(FatalError); + } + + + // Update sampled data fields. + + if (lo != startSampleTime_) + { + startSampleTime_ = lo; + + if (startSampleTime_ == endSampleTime_) + { + // No need to reread since are end values + if (debug) + { + Pout<< "checkTable : Setting startValues to (already read) " + << "boundaryData" + /this->patch().name() + /sampleTimes_[startSampleTime_].name() + << endl; + } + startSampledValues_ = endSampledValues_; + startAverage_ = endAverage_; + } + else + { + if (debug) + { + Pout<< "checkTable : Reading startValues from " + << "boundaryData" + /this->patch().name() + /sampleTimes_[lo].name() + << endl; + } + + // Reread values and interpolate + AverageIOField vals + ( + IOobject + ( + fieldTableName_, + this->db().time().constant(), + "boundaryData" + /this->patch().name() + /sampleTimes_[startSampleTime_].name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + + startAverage_ = vals.average(); + startSampledValues_ = mapperPtr_().interpolate(vals); + } + } + + if (hi != endSampleTime_) + { + endSampleTime_ = hi; + + if (endSampleTime_ == -1) + { + // endTime no longer valid. Might as well clear endValues. + if (debug) + { + Pout<< "checkTable : Clearing endValues" << endl; + } + endSampledValues_.clear(); + } + else + { + if (debug) + { + Pout<< "checkTable : Reading endValues from " + << "boundaryData" + /this->patch().name() + /sampleTimes_[endSampleTime_].name() + << endl; + } + // Reread values and interpolate + AverageIOField vals + ( + IOobject + ( + fieldTableName_, + this->db().time().constant(), + "boundaryData" + /this->patch().name() + /sampleTimes_[endSampleTime_].name(), + this->db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE, + false + ) + ); + endAverage_ = vals.average(); + endSampledValues_ = mapperPtr_().interpolate(vals); + } + } +} + + +template +void Foam::timeVaryingMappedFixedValuePointPatchField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + checkTable(); + + // Interpolate between the sampled data + + Type wantedAverage; + + if (endSampleTime_ == -1) + { + // only start value + if (debug) + { + Pout<< "updateCoeffs : Sampled, non-interpolated values" + << " from start time:" + << sampleTimes_[startSampleTime_].name() << nl; + } + + this->operator==(startSampledValues_); + wantedAverage = startAverage_; + } + else + { + scalar start = sampleTimes_[startSampleTime_].value(); + scalar end = sampleTimes_[endSampleTime_].value(); + + scalar s = (this->db().time().value()-start)/(end-start); + + if (debug) + { + Pout<< "updateCoeffs : Sampled, interpolated values" + << " between start time:" + << sampleTimes_[startSampleTime_].name() + << " and end time:" << sampleTimes_[endSampleTime_].name() + << " with weight:" << s << endl; + } + + this->operator==((1-s)*startSampledValues_ + s*endSampledValues_); + wantedAverage = (1-s)*startAverage_ + s*endAverage_; + } + + // Enforce average. Either by scaling (if scaling factor > 0.5) or by + // offsetting. + if (setAverage_) + { + const Field& fld = *this; + + Type averagePsi = gAverage(fld); + + if (debug) + { + Pout<< "updateCoeffs :" + << " actual average:" << averagePsi + << " wanted average:" << wantedAverage + << endl; + } + + if (mag(averagePsi) < VSMALL) + { + // Field too small to scale. Offset instead. + const Type offset = wantedAverage - averagePsi; + if (debug) + { + Pout<< "updateCoeffs :" + << " offsetting with:" << offset << endl; + } + this->operator==(fld+offset); + } + else + { + const scalar scale = mag(wantedAverage)/mag(averagePsi); + + if (debug) + { + Pout<< "updateCoeffs :" + << " scaling with:" << scale << endl; + } + this->operator==(scale*fld); + } + } + + if (debug) + { + Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this) + << " max:" << gMax(*this) << endl; + } + + fixedValuePointPatchField::updateCoeffs(); +} + + +template +void Foam::timeVaryingMappedFixedValuePointPatchField::write +( + Ostream& os +) const +{ + fixedValuePointPatchField::write(os); + os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl; + os.writeKeyword("perturb") << perturb_ << token::END_STATEMENT << nl; + + if (fieldTableName_ != this->dimensionedInternalField().name()) + { + os.writeKeyword("fieldTableName") << fieldTableName_ + << token::END_STATEMENT << nl; + } +} + + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H new file mode 100644 index 0000000000..108ec6786e --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H @@ -0,0 +1,195 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +Class + Foam::timeVaryingMappedFixedValuePointPatchField + +Description + A time-varying form of a mapped fixed value boundary condition. + +See Also + Foam::timeVaryingMappedFixedValueFvPatchField + +SourceFiles + timeVaryingMappedFixedValuePointPatchField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingMappedFixedValuePointPatchField_H +#define timeVaryingMappedFixedValuePointPatchField_H + +#include "fixedValuePointPatchField.H" +#include "instantList.H" +#include "pointToPointPlanarInterpolation.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class timeVaryingMappedFixedValuePointPatchField Declaration +\*---------------------------------------------------------------------------*/ + +template +class timeVaryingMappedFixedValuePointPatchField +: + public fixedValuePointPatchField +{ + // Private data + + //- Name of the field data table, defaults to the name of the field + word fieldTableName_; + + //- If true adjust the mapped field to maintain average value + bool setAverage_; + + //- Fraction of perturbation (fraction of bounding box) to add + scalar perturb_; + + //- 2D interpolation + autoPtr mapperPtr_; + + //- List of boundaryData time directories + instantList sampleTimes_; + + //- Current starting index in sampleTimes + label startSampleTime_; + + //- Interpolated values from startSampleTime + Field startSampledValues_; + + //- If setAverage: starting average value + Type startAverage_; + + //- Current end index in sampleTimes + label endSampleTime_; + + //- Interpolated values from endSampleTime + Field endSampledValues_; + + //- If setAverage: end average value + Type endAverage_; + + +public: + + //- Runtime type information + TypeName("timeVaryingMappedFixedValue"); + + + // Constructors + + //- Construct from patch and internal field + timeVaryingMappedFixedValuePointPatchField + ( + const pointPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + timeVaryingMappedFixedValuePointPatchField + ( + const pointPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given patch field onto a new patch + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField&, + const pointPatch&, + const DimensionedField&, + const pointPatchFieldMapper& + ); + + //- Construct as copy + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField& + ); + + //- Construct and return a clone + virtual autoPtr > clone() const + { + return autoPtr > + ( + new timeVaryingMappedFixedValuePointPatchField(*this) + ); + } + + //- Construct as copy setting internal field reference + timeVaryingMappedFixedValuePointPatchField + ( + const timeVaryingMappedFixedValuePointPatchField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr > clone + ( + const DimensionedField& iF + ) const + { + return autoPtr > + ( + new timeVaryingMappedFixedValuePointPatchField(*this, iF) + ); + } + + + // Member functions + + // Utility functions + + //- Find boundary data inbetween current time and interpolate + void checkTable(); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "timeVaryingMappedFixedValuePointPatchField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C new file mode 100644 index 0000000000..ed847ab52d --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C @@ -0,0 +1,43 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "timeVaryingMappedFixedValuePointPatchFields.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makePointPatchFields(timeVaryingMappedFixedValue); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H new file mode 100644 index 0000000000..ad06b48066 --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +InClass + Foam::timeVaryingMappedFixedValuePointPatchFields + +Description + +SourceFiles + timeVaryingMappedFixedValuePointPatchFields.C + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingMappedFixedValuePointPatchFields_H +#define timeVaryingMappedFixedValuePointPatchFields_H + +#include "timeVaryingMappedFixedValuePointPatchField.H" +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchFieldTypedefs(timeVaryingMappedFixedValue); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 5912ea12e4..3329c2b28d 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -158,6 +158,7 @@ triSurface/triangleFuncs/triangleFuncs.C triSurface/surfaceFeatures/surfaceFeatures.C triSurface/triSurfaceTools/triSurfaceTools.C triSurface/triSurfaceTools/geompack/geompack.C +triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C twoDPointCorrector/twoDPointCorrector.C diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C new file mode 100644 index 0000000000..c60be140af --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C @@ -0,0 +1,322 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "pointToPointPlanarInterpolation.H" +#include "boundBox.H" +#include "Random.H" +#include "vector2D.H" +#include "triSurface.H" +#include "triSurfaceTools.H" +#include "OFstream.H" +#include "Time.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::pointToPointPlanarInterpolation, 0); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::coordinateSystem +Foam::pointToPointPlanarInterpolation::calcCoordinateSystem +( + const pointField& points +) const +{ + if (points.size() < 3) + { + FatalErrorIn + ( + "pointToPointPlanarInterpolation::calcCoordinateSystem" + "(const pointField&)" + ) << "Only " << points.size() << " provided." << nl + << "Need at least three non-colinear points" + << " to be able to interpolate." + << exit(FatalError); + } + + const point& p0 = points[0]; + + // Find furthest away point + vector e1; + label index1 = -1; + scalar maxDist = -GREAT; + + for (label i = 1; i < points.size(); i++) + { + const vector d = points[i] - p0; + scalar magD = mag(d); + + if (magD > maxDist) + { + e1 = d/magD; + index1 = i; + maxDist = magD; + } + } + // Find point that is furthest away from line p0-p1 + const point& p1 = points[index1]; + + label index2 = -1; + maxDist = -GREAT; + for (label i = 1; i < points.size(); i++) + { + if (i != index1) + { + const point& p2 = points[i]; + vector e2(p2 - p0); + e2 -= (e2&e1)*e1; + scalar magE2 = mag(e2); + + if (magE2 > maxDist) + { + index2 = i; + maxDist = magE2; + } + } + } + if (index2 == -1) + { + FatalErrorIn + ( + "pointToPointPlanarInterpolation::calcCoordinateSystem" + "(const pointField&)" + ) << "Cannot find points that make valid normal." << nl + << "Have so far points " << p0 << " and " << p1 + << "Need at least three points which are not in a line." + << exit(FatalError); + } + + vector n = e1^(points[index2]-p0); + n /= mag(n); + + if (debug) + { + Info<< "pointToPointPlanarInterpolation::calcCoordinateSystem :" + << " Used points " << p0 << ' ' << points[index1] + << ' ' << points[index2] + << " to define coordinate system with normal " << n << endl; + } + + return coordinateSystem + ( + "reference", + p0, // origin + n, // normal + e1 // 0-axis + ); +} + + +void Foam::pointToPointPlanarInterpolation::calcWeights +( + const pointField& sourcePoints, + const pointField& destPoints +) +{ + tmp tlocalVertices + ( + referenceCS_.localPosition(sourcePoints) + ); + vectorField& localVertices = tlocalVertices(); + + const boundBox bb(localVertices, true); + const point bbMid(bb.midpoint()); + + if (debug) + { + Info<< "pointToPointPlanarInterpolation::readData :" + << " Perturbing points with " << perturb_ + << " fraction of a random position inside " << bb + << " to break any ties on regular meshes." + << nl << endl; + } + + Random rndGen(123456); + forAll(localVertices, i) + { + localVertices[i] += + perturb_ + *(rndGen.position(bb.min(), bb.max())-bbMid); + } + + // Determine triangulation + List localVertices2D(localVertices.size()); + forAll(localVertices, i) + { + localVertices2D[i][0] = localVertices[i][0]; + localVertices2D[i][1] = localVertices[i][1]; + } + + triSurface s(triSurfaceTools::delaunay2D(localVertices2D)); + + tmp tlocalFaceCentres + ( + referenceCS_.localPosition + ( + destPoints + ) + ); + const pointField& localFaceCentres = tlocalFaceCentres(); + + if (debug) + { + Pout<< "pointToPointPlanarInterpolation::readData :" + <<" Dumping triangulated surface to triangulation.stl" << endl; + s.write("triangulation.stl"); + + OFstream str("localFaceCentres.obj"); + Pout<< "readSamplePoints :" + << " Dumping face centres to " << str.name() << endl; + + forAll(localFaceCentres, i) + { + const point& p = localFaceCentres[i]; + str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl; + } + } + + // Determine interpolation onto face centres. + triSurfaceTools::calcInterpolationWeights + ( + s, + localFaceCentres, // points to interpolate to + nearestVertex_, + nearestVertexWeight_ + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation +( + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb +) +: + perturb_(perturb), + referenceCS_(calcCoordinateSystem(sourcePoints)) + +{ + calcWeights(sourcePoints, destPoints); +} + + +Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation +( + const coordinateSystem& referenceCS, + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb +) +: + perturb_(perturb), + referenceCS_(referenceCS) +{ + calcWeights(sourcePoints, destPoints); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::wordList Foam::pointToPointPlanarInterpolation::timeNames +( + const instantList& times +) +{ + wordList names(times.size()); + + forAll(times, i) + { + names[i] = times[i].name(); + } + return names; +} + + +bool Foam::pointToPointPlanarInterpolation::findTime +( + const instantList& times, + const label startSampleTime, + const scalar timeVal, + label& lo, + label& hi +) +{ + lo = startSampleTime; + hi = -1; + + for (label i = startSampleTime+1; i < times.size(); i++) + { + if (times[i].value() > timeVal) + { + break; + } + else + { + lo = i; + } + } + + if (lo == -1) + { + //FatalErrorIn("findTime(..)") + // << "Cannot find starting sampling values for current time " + // << timeVal << nl + // << "Have sampling values for times " + // << timeNames(times) << nl + // << exit(FatalError); + return false; + } + + if (lo < times.size()-1) + { + hi = lo+1; + } + + + if (debug) + { + if (hi == -1) + { + Pout<< "findTime : Found time " << timeVal << " after" + << " index:" << lo << " time:" << times[lo].value() + << endl; + } + else + { + Pout<< "findTime : Found time " << timeVal << " inbetween" + << " index:" << lo << " time:" << times[lo].value() + << " and index:" << hi << " time:" << times[hi].value() + << endl; + } + } + return true; +} + + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H new file mode 100644 index 0000000000..daec2e4688 --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H @@ -0,0 +1,164 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +Class + Foam::pointToPointPlanarInterpolation + +Description + Interpolates between two sets of unstructured points using 2D Delaunay + triangulation. Used in e.g. timeVaryingMapped bcs. + +SourceFiles + pointToPointPlanarInterpolation.C + +\*---------------------------------------------------------------------------*/ + +#ifndef pointToPointPlanarInterpolation_H +#define pointToPointPlanarInterpolation_H + +#include "FixedList.H" +#include "coordinateSystem.H" +#include "instantList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pointToPointPlanarInterpolation Declaration +\*---------------------------------------------------------------------------*/ + +class pointToPointPlanarInterpolation +{ + // Private data + + //- Perturbation factor + const scalar perturb_; + + //- Coordinate system + coordinateSystem referenceCS_; + + //- Current interpolation addressing to face centres of underlying + // patch + List > nearestVertex_; + + //- Current interpolation factors to face centres of underlying + // patch + List > nearestVertexWeight_; + + // Private Member Functions + + //- Calculate a local coordinate system from set of points + coordinateSystem calcCoordinateSystem(const pointField&) const; + + //- Calculate addressing and weights + void calcWeights + ( + const pointField& sourcePoints, + const pointField& destPoints + ); + +public: + + // Declare name of the class and its debug switch + ClassName("pointToPointPlanarInterpolation"); + + + // Constructors + + //- Construct from 3D locations. Determines local coordinate system + // from sourcePoints and maps onto that. + pointToPointPlanarInterpolation + ( + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb + ); + + //- Construct from coordinate system and locations. + pointToPointPlanarInterpolation + ( + const coordinateSystem& referenceCS, + const pointField& sourcePoints, + const pointField& destPoints, + const scalar perturb + ); + + + // Member Functions + + //- Return the coordinateSystem + const coordinateSystem& referenceCS() const + { + return referenceCS_; + } + + // patch + const List >& nearestVertex() const + { + return nearestVertex_; + } + + //- Current interpolation factors to face centres of underlying + // patch + const List >& nearestVertexWeight() const + { + return nearestVertexWeight_; + } + + //- Helper: extract words of times + static wordList timeNames(const instantList&); + + //- Helper: find time. Return true if succesful. + static bool findTime + ( + const instantList& times, + const label startSampleTime, + const scalar timeVal, + label& lo, + label& hi + ); + + //- Interpolate from field on source points to dest points + template + tmp > interpolate(const Field& sourceFld) const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "pointToPointPlanarInterpolationTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C new file mode 100644 index 0000000000..9159a03edd --- /dev/null +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolationTemplates.C @@ -0,0 +1,71 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "pointToPointPlanarInterpolation.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::tmp > Foam::pointToPointPlanarInterpolation::interpolate +( + const Field& sourceFld +) const +{ + tmp > tfld(new Field(nearestVertex_.size())); + Field& fld = tfld(); + + forAll(fld, i) + { + const FixedList& verts = nearestVertex_[i]; + const FixedList& w = nearestVertexWeight_[i]; + + if (verts[2] == -1) + { + if (verts[1] == -1) + { + // Use vertex0 only + fld[i] = sourceFld[verts[0]]; + } + else + { + // Use vertex 0,1 + fld[i] = + w[0]*sourceFld[verts[0]] + + w[1]*sourceFld[verts[1]]; + } + } + else + { + fld[i] = + w[0]*sourceFld[verts[0]] + + w[1]*sourceFld[verts[1]] + + w[2]*sourceFld[verts[2]]; + } + } + return tfld; +} + + +// ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H index 2328aa8284..ede548fd83 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,7 @@ Class Foam::triSurfaceTools Description - A collection of tools for triSurfaceMesh + A collection of tools for triSurface. SourceFiles triSurfaceTools.C From 25dafe92eabc66b731877ba2f13dd28c9fc8fd21 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 16 Apr 2012 11:36:13 +0100 Subject: [PATCH 02/82] ENH: decomposePar: add -time option. --- .../decomposePar/decomposePar.C | 1116 ++++++++--------- 1 file changed, 557 insertions(+), 559 deletions(-) diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C index 02acd59236..c7b96b44b3 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C +++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,7 +43,10 @@ Usage Copy any \a uniform directories too. \param -constant \n - Override controlDict settings and use constant directory. + \param -time xxx:yyy \n + Override controlDict settings and decompose selected times. Does not + re-decompose the mesh i.e. does not handle moving mesh or changing + mesh cases. \param -fields \n Use existing geometry decomposition and convert fields only. @@ -123,11 +126,9 @@ int main(int argc, char *argv[]) "ifRequired", "only decompose geometry if the number of domains has changed" ); - argList::addBoolOption - ( - "constant", - "include the 'constant/' dir in the times list" - ); + + // Include explicit constant options, have zero from time range + timeSelector::addOptions(true, false); #include "setRootCase.H" @@ -146,23 +147,10 @@ int main(int argc, char *argv[]) bool forceOverwrite = args.optionFound("force"); bool ifRequiredDecomposition = args.optionFound("ifRequired"); + // Set time from database #include "createTime.H" - - // Allow -constant to override controlDict settings. - if (args.optionFound("constant")) - { - instantList timeDirs = timeSelector::select0(runTime, args); - if (runTime.timeName() != runTime.constant()) - { - FatalErrorIn(args.executable()) - << "No '" << runTime.constant() << "' time present." << endl - << "Valid times are " << runTime.times() - << exit(FatalError); - } - } - - - Info<< "Time = " << runTime.timeName() << endl; + // Allow override of time + instantList times = timeSelector::selectIfPresent(runTime, args); // determine the existing processor count directly label nProcs = 0; @@ -334,475 +322,394 @@ int main(int argc, char *argv[]) } - // Search for list of objects for this time - IOobjectList objects(mesh, runTime.timeName()); - - // Construct the vol fields - // ~~~~~~~~~~~~~~~~~~~~~~~~ - PtrList volScalarFields; - readFields(mesh, objects, volScalarFields); - - PtrList volVectorFields; - readFields(mesh, objects, volVectorFields); - - PtrList volSphericalTensorFields; - readFields(mesh, objects, volSphericalTensorFields); - - PtrList volSymmTensorFields; - readFields(mesh, objects, volSymmTensorFields); - - PtrList volTensorFields; - readFields(mesh, objects, volTensorFields); - - - // Construct the dimensioned fields - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - PtrList > dimScalarFields; - readFields(mesh, objects, dimScalarFields); - - PtrList > dimVectorFields; - readFields(mesh, objects, dimVectorFields); - - PtrList > - dimSphericalTensorFields; - readFields(mesh, objects, dimSphericalTensorFields); - - PtrList > dimSymmTensorFields; - readFields(mesh, objects, dimSymmTensorFields); - - PtrList > dimTensorFields; - readFields(mesh, objects, dimTensorFields); - - - // Construct the surface fields - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - PtrList surfaceScalarFields; - readFields(mesh, objects, surfaceScalarFields); - PtrList surfaceVectorFields; - readFields(mesh, objects, surfaceVectorFields); - PtrList surfaceSphericalTensorFields; - readFields(mesh, objects, surfaceSphericalTensorFields); - PtrList surfaceSymmTensorFields; - readFields(mesh, objects, surfaceSymmTensorFields); - PtrList surfaceTensorFields; - readFields(mesh, objects, surfaceTensorFields); - - - // Construct the point fields - // ~~~~~~~~~~~~~~~~~~~~~~~~~~ - pointMesh pMesh(mesh); - - PtrList pointScalarFields; - readFields(pMesh, objects, pointScalarFields); - - PtrList pointVectorFields; - readFields(pMesh, objects, pointVectorFields); - - PtrList pointSphericalTensorFields; - readFields(pMesh, objects, pointSphericalTensorFields); - - PtrList pointSymmTensorFields; - readFields(pMesh, objects, pointSymmTensorFields); - - PtrList pointTensorFields; - readFields(pMesh, objects, pointTensorFields); - - - // Construct the Lagrangian fields - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - fileNameList cloudDirs - ( - readDir(runTime.timePath()/cloud::prefix, fileName::DIRECTORY) - ); - - // Particles - PtrList > lagrangianPositions(cloudDirs.size()); - // Particles per cell - PtrList< List*> > cellParticles(cloudDirs.size()); - - PtrList > lagrangianLabelFields(cloudDirs.size()); - PtrList > lagrangianLabelFieldFields - ( - cloudDirs.size() - ); - PtrList > lagrangianScalarFields(cloudDirs.size()); - PtrList > lagrangianScalarFieldFields - ( - cloudDirs.size() - ); - PtrList > lagrangianVectorFields(cloudDirs.size()); - PtrList > lagrangianVectorFieldFields - ( - cloudDirs.size() - ); - PtrList > lagrangianSphericalTensorFields - ( - cloudDirs.size() - ); - PtrList > - lagrangianSphericalTensorFieldFields(cloudDirs.size()); - PtrList > lagrangianSymmTensorFields - ( - cloudDirs.size() - ); - PtrList > - lagrangianSymmTensorFieldFields - ( - cloudDirs.size() - ); - PtrList > lagrangianTensorFields - ( - cloudDirs.size() - ); - PtrList > lagrangianTensorFieldFields - ( - cloudDirs.size() - ); - - label cloudI = 0; - - forAll(cloudDirs, i) + // Loop over all times + forAll(times, timeI) { - IOobjectList sprayObjs + runTime.setTime(times[timeI], timeI); + + Info<< "Time = " << runTime.timeName() << endl; + + // Search for list of objects for this time + IOobjectList objects(mesh, runTime.timeName()); + + + // Construct the vol fields + // ~~~~~~~~~~~~~~~~~~~~~~~~ + PtrList volScalarFields; + readFields(mesh, objects, volScalarFields); + PtrList volVectorFields; + readFields(mesh, objects, volVectorFields); + PtrList volSphericalTensorFields; + readFields(mesh, objects, volSphericalTensorFields); + PtrList volSymmTensorFields; + readFields(mesh, objects, volSymmTensorFields); + PtrList volTensorFields; + readFields(mesh, objects, volTensorFields); + + + // Construct the dimensioned fields + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + PtrList > dimScalarFields; + readFields(mesh, objects, dimScalarFields); + PtrList > dimVectorFields; + readFields(mesh, objects, dimVectorFields); + PtrList > + dimSphericalTensorFields; + readFields(mesh, objects, dimSphericalTensorFields); + PtrList > dimSymmTensorFields; + readFields(mesh, objects, dimSymmTensorFields); + PtrList > dimTensorFields; + readFields(mesh, objects, dimTensorFields); + + + // Construct the surface fields + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + PtrList surfaceScalarFields; + readFields(mesh, objects, surfaceScalarFields); + PtrList surfaceVectorFields; + readFields(mesh, objects, surfaceVectorFields); + PtrList surfaceSphericalTensorFields; + readFields(mesh, objects, surfaceSphericalTensorFields); + PtrList surfaceSymmTensorFields; + readFields(mesh, objects, surfaceSymmTensorFields); + PtrList surfaceTensorFields; + readFields(mesh, objects, surfaceTensorFields); + + + // Construct the point fields + // ~~~~~~~~~~~~~~~~~~~~~~~~~~ + pointMesh pMesh(mesh); + + PtrList pointScalarFields; + readFields(pMesh, objects, pointScalarFields); + PtrList pointVectorFields; + readFields(pMesh, objects, pointVectorFields); + PtrList pointSphericalTensorFields; + readFields(pMesh, objects, pointSphericalTensorFields); + PtrList pointSymmTensorFields; + readFields(pMesh, objects, pointSymmTensorFields); + PtrList pointTensorFields; + readFields(pMesh, objects, pointTensorFields); + + + // Construct the Lagrangian fields + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + fileNameList cloudDirs ( - mesh, - runTime.timeName(), - cloud::prefix/cloudDirs[i] + readDir(runTime.timePath()/cloud::prefix, fileName::DIRECTORY) ); - IOobject* positionsPtr = sprayObjs.lookup("positions"); + // Particles + PtrList > lagrangianPositions(cloudDirs.size()); + // Particles per cell + PtrList< List*> > cellParticles + ( + cloudDirs.size() + ); - if (positionsPtr) + PtrList > lagrangianLabelFields + ( + cloudDirs.size() + ); + PtrList > lagrangianLabelFieldFields + ( + cloudDirs.size() + ); + PtrList > lagrangianScalarFields + ( + cloudDirs.size() + ); + PtrList > lagrangianScalarFieldFields + ( + cloudDirs.size() + ); + PtrList > lagrangianVectorFields + ( + cloudDirs.size() + ); + PtrList > lagrangianVectorFieldFields + ( + cloudDirs.size() + ); + PtrList > + lagrangianSphericalTensorFields + ( + cloudDirs.size() + ); + PtrList > + lagrangianSphericalTensorFieldFields(cloudDirs.size()); + PtrList > lagrangianSymmTensorFields + ( + cloudDirs.size() + ); + PtrList > + lagrangianSymmTensorFieldFields + ( + cloudDirs.size() + ); + PtrList > lagrangianTensorFields + ( + cloudDirs.size() + ); + PtrList > lagrangianTensorFieldFields + ( + cloudDirs.size() + ); + + label cloudI = 0; + + forAll(cloudDirs, i) { - // Read lagrangian particles - // ~~~~~~~~~~~~~~~~~~~~~~~~~ - - Info<< "Identified lagrangian data set: " << cloudDirs[i] << endl; - - lagrangianPositions.set - ( - cloudI, - new Cloud - ( - mesh, - cloudDirs[i], - false - ) - ); - - - // Sort particles per cell - // ~~~~~~~~~~~~~~~~~~~~~~~ - - cellParticles.set - ( - cloudI, - new List*> - ( - mesh.nCells(), - static_cast*>(NULL) - ) - ); - - label i = 0; - - forAllIter - ( - Cloud, - lagrangianPositions[cloudI], - iter - ) - { - iter().index() = i++; - - label celli = iter().cell(); - - // Check - if (celli < 0 || celli >= mesh.nCells()) - { - FatalErrorIn(args.executable()) - << "Illegal cell number " << celli - << " for particle with index " << iter().index() - << " at position " << iter().position() << nl - << "Cell number should be between 0 and " - << mesh.nCells()-1 << nl - << "On this mesh the particle should be in cell " - << mesh.findCell(iter().position()) - << exit(FatalError); - } - - if (!cellParticles[cloudI][celli]) - { - cellParticles[cloudI][celli] = new SLList - (); - } - - cellParticles[cloudI][celli]->append(&iter()); - } - - // Read fields - // ~~~~~~~~~~~ - - IOobjectList lagrangianObjects + IOobjectList sprayObjs ( mesh, runTime.timeName(), - cloud::prefix/cloudDirs[cloudI] + cloud::prefix/cloudDirs[i] ); - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianLabelFields - ); + IOobject* positionsPtr = sprayObjs.lookup("positions"); - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianLabelFieldFields - ); + if (positionsPtr) + { + // Read lagrangian particles + // ~~~~~~~~~~~~~~~~~~~~~~~~~ - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianScalarFields - ); + Info<< "Identified lagrangian data set: " << cloudDirs[i] + << endl; - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianScalarFieldFields - ); - - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianVectorFields - ); - - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianVectorFieldFields - ); - - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianSphericalTensorFields - ); - - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianSphericalTensorFieldFields - ); - - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianSymmTensorFields - ); - - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianSymmTensorFieldFields - ); - - lagrangianFieldDecomposer::readFields - ( - cloudI, - lagrangianObjects, - lagrangianTensorFields - ); - - lagrangianFieldDecomposer::readFieldFields - ( - cloudI, - lagrangianObjects, - lagrangianTensorFieldFields - ); - - cloudI++; - } - } - - lagrangianPositions.setSize(cloudI); - cellParticles.setSize(cloudI); - lagrangianLabelFields.setSize(cloudI); - lagrangianLabelFieldFields.setSize(cloudI); - lagrangianScalarFields.setSize(cloudI); - lagrangianScalarFieldFields.setSize(cloudI); - lagrangianVectorFields.setSize(cloudI); - lagrangianVectorFieldFields.setSize(cloudI); - lagrangianSphericalTensorFields.setSize(cloudI); - lagrangianSphericalTensorFieldFields.setSize(cloudI); - lagrangianSymmTensorFields.setSize(cloudI); - lagrangianSymmTensorFieldFields.setSize(cloudI); - lagrangianTensorFields.setSize(cloudI); - lagrangianTensorFieldFields.setSize(cloudI); + lagrangianPositions.set + ( + cloudI, + new Cloud + ( + mesh, + cloudDirs[i], + false + ) + ); - // Any uniform data to copy/link? - fileName uniformDir("uniform"); + // Sort particles per cell + // ~~~~~~~~~~~~~~~~~~~~~~~ - if (isDir(runTime.timePath()/uniformDir)) - { - Info<< "Detected additional non-decomposed files in " - << runTime.timePath()/uniformDir - << endl; - } - else - { - uniformDir.clear(); - } + cellParticles.set + ( + cloudI, + new List*> + ( + mesh.nCells(), + static_cast*>(NULL) + ) + ); - Info<< endl; + label i = 0; - // split the fields over processors - for (label procI = 0; procI < mesh.nProcs(); procI++) - { - Info<< "Processor " << procI << ": field transfer" << endl; + forAllIter + ( + Cloud, + lagrangianPositions[cloudI], + iter + ) + { + iter().index() = i++; - // open the database - Time processorDb - ( - Time::controlDictName, - args.rootPath(), - args.caseName()/fileName(word("processor") + name(procI)) - ); + label celli = iter().cell(); - processorDb.setTime(runTime); + // Check + if (celli < 0 || celli >= mesh.nCells()) + { + FatalErrorIn(args.executable()) + << "Illegal cell number " << celli + << " for particle with index " << iter().index() + << " at position " << iter().position() << nl + << "Cell number should be between 0 and " + << mesh.nCells()-1 << nl + << "On this mesh the particle should be in cell " + << mesh.findCell(iter().position()) + << exit(FatalError); + } - // remove files remnants that can cause horrible problems - // - mut and nut are used to mark the new turbulence models, - // their existence prevents old models from being upgraded - { - fileName timeDir(processorDb.path()/processorDb.timeName()); + if (!cellParticles[cloudI][celli]) + { + cellParticles[cloudI][celli] = + new SLList(); + } - rm(timeDir/"mut"); - rm(timeDir/"nut"); + cellParticles[cloudI][celli]->append(&iter()); + } + + // Read fields + // ~~~~~~~~~~~ + + IOobjectList lagrangianObjects + ( + mesh, + runTime.timeName(), + cloud::prefix/cloudDirs[cloudI] + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianLabelFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianLabelFieldFields + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianScalarFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianScalarFieldFields + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianVectorFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianVectorFieldFields + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianSphericalTensorFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianSphericalTensorFieldFields + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianSymmTensorFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianSymmTensorFieldFields + ); + + lagrangianFieldDecomposer::readFields + ( + cloudI, + lagrangianObjects, + lagrangianTensorFields + ); + + lagrangianFieldDecomposer::readFieldFields + ( + cloudI, + lagrangianObjects, + lagrangianTensorFieldFields + ); + + cloudI++; + } } - // read the mesh - fvMesh procMesh - ( - IOobject - ( - regionName, - processorDb.timeName(), - processorDb - ) - ); + lagrangianPositions.setSize(cloudI); + cellParticles.setSize(cloudI); + lagrangianLabelFields.setSize(cloudI); + lagrangianLabelFieldFields.setSize(cloudI); + lagrangianScalarFields.setSize(cloudI); + lagrangianScalarFieldFields.setSize(cloudI); + lagrangianVectorFields.setSize(cloudI); + lagrangianVectorFieldFields.setSize(cloudI); + lagrangianSphericalTensorFields.setSize(cloudI); + lagrangianSphericalTensorFieldFields.setSize(cloudI); + lagrangianSymmTensorFields.setSize(cloudI); + lagrangianSymmTensorFieldFields.setSize(cloudI); + lagrangianTensorFields.setSize(cloudI); + lagrangianTensorFieldFields.setSize(cloudI); - labelIOList faceProcAddressing - ( - IOobject - ( - "faceProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - labelIOList cellProcAddressing - ( - IOobject - ( - "cellProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); + // Any uniform data to copy/link? + fileName uniformDir("uniform"); - labelIOList boundaryProcAddressing - ( - IOobject - ( - "boundaryProcAddressing", - procMesh.facesInstance(), - procMesh.meshSubDir, - procMesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - - // FV fields + if (isDir(runTime.timePath()/uniformDir)) { - fvFieldDecomposer fieldDecomposer - ( - mesh, - procMesh, - faceProcAddressing, - cellProcAddressing, - boundaryProcAddressing - ); - - fieldDecomposer.decomposeFields(volScalarFields); - fieldDecomposer.decomposeFields(volVectorFields); - fieldDecomposer.decomposeFields(volSphericalTensorFields); - fieldDecomposer.decomposeFields(volSymmTensorFields); - fieldDecomposer.decomposeFields(volTensorFields); - - fieldDecomposer.decomposeFields(surfaceScalarFields); - fieldDecomposer.decomposeFields(surfaceVectorFields); - fieldDecomposer.decomposeFields(surfaceSphericalTensorFields); - fieldDecomposer.decomposeFields(surfaceSymmTensorFields); - fieldDecomposer.decomposeFields(surfaceTensorFields); + Info<< "Detected additional non-decomposed files in " + << runTime.timePath()/uniformDir + << endl; + } + else + { + uniformDir.clear(); } - // Dimensioned fields + Info<< endl; + + // split the fields over processors + for (label procI = 0; procI < mesh.nProcs(); procI++) { - dimFieldDecomposer fieldDecomposer + Info<< "Processor " << procI << ": field transfer" << endl; + + // open the database + Time processorDb ( - mesh, - procMesh, - faceProcAddressing, - cellProcAddressing + Time::controlDictName, + args.rootPath(), + args.caseName()/fileName(word("processor") + name(procI)) ); - fieldDecomposer.decomposeFields(dimScalarFields); - fieldDecomposer.decomposeFields(dimVectorFields); - fieldDecomposer.decomposeFields(dimSphericalTensorFields); - fieldDecomposer.decomposeFields(dimSymmTensorFields); - fieldDecomposer.decomposeFields(dimTensorFields); - } + processorDb.setTime(runTime); + // remove files remnants that can cause horrible problems + // - mut and nut are used to mark the new turbulence models, + // their existence prevents old models from being upgraded + { + fileName timeDir(processorDb.path()/processorDb.timeName()); - // Point fields - if - ( - pointScalarFields.size() - || pointVectorFields.size() - || pointSphericalTensorFields.size() - || pointSymmTensorFields.size() - || pointTensorFields.size() - ) - { - labelIOList pointProcAddressing + rm(timeDir/"mut"); + rm(timeDir/"nut"); + } + + // read the mesh + fvMesh procMesh ( IOobject ( - "pointProcAddressing", + regionName, + processorDb.timeName(), + processorDb + ) + ); + + labelIOList faceProcAddressing + ( + IOobject + ( + "faceProcAddressing", procMesh.facesInstance(), procMesh.meshSubDir, procMesh, @@ -811,138 +718,229 @@ int main(int argc, char *argv[]) ) ); - pointMesh procPMesh(procMesh); - - pointFieldDecomposer fieldDecomposer + labelIOList cellProcAddressing ( - pMesh, - procPMesh, - pointProcAddressing, - boundaryProcAddressing + IOobject + ( + "cellProcAddressing", + procMesh.facesInstance(), + procMesh.meshSubDir, + procMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) ); - fieldDecomposer.decomposeFields(pointScalarFields); - fieldDecomposer.decomposeFields(pointVectorFields); - fieldDecomposer.decomposeFields(pointSphericalTensorFields); - fieldDecomposer.decomposeFields(pointSymmTensorFields); - fieldDecomposer.decomposeFields(pointTensorFields); - } + labelIOList boundaryProcAddressing + ( + IOobject + ( + "boundaryProcAddressing", + procMesh.facesInstance(), + procMesh.meshSubDir, + procMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); - - // If there is lagrangian data write it out - forAll(lagrangianPositions, cloudI) - { - if (lagrangianPositions[cloudI].size()) + // FV fields { - lagrangianFieldDecomposer fieldDecomposer + fvFieldDecomposer fieldDecomposer ( mesh, procMesh, faceProcAddressing, cellProcAddressing, - cloudDirs[cloudI], - lagrangianPositions[cloudI], - cellParticles[cloudI] + boundaryProcAddressing ); - // Lagrangian fields + fieldDecomposer.decomposeFields(volScalarFields); + fieldDecomposer.decomposeFields(volVectorFields); + fieldDecomposer.decomposeFields(volSphericalTensorFields); + fieldDecomposer.decomposeFields(volSymmTensorFields); + fieldDecomposer.decomposeFields(volTensorFields); + + fieldDecomposer.decomposeFields(surfaceScalarFields); + fieldDecomposer.decomposeFields(surfaceVectorFields); + fieldDecomposer.decomposeFields(surfaceSphericalTensorFields); + fieldDecomposer.decomposeFields(surfaceSymmTensorFields); + fieldDecomposer.decomposeFields(surfaceTensorFields); + } + + // Dimensioned fields + { + dimFieldDecomposer fieldDecomposer + ( + mesh, + procMesh, + faceProcAddressing, + cellProcAddressing + ); + + fieldDecomposer.decomposeFields(dimScalarFields); + fieldDecomposer.decomposeFields(dimVectorFields); + fieldDecomposer.decomposeFields(dimSphericalTensorFields); + fieldDecomposer.decomposeFields(dimSymmTensorFields); + fieldDecomposer.decomposeFields(dimTensorFields); + } + + + // Point fields + if + ( + pointScalarFields.size() + || pointVectorFields.size() + || pointSphericalTensorFields.size() + || pointSymmTensorFields.size() + || pointTensorFields.size() + ) + { + labelIOList pointProcAddressing + ( + IOobject + ( + "pointProcAddressing", + procMesh.facesInstance(), + procMesh.meshSubDir, + procMesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + pointMesh procPMesh(procMesh); + + pointFieldDecomposer fieldDecomposer + ( + pMesh, + procPMesh, + pointProcAddressing, + boundaryProcAddressing + ); + + fieldDecomposer.decomposeFields(pointScalarFields); + fieldDecomposer.decomposeFields(pointVectorFields); + fieldDecomposer.decomposeFields(pointSphericalTensorFields); + fieldDecomposer.decomposeFields(pointSymmTensorFields); + fieldDecomposer.decomposeFields(pointTensorFields); + } + + + // If there is lagrangian data write it out + forAll(lagrangianPositions, cloudI) + { + if (lagrangianPositions[cloudI].size()) { - fieldDecomposer.decomposeFields + lagrangianFieldDecomposer fieldDecomposer ( + mesh, + procMesh, + faceProcAddressing, + cellProcAddressing, cloudDirs[cloudI], - lagrangianLabelFields[cloudI] + lagrangianPositions[cloudI], + cellParticles[cloudI] ); - fieldDecomposer.decomposeFieldFields + + // Lagrangian fields + { + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianLabelFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianLabelFieldFields[cloudI] + ); + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianScalarFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianScalarFieldFields[cloudI] + ); + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianVectorFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianVectorFieldFields[cloudI] + ); + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianSphericalTensorFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianSphericalTensorFieldFields[cloudI] + ); + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianSymmTensorFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianSymmTensorFieldFields[cloudI] + ); + fieldDecomposer.decomposeFields + ( + cloudDirs[cloudI], + lagrangianTensorFields[cloudI] + ); + fieldDecomposer.decomposeFieldFields + ( + cloudDirs[cloudI], + lagrangianTensorFieldFields[cloudI] + ); + } + } + } + + + // Any non-decomposed data to copy? + if (uniformDir.size()) + { + const fileName timePath = processorDb.timePath(); + + if (copyUniform || mesh.distributed()) + { + cp ( - cloudDirs[cloudI], - lagrangianLabelFieldFields[cloudI] + runTime.timePath()/uniformDir, + timePath/uniformDir ); - fieldDecomposer.decomposeFields + } + else + { + // link with relative paths + const string parentPath = string("..")/".."; + + fileName currentDir(cwd()); + chDir(timePath); + ln ( - cloudDirs[cloudI], - lagrangianScalarFields[cloudI] - ); - fieldDecomposer.decomposeFieldFields - ( - cloudDirs[cloudI], - lagrangianScalarFieldFields[cloudI] - ); - fieldDecomposer.decomposeFields - ( - cloudDirs[cloudI], - lagrangianVectorFields[cloudI] - ); - fieldDecomposer.decomposeFieldFields - ( - cloudDirs[cloudI], - lagrangianVectorFieldFields[cloudI] - ); - fieldDecomposer.decomposeFields - ( - cloudDirs[cloudI], - lagrangianSphericalTensorFields[cloudI] - ); - fieldDecomposer.decomposeFieldFields - ( - cloudDirs[cloudI], - lagrangianSphericalTensorFieldFields[cloudI] - ); - fieldDecomposer.decomposeFields - ( - cloudDirs[cloudI], - lagrangianSymmTensorFields[cloudI] - ); - fieldDecomposer.decomposeFieldFields - ( - cloudDirs[cloudI], - lagrangianSymmTensorFieldFields[cloudI] - ); - fieldDecomposer.decomposeFields - ( - cloudDirs[cloudI], - lagrangianTensorFields[cloudI] - ); - fieldDecomposer.decomposeFieldFields - ( - cloudDirs[cloudI], - lagrangianTensorFieldFields[cloudI] + parentPath/runTime.timeName()/uniformDir, + uniformDir ); + chDir(currentDir); } } } - - - // Any non-decomposed data to copy? - if (uniformDir.size()) - { - const fileName timePath = processorDb.timePath(); - - if (copyUniform || mesh.distributed()) - { - cp - ( - runTime.timePath()/uniformDir, - timePath/uniformDir - ); - } - else - { - // link with relative paths - const string parentPath = string("..")/".."; - - fileName currentDir(cwd()); - chDir(timePath); - ln - ( - parentPath/runTime.timeName()/uniformDir, - uniformDir - ); - chDir(currentDir); - } - } } - Info<< "\nEnd.\n" << endl; return 0; From aa7c66a3ffad756e57f3eaa3dd8de2411a6f6c57 Mon Sep 17 00:00:00 2001 From: Henry Date: Mon, 16 Apr 2012 11:46:35 +0100 Subject: [PATCH 03/82] twoPhaseMixture: Added support for named phases e.g. in transportProperties phases (air water); --- .../multiphase/cavitatingFoam/createFields.H | 53 ++++++++----------- .../multiphase/interFoam/createFields.H | 16 +----- .../twoLiquidMixingFoam/createFields.H | 16 +----- .../twoPhaseMixture.C | 19 +++++-- .../twoPhaseMixture.H | 16 +++++- 5 files changed, 54 insertions(+), 66 deletions(-) diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H index 49c7de1473..dbacf1dbd3 100644 --- a/applications/solvers/multiphase/cavitatingFoam/createFields.H +++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H @@ -25,38 +25,6 @@ mesh ); - volScalarField gamma - ( - IOobject - ( - "gamma", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0)) - ); - gamma.oldTime(); - - Info<< "Creating compressibilityModel\n" << endl; - autoPtr psiModel = - barotropicCompressibilityModel::New - ( - thermodynamicProperties, - gamma - ); - - const volScalarField& psi = psiModel->psi(); - - rho == max - ( - psi*p - + (1.0 - gamma)*rhol0 - + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, - rhoMin - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -78,6 +46,27 @@ twoPhaseMixture twoPhaseProperties(U, phiv, "gamma"); + volScalarField& gamma(twoPhaseProperties.alpha1()); + gamma.oldTime(); + + Info<< "Creating compressibilityModel\n" << endl; + autoPtr psiModel = + barotropicCompressibilityModel::New + ( + thermodynamicProperties, + gamma + ); + + const volScalarField& psi = psiModel->psi(); + + rho == max + ( + psi*p + + (1.0 - gamma)*rhol0 + + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, + rhoMin + ); + // Create incompressible turbulence model autoPtr turbulence ( diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H index 34dea334e8..68765262db 100644 --- a/applications/solvers/multiphase/interFoam/createFields.H +++ b/applications/solvers/multiphase/interFoam/createFields.H @@ -12,20 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -46,6 +32,8 @@ Info<< "Reading transportProperties\n" << endl; twoPhaseMixture twoPhaseProperties(U, phi); + volScalarField& alpha1(twoPhaseProperties.alpha1()); + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H index 73af502ce2..0d01b9a9e5 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H @@ -12,20 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -45,6 +31,8 @@ Info<< "Reading transportProperties\n" << endl; twoPhaseMixture twoPhaseProperties(U, phi); + volScalarField& alpha1(twoPhaseProperties.alpha1()); + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C index 43f1ca8733..392704b940 100644 --- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C +++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,8 +59,8 @@ Foam::twoPhaseMixture::twoPhaseMixture : transportModel(U, phi), - phase1Name_("phase1"), - phase2Name_("phase2"), + phase1Name_(found("phases") ? wordList(lookup("phases"))[0] : "phase1"), + phase2Name_(found("phases") ? wordList(lookup("phases"))[1] : "phase2"), nuModel1_ ( @@ -89,7 +89,18 @@ Foam::twoPhaseMixture::twoPhaseMixture U_(U), phi_(phi), - alpha1_(U_.db().lookupObject (alpha1Name)), + alpha1_ + ( + IOobject + ( + found("phases") ? word("alpha" + phase1Name_) : alpha1Name, + U_.time().timeName(), + U_.db(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + U_.mesh() + ), nu_ ( diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H index 5e5eaaf10b..025d95f701 100644 --- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H +++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,7 +69,7 @@ protected: const volVectorField& U_; const surfaceScalarField& phi_; - const volScalarField& alpha1_; + volScalarField alpha1_; volScalarField nu_; @@ -110,6 +110,18 @@ public: return phase2Name_; } + //- Return the phase-fraction of phase 1 + const volScalarField& alpha1() const + { + return alpha1_; + } + + //- Return the phase-fraction of phase 1 + volScalarField& alpha1() + { + return alpha1_; + } + //- Return const-access to phase1 viscosityModel const viscosityModel& nuModel1() const { From 8e3e1808ec0ef83ed1229a8af70b6a3c46ab8258 Mon Sep 17 00:00:00 2001 From: Henry Date: Mon, 16 Apr 2012 11:47:28 +0100 Subject: [PATCH 04/82] compressibleInterFoam: Added thermal support --- .../compressibleInterFoam/Make/files | 1 + .../multiphase/compressibleInterFoam/TEqn.H | 20 ++ .../compressibleInterDyMFoam.C | 8 +- .../compressibleInterDyMFoam/pEqn.H | 88 -------- .../compressibleInterFoam.C | 5 +- .../compressibleInterFoam/createFields.H | 105 ++++++++-- .../wallHeatTransferFvPatchScalarField.C | 184 +++++++++++++++++ .../wallHeatTransferFvPatchScalarField.H | 194 ++++++++++++++++++ .../multiphase/compressibleInterFoam/pEqn.H | 61 +++--- .../les/depthCharge2D/0/T.org | 34 +++ .../0/{alpha1.org => alphawater.org} | 2 +- .../les/depthCharge2D/Allclean | 3 +- .../les/depthCharge2D/Allrun | 3 +- .../constant/polyMesh/blockMeshDict | 34 ++- .../constant/transportProperties | 14 +- .../les/depthCharge2D/system/controlDict | 2 +- .../les/depthCharge2D/system/fvSchemes | 5 +- .../les/depthCharge2D/system/fvSolution | 2 +- .../les/depthCharge2D/system/setFieldsDict | 8 +- .../les/depthCharge3D/0/T.org | 29 +++ .../0/{alpha1.org => alphawater.org} | 2 +- .../les/depthCharge3D/Allclean | 2 +- .../les/depthCharge3D/Allrun | 3 +- .../constant/transportProperties | 14 +- .../les/depthCharge3D/system/fvSchemes | 5 +- .../les/depthCharge3D/system/fvSolution | 2 +- .../les/depthCharge3D/system/setFieldsDict | 8 +- 27 files changed, 651 insertions(+), 187 deletions(-) create mode 100644 applications/solvers/multiphase/compressibleInterFoam/TEqn.H delete mode 100644 applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H create mode 100644 applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C create mode 100644 applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H create mode 100644 tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/T.org rename tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/{alpha1.org => alphawater.org} (97%) create mode 100644 tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/T.org rename tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/{alpha1.org => alphawater.org} (97%) diff --git a/applications/solvers/multiphase/compressibleInterFoam/Make/files b/applications/solvers/multiphase/compressibleInterFoam/Make/files index de5437219c..0e009f09bc 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/Make/files +++ b/applications/solvers/multiphase/compressibleInterFoam/Make/files @@ -1,3 +1,4 @@ +derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C compressibleInterFoam.C EXE = $(FOAM_APPBIN)/compressibleInterFoam diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H new file mode 100644 index 0000000000..34b6bd35dd --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H @@ -0,0 +1,20 @@ +{ + volScalarField kByCv + ( + "kByCv", + (alpha1*k1/Cv1 + alpha2*k2/Cv2) + + (alpha1*rho1 + alpha2*rho2)*turbulence->nut() + ); + + solve + ( + fvm::ddt(rho, T) + + fvm::div(rhoPhi, T) + - fvm::laplacian(kByCv, T) + + p*fvc::div(phi)*(alpha1/Cv1 + alpha2/Cv2) + ); + + // Update compressibilities + psi1 = 1.0/(R1*T); + psi2 = 1.0/(R2*T); +} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C index 1ac1596c4d..a4cfc215fc 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,7 @@ Application compressibleInterDyMFoam Description - Solver for 2 compressible, isothermal immiscible fluids using a VOF + Solver for 2 compressible, non-isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach, with optional mesh motion and mesh topology changes including adaptive re-meshing. @@ -124,11 +124,15 @@ int main(int argc, char *argv[]) solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" + #include "TEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); } } diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H deleted file mode 100644 index 26666c4120..0000000000 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ /dev/null @@ -1,88 +0,0 @@ -{ - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); - - tmp p_rghEqnComp; - - if (pimple.transonic()) - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvm::div(phi, p_rgh) - - fvm::Sp(fvc::div(phi), p_rgh) - ); - } - else - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvc::div(phi, p_rgh) - - fvc::Sp(fvc::div(phi), p_rgh) - ); - } - - - U = rAU*UEqn.H(); - - surfaceScalarField phiU - ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU, rho, U, phi) - ); - - phi = phiU + - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix p_rghEqnIncomp - ( - fvc::div(phi) - - fvm::laplacian(rAUf, p_rgh) - ); - - solve - ( - ( - max(alpha1, scalar(0))*(psi1/rho1) - + max(alpha2, scalar(0))*(psi2/rho2) - ) - *p_rghEqnComp() - + p_rghEqnIncomp, - mesh.solver(p_rgh.select(pimple.finalInnerIter())) - ); - - if (pimple.finalNonOrthogonalIter()) - { - dgdt = - (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) - *(p_rghEqnComp & p_rgh); - phi += p_rghEqnIncomp.flux(); - } - } - - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); - U.correctBoundaryConditions(); - - p = max - ( - (p_rgh + gh*(alpha1*rho10 + alpha2*rho20)) - /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)), - pMin - ); - - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; - - Info<< "max(U) " << max(mag(U)).value() << endl; - Info<< "min(p_rgh) " << min(p_rgh).value() << endl; - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); -} diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index 4ad1b3d01d..58dd7e6866 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,7 @@ Application compressibleInterFoam Description - Solver for 2 compressible, isothermal immiscible fluids using a VOF + Solver for 2 compressible, non-isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the "mixture" and a single @@ -82,6 +82,7 @@ int main(int argc, char *argv[]) solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" + #include "TEqn.H" // --- Pressure corrector loop while (pimple.correct()) diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H index c598cb75ce..4930743887 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H @@ -12,23 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< "Calculating field alpha1\n" << endl; - volScalarField alpha2("alpha2", scalar(1) - alpha1); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -45,10 +28,19 @@ #include "createPhi.H" - Info<< "Reading transportProperties\n" << endl; twoPhaseMixture twoPhaseProperties(U, phi); + volScalarField& alpha1(twoPhaseProperties.alpha1()); + + Info<< "Calculating phase-fraction alpha" << twoPhaseProperties.phase2Name() + << nl << endl; + volScalarField alpha2 + ( + "alpha" + twoPhaseProperties.phase2Name(), + scalar(1) - alpha1 + ); + dimensionedScalar rho10 ( twoPhaseProperties.subDict @@ -65,22 +57,91 @@ ).lookup("rho0") ); - dimensionedScalar psi1 + dimensionedScalar k1 ( twoPhaseProperties.subDict ( twoPhaseProperties.phase1Name() - ).lookup("psi") + ).lookup("k") ); - dimensionedScalar psi2 + dimensionedScalar k2 ( twoPhaseProperties.subDict ( twoPhaseProperties.phase2Name() - ).lookup("psi") + ).lookup("k") ); + dimensionedScalar Cv1 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ).lookup("Cv") + ); + + dimensionedScalar Cv2 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("Cv") + ); + + dimensionedScalar R1 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ).lookup("R") + ); + + dimensionedScalar R2 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("R") + ); + + Info<< "Reading field T\n" << endl; + volScalarField T + ( + IOobject + ( + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + volScalarField psi1 + ( + IOobject + ( + "psi1", + runTime.timeName(), + mesh + ), + 1.0/(R1*T) + ); + + volScalarField psi2 + ( + IOobject + ( + "psi2", + runTime.timeName(), + mesh + ), + 1.0/(R2*T) + ); + psi2.oldTime(); + dimensionedScalar pMin(twoPhaseProperties.lookup("pMin")); Info<< "Calculating field g.h\n" << endl; diff --git a/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C new file mode 100644 index 0000000000..e6782e8b3a --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C @@ -0,0 +1,184 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +#include "wallHeatTransferFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(p, iF), + Tinf_(p.size(), 0.0), + alphaWall_(p.size(), 0.0) +{ + refValue() = 0.0; + refGrad() = 0.0; + valueFraction() = 0.0; +} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + mixedFvPatchScalarField(ptf, p, iF, mapper), + Tinf_(ptf.Tinf_, mapper), + alphaWall_(ptf.alphaWall_, mapper) +{} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + mixedFvPatchScalarField(p, iF), + Tinf_("Tinf", dict, p.size()), + alphaWall_("alphaWall", dict, p.size()) +{ + refValue() = Tinf_; + refGrad() = 0.0; + valueFraction() = 0.0; + + if (dict.found("value")) + { + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); + } + else + { + evaluate(); + } +} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& tppsf +) +: + mixedFvPatchScalarField(tppsf), + Tinf_(tppsf.Tinf_), + alphaWall_(tppsf.alphaWall_) +{} + + +Foam::wallHeatTransferFvPatchScalarField::wallHeatTransferFvPatchScalarField +( + const wallHeatTransferFvPatchScalarField& tppsf, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(tppsf, iF), + Tinf_(tppsf.Tinf_), + alphaWall_(tppsf.alphaWall_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::wallHeatTransferFvPatchScalarField::autoMap +( + const fvPatchFieldMapper& m +) +{ + scalarField::autoMap(m); + Tinf_.autoMap(m); + alphaWall_.autoMap(m); +} + + +void Foam::wallHeatTransferFvPatchScalarField::rmap +( + const fvPatchScalarField& ptf, + const labelList& addr +) +{ + mixedFvPatchScalarField::rmap(ptf, addr); + + const wallHeatTransferFvPatchScalarField& tiptf = + refCast(ptf); + + Tinf_.rmap(tiptf.Tinf_, addr); + alphaWall_.rmap(tiptf.alphaWall_, addr); +} + + +void Foam::wallHeatTransferFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const fvPatchScalarField& Cpw = + patch().lookupPatchField("Cp"); + + const fvPatchScalarField& kByCpw = + patch().lookupPatchField("kByCp"); + + valueFraction() = + 1.0/ + ( + 1.0 + + Cpw*kByCpw*patch().deltaCoeffs()/alphaWall_ + ); + + mixedFvPatchScalarField::updateCoeffs(); +} + + +void Foam::wallHeatTransferFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + Tinf_.writeEntry("Tinf", os); + alphaWall_.writeEntry("alphaWall", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField(fvPatchScalarField, wallHeatTransferFvPatchScalarField); +} + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H new file mode 100644 index 0000000000..cf5284a087 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.H @@ -0,0 +1,194 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +Class + Foam::wallHeatTransferFvPatchScalarField + +Description + Enthalpy boundary conditions for wall heat transfer + +SourceFiles + wallHeatTransferFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef wallHeatTransferFvPatchScalarField_H +#define wallHeatTransferFvPatchScalarField_H + +#include "mixedFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class wallHeatTransferFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class wallHeatTransferFvPatchScalarField +: + public mixedFvPatchScalarField +{ + // Private data + + //- Tinf + scalarField Tinf_; + + //- alphaWall + scalarField alphaWall_; + + +public: + + //- Runtime type information + TypeName("wallHeatTransfer"); + + + // Constructors + + //- Construct from patch and internal field + wallHeatTransferFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + wallHeatTransferFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given wallHeatTransferFvPatchScalarField + // onto a new patch + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new wallHeatTransferFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + wallHeatTransferFvPatchScalarField + ( + const wallHeatTransferFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new wallHeatTransferFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Access + + //- Return Tinf + const scalarField& Tinf() const + { + return Tinf_; + } + + //- Return reference to Tinf to allow adjustment + scalarField& Tinf() + { + return Tinf_; + } + + //- Return alphaWall + const scalarField& alphaWall() const + { + return alphaWall_; + } + + //- Return reference to alphaWall to allow adjustment + scalarField& alphaWall() + { + return alphaWall_; + } + + + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchScalarField&, + const labelList& + ); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 035e8e237d..b803b1ac18 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -1,28 +1,31 @@ { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rAUf(fvc::interpolate(rAU)); + rho1 = rho10 + psi1*p; + rho2 = rho20 + psi2*p; - tmp p_rghEqnComp; + volScalarField rAU = 1.0/UEqn.A(); + surfaceScalarField rAUf = fvc::interpolate(rAU); - if (pimple.transonic()) + tmp p_rghEqnComp1; + tmp p_rghEqnComp2; + + //if (transonic) + //{ + //} + //else { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvm::div(phi, p_rgh) - - fvm::Sp(fvc::div(phi), p_rgh) - ); - } - else - { - p_rghEqnComp = - ( - fvm::ddt(p_rgh) - + fvc::div(phi, p_rgh) - - fvc::Sp(fvc::div(phi), p_rgh) - ); - } + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); + surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); + p_rghEqnComp1 = + fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh)) + + fvc::div(phid1, p_rgh) + - fvc::Sp(fvc::div(phid1), p_rgh); + + p_rghEqnComp2 = + fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh)) + + fvc::div(phid2, p_rgh) + - fvc::Sp(fvc::div(phid2), p_rgh); + } U = rAU*UEqn.H(); @@ -39,6 +42,10 @@ - ghf*fvc::snGrad(rho) )*rAUf*mesh.magSf(); + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + //thermo.rho() -= psi*p_rgh; + while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqnIncomp @@ -50,19 +57,23 @@ solve ( ( - max(alpha1, scalar(0))*(psi1/rho1) - + max(alpha2, scalar(0))*(psi2/rho2) + (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1() + + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2() ) - *p_rghEqnComp() + p_rghEqnIncomp, mesh.solver(p_rgh.select(pimple.finalInnerIter())) ); if (pimple.finalNonOrthogonalIter()) { + // Second part of thermodynamic density update + //thermo.rho() += psi*p_rgh; + dgdt = - (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) - *(p_rghEqnComp & p_rgh); + ( + pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 + - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 + ); phi += p_rghEqnIncomp.flux(); } } diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/T.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/T.org new file mode 100644 index 0000000000..e5ac2eeb27 --- /dev/null +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/T.org @@ -0,0 +1,34 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object alpha1; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + walls + { + type zeroGradient; + } + + defaultFaces + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alpha1.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alphawater.org similarity index 97% rename from tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alpha1.org rename to tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alphawater.org index b55820f9a9..d90b720092 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alpha1.org +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/0/alphawater.org @@ -10,7 +10,7 @@ FoamFile version 2.0; format ascii; class volScalarField; - object alpha1; + object alphawater; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allclean b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allclean index 1f81d5f86a..d6e8dce5df 100755 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allclean +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allclean @@ -1,5 +1,6 @@ #!/bin/sh +cd ${0%/*} || exit 1 # run from this directory foamCleanTutorials cases rm -rf processor* -rm -rf 0/p_rgh 0/p_rgh.gz 0/alpha1 0/alpha1.gz +rm -rf 0/p_rgh.gz 0/alphawater.gz 0/T.gz diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allrun b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allrun index 5fd844dc71..76547516b7 100755 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allrun +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/Allrun @@ -8,8 +8,9 @@ cd ${0%/*} || exit 1 # run from this directory application=`getApplication` runApplication blockMesh -cp 0/alpha1.org 0/alpha1 +cp 0/alphawater.org 0/alphawater cp 0/p_rgh.org 0/p_rgh +cp 0/T.org 0/T runApplication setFields runApplication $application diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/blockMeshDict b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/blockMeshDict index 22d11441f9..9ed9a630b5 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/blockMeshDict +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/polyMesh/blockMeshDict @@ -33,28 +33,20 @@ blocks hex (0 1 2 3 4 5 6 7) (80 160 1) simpleGrading (1 1 1) ); -boundary +patches ( - walls - { - type wall; - faces - ( - (3 7 6 2) - (0 4 7 3) - (2 6 5 1) - (1 5 4 0) - ); - } - frontAndBack - { - type empty; - faces - ( - (0 3 2 1) - (4 5 6 7) - ); - } + wall walls + ( + (3 7 6 2) + (0 4 7 3) + (2 6 5 1) + (1 5 4 0) + ) + empty frontAndBack + ( + (0 3 2 1) + (4 5 6 7) + ) ); // ************************************************************************* // diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/transportProperties b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/transportProperties index 67e66cd27d..1803697ff7 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/transportProperties +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/constant/transportProperties @@ -15,22 +15,28 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -phase1 +phases (water air); + +water { transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [1 1 -3 -1 0 0 0] 0; //0.613; rho rho [ 1 -3 0 0 0 0 0 ] 1000; rho0 rho0 [ 1 -3 0 0 0 0 0 ] 1000; - psi psi [ 0 -2 2 0 0 ] 1e-05; + R R [0 2 -2 -1 0] 3000; + Cv Cv [0 2 -2 -1 0] 4179; } -phase2 +air { transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1.589e-05; + k k [1 1 -3 -1 0 0 0] 0; //2.63e-2; rho rho [ 1 -3 0 0 0 0 0 ] 1; rho0 rho0 [ 1 -3 0 0 0 0 0 ] 0; - psi psi [ 0 -2 2 0 0 ] 1e-05; + R R [0 2 -2 -1 0] 287; + Cv Cv [0 2 -2 -1 0] 721; } pMin pMin [ 1 -1 -2 0 0 0 0 ] 10000; diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/controlDict b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/controlDict index 4c225df9c7..933d1de450 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/controlDict +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/controlDict @@ -17,7 +17,7 @@ FoamFile application compressibleInterFoam; -startFrom latestTime; +startFrom startTime; startTime 0; diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes index 5ab2bdc2fd..c35af6cc4f 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes @@ -30,7 +30,9 @@ divSchemes div(rho*phi,U) Gauss upwind; div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression 1; - div(phi,p_rgh) Gauss upwind; + div(phid1,p_rgh) Gauss upwind; + div(phid2,p_rgh) Gauss upwind; + div(rho*phi,T) Gauss upwind; div(phi,k) Gauss vanLeer; div((nuEff*dev(T(grad(U))))) Gauss linear; } @@ -55,7 +57,6 @@ fluxRequired default no; p_rgh; pcorr; - gamma; } diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSolution index 4a7f094cec..4618974546 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSolution +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSolution @@ -91,7 +91,7 @@ solvers nSweeps 1; } - "(k|B|nuTilda)" + "(T|k|B|nuTilda).*" { solver PBiCG; preconditioner DILU; diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/setFieldsDict b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/setFieldsDict index 37f90fd1d5..88198a4d24 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/setFieldsDict +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/setFieldsDict @@ -17,8 +17,9 @@ FoamFile defaultFieldValues ( - volScalarFieldValue alpha1 1 + volScalarFieldValue alphawater 1 volScalarFieldValue p_rgh 1e5 + volScalarFieldValue T 300 ); regions @@ -29,8 +30,9 @@ regions radius 0.1; fieldValues ( - volScalarFieldValue alpha1 0 + volScalarFieldValue alphawater 0 volScalarFieldValue p_rgh 1e6 + volScalarFieldValue T 578 ); } boxToCell @@ -38,7 +40,7 @@ regions box (-10 1 -1) (10 10 1); fieldValues ( - volScalarFieldValue alpha1 0 + volScalarFieldValue alphawater 0 ); } ); diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/T.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/T.org new file mode 100644 index 0000000000..b40cb08478 --- /dev/null +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/T.org @@ -0,0 +1,29 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object alpha1; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + walls + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alpha1.org b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alphawater.org similarity index 97% rename from tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alpha1.org rename to tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alphawater.org index b15de21b75..62be61f403 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alpha1.org +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/0/alphawater.org @@ -10,7 +10,7 @@ FoamFile version 2.0; format ascii; class volScalarField; - object alpha1.org; + object alphawater; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allclean b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allclean index 902055c7f7..539c772166 100755 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allclean +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allclean @@ -3,6 +3,6 @@ cd ${0%/*} || exit 1 # run from this directory foamCleanTutorials cases rm -rf processor* -rm -rf 0/p_rgh 0/p_rgh.gz 0/alpha1 0/alpha1.gz +rm -rf 0/p_rgh 0/p_rgh.gz 0/alphawater 0/alphawater.gz 0/T.gz # ----------------------------------------------------------------- end-of-file diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allrun b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allrun index ac00d007bc..0941ac4aa1 100755 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allrun +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/Allrun @@ -8,8 +8,9 @@ cd ${0%/*} || exit 1 # run from this directory application=`getApplication` runApplication blockMesh -cp 0/alpha1.org 0/alpha1 +cp 0/alphawater.org 0/alphawater cp 0/p_rgh.org 0/p_rgh +cp 0/T.org 0/T runApplication setFields runApplication decomposePar runParallel $application 4 diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/constant/transportProperties b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/constant/transportProperties index 67e66cd27d..1803697ff7 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/constant/transportProperties +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/constant/transportProperties @@ -15,22 +15,28 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -phase1 +phases (water air); + +water { transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1e-06; + k k [1 1 -3 -1 0 0 0] 0; //0.613; rho rho [ 1 -3 0 0 0 0 0 ] 1000; rho0 rho0 [ 1 -3 0 0 0 0 0 ] 1000; - psi psi [ 0 -2 2 0 0 ] 1e-05; + R R [0 2 -2 -1 0] 3000; + Cv Cv [0 2 -2 -1 0] 4179; } -phase2 +air { transportModel Newtonian; nu nu [ 0 2 -1 0 0 0 0 ] 1.589e-05; + k k [1 1 -3 -1 0 0 0] 0; //2.63e-2; rho rho [ 1 -3 0 0 0 0 0 ] 1; rho0 rho0 [ 1 -3 0 0 0 0 0 ] 0; - psi psi [ 0 -2 2 0 0 ] 1e-05; + R R [0 2 -2 -1 0] 287; + Cv Cv [0 2 -2 -1 0] 721; } pMin pMin [ 1 -1 -2 0 0 0 0 ] 10000; diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes index 5ab2bdc2fd..c35af6cc4f 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes @@ -30,7 +30,9 @@ divSchemes div(rho*phi,U) Gauss upwind; div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression 1; - div(phi,p_rgh) Gauss upwind; + div(phid1,p_rgh) Gauss upwind; + div(phid2,p_rgh) Gauss upwind; + div(rho*phi,T) Gauss upwind; div(phi,k) Gauss vanLeer; div((nuEff*dev(T(grad(U))))) Gauss linear; } @@ -55,7 +57,6 @@ fluxRequired default no; p_rgh; pcorr; - gamma; } diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSolution index 4a7f094cec..4618974546 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSolution +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSolution @@ -91,7 +91,7 @@ solvers nSweeps 1; } - "(k|B|nuTilda)" + "(T|k|B|nuTilda).*" { solver PBiCG; preconditioner DILU; diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/setFieldsDict b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/setFieldsDict index 07f675b5de..ec31deae03 100644 --- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/setFieldsDict +++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/setFieldsDict @@ -17,8 +17,9 @@ FoamFile defaultFieldValues ( - volScalarFieldValue alpha1 1 + volScalarFieldValue alphawater 1 volScalarFieldValue p_rgh 1e5 + volScalarFieldValue T 300 ); regions @@ -29,8 +30,9 @@ regions radius 0.1; fieldValues ( - volScalarFieldValue alpha1 0 + volScalarFieldValue alphawater 0 volScalarFieldValue p_rgh 1e6 + volScalarFieldValue T 578 ); } boxToCell @@ -38,7 +40,7 @@ regions box (-10 1 -1) (10 10 1); fieldValues ( - volScalarFieldValue alpha1 0 + volScalarFieldValue alphawater 0 ); } ); From 96f8104ae5fb46ef98181dd16795d25f102f69e4 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 16 Apr 2012 12:37:01 +0100 Subject: [PATCH 05/82] ENH: decomposePar: cache decomposeers if running with multiple times. --- .../decomposePar/decomposePar.C | 319 +++++++++++++----- 1 file changed, 227 insertions(+), 92 deletions(-) diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C index c7b96b44b3..49304aadcd 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C +++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C @@ -322,6 +322,23 @@ int main(int argc, char *argv[]) } + + // Caches + // ~~~~~~ + // Cached processor meshes and maps. These are only preserved if running + // with multiple times. + PtrList