From 7a2c87c10f347698ac4beb8a356d88c7d3fe50d0 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 1 Jun 2017 14:49:38 +0100 Subject: [PATCH] MPPIC: Optimised the averaging methods The averaging methods now take the particle barycentric coordinates as inputs rather than global positions. This change significantly optimises Dual averaging, which is the most commonly used method. The run time of the lagrangian/MPPICFoam/Goldschmidt tutorial has been reduced by a factor of about two. --- .../MPPICParcel/MPPICParcelTrackingDataI.H | 24 +++---- .../AveragingMethod/AveragingMethod.C | 20 ++++-- .../AveragingMethod/AveragingMethod.H | 9 +-- .../MPPIC/AveragingMethods/Basic/Basic.C | 8 +-- .../MPPIC/AveragingMethods/Basic/Basic.H | 8 +-- .../MPPIC/AveragingMethods/Dual/Dual.C | 62 +++++++------------ .../MPPIC/AveragingMethods/Dual/Dual.H | 21 ++----- .../MPPIC/AveragingMethods/Moment/Moment.C | 31 ++++++---- .../MPPIC/AveragingMethods/Moment/Moment.H | 8 +-- .../DampingModels/Relaxation/Relaxation.C | 4 +- .../IsotropyModels/Stochastic/Stochastic.C | 23 ++++--- .../MPPIC/PackingModels/Explicit/Explicit.C | 8 +-- .../MPPIC/PackingModels/Implicit/Implicit.C | 5 +- 13 files changed, 107 insertions(+), 124 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H index 1cebc01a2..bd26b231b 100644 --- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H +++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H @@ -180,10 +180,10 @@ Foam::MPPICParcel::TrackingData::updateAverages const scalar m = p.nParticle()*p.mass(); - volumeAverage_->add(p.position(), tetIs, p.nParticle()*p.volume()); - rhoAverage_->add(p.position(), tetIs, m*p.rho()); - uAverage_->add(p.position(), tetIs, m*p.U()); - massAverage_->add(p.position(), tetIs, m); + volumeAverage_->add(p.coordinates(), tetIs, p.nParticle()*p.volume()); + rhoAverage_->add(p.coordinates(), tetIs, m*p.rho()); + uAverage_->add(p.coordinates(), tetIs, m*p.U()); + massAverage_->add(p.coordinates(), tetIs, m); } volumeAverage_->average(); massAverage_->average(); @@ -196,11 +196,11 @@ Foam::MPPICParcel::TrackingData::updateAverages const typename CloudType::parcelType& p = iter(); const tetIndices tetIs = p.currentTetIndices(); - const vector u = uAverage_->interpolate(p.position(), tetIs); + const vector u = uAverage_->interpolate(p.coordinates(), tetIs); uSqrAverage_->add ( - p.position(), + p.coordinates(), tetIs, p.nParticle()*p.mass()*magSqr(p.U() - u) ); @@ -217,7 +217,7 @@ Foam::MPPICParcel::TrackingData::updateAverages weightAverage.add ( - p.position(), + p.coordinates(), tetIs, p.nParticle()*pow(p.volume(), 2.0/3.0) ); @@ -232,15 +232,15 @@ Foam::MPPICParcel::TrackingData::updateAverages const typename CloudType::parcelType& p = iter(); const tetIndices tetIs = p.currentTetIndices(); - const scalar a = volumeAverage_->interpolate(p.position(), tetIs); - const scalar r = radiusAverage_->interpolate(p.position(), tetIs); - const vector u = uAverage_->interpolate(p.position(), tetIs); + const scalar a = volumeAverage_->interpolate(p.coordinates(), tetIs); + const scalar r = radiusAverage_->interpolate(p.coordinates(), tetIs); + const vector u = uAverage_->interpolate(p.coordinates(), tetIs); const scalar f = 0.75*a/pow3(r)*sqr(0.5*p.d() + r)*mag(p.U() - u); - frequencyAverage_->add(p.position(), tetIs, p.nParticle()*f*f); + frequencyAverage_->add(p.coordinates(), tetIs, p.nParticle()*f*f); - weightAverage.add(p.position(), tetIs, p.nParticle()*f); + weightAverage.add(p.coordinates(), tetIs, p.nParticle()*f); } frequencyAverage_->average(weightAverage); } diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C index cb830f295..20c054e91 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C @@ -195,6 +195,16 @@ bool Foam::AveragingMethod::write() const dimensioned("zero", dimless, Zero) ); + // Barycentric coordinates of the tet vertices + const FixedList + tetCrds + ({ + barycentric(1, 0, 0, 0), + barycentric(0, 1, 0, 0), + barycentric(0, 0, 1, 0), + barycentric(0, 0, 0, 1) + }); + // tet-volume weighted sums forAll(mesh_.C(), celli) { @@ -207,18 +217,16 @@ bool Foam::AveragingMethod::write() const const triFace triIs = tetIs.faceTriIs(mesh_); const scalar v = tetIs.tet(mesh_).mag(); - cellValue[celli] += v*interpolate(mesh_.C()[celli], tetIs); - cellGrad[celli] += v*interpolateGrad(mesh_.C()[celli], tetIs); + cellValue[celli] += v*interpolate(tetCrds[0], tetIs); + cellGrad[celli] += v*interpolateGrad(tetCrds[0], tetIs); forAll(triIs, vertexI) { const label pointi = triIs[vertexI]; pointVolume[pointi] += v; - pointValue[pointi] += - v*interpolate(mesh_.points()[pointi], tetIs); - pointGrad[pointi] += - v*interpolateGrad(mesh_.points()[pointi], tetIs); + pointValue[pointi] += v*interpolate(tetCrds[vertexI], tetIs); + pointGrad[pointi] += v*interpolateGrad(tetCrds[vertexI], tetIs); } } } diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H index 42e8ce7cc..cdaea9863 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -38,6 +38,7 @@ SourceFiles #include "IOdictionary.H" #include "autoPtr.H" +#include "barycentric.H" #include "runTimeSelectionTables.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -134,7 +135,7 @@ public: //- Add point value to interpolation virtual void add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ) = 0; @@ -142,14 +143,14 @@ public: //- Interpolate virtual Type interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const = 0; //- Interpolate gradient virtual TypeGrad interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const = 0; diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C index 51dfe78ac..49258aeda 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,7 +91,7 @@ void Foam::AveragingMethods::Basic::updateGrad() template void Foam::AveragingMethods::Basic::add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ) @@ -103,7 +103,7 @@ void Foam::AveragingMethods::Basic::add template Type Foam::AveragingMethods::Basic::interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { @@ -115,7 +115,7 @@ template typename Foam::AveragingMethods::Basic::TypeGrad Foam::AveragingMethods::Basic::interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H index 7c13d475e..2adb016c4 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -125,7 +125,7 @@ public: //- Add point value to interpolation void add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ); @@ -133,14 +133,14 @@ public: //- Interpolate Type interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; //- Interpolate gradient TypeGrad interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C index 17fc6ba60..bbc44a98b 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C @@ -91,9 +91,7 @@ Foam::AveragingMethods::Dual::Dual volumeCell_(am.volumeCell_), volumeDual_(am.volumeDual_), dataCell_(FieldField::operator[](0)), - dataDual_(FieldField::operator[](1)), - tetVertices_(am.tetVertices_), - tetCoordinates_(am.tetCoordinates_) + dataDual_(FieldField::operator[](1)) {} @@ -106,24 +104,6 @@ Foam::AveragingMethods::Dual::~Dual() // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -template -void Foam::AveragingMethods::Dual::tetGeometry -( - const point position, - const tetIndices& tetIs -) const -{ - tetVertices_ = tetIs.faceTriIs(this->mesh_); - - tetIs.tet(this->mesh_).barycentric(position, tetCoordinates_); - - forAll(tetCoordinates_, i) - { - tetCoordinates_[i] = max(tetCoordinates_[i], scalar(0)); - } -} - - template void Foam::AveragingMethods::Dual::syncDualData() { @@ -141,22 +121,22 @@ void Foam::AveragingMethods::Dual::syncDualData() template void Foam::AveragingMethods::Dual::add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ) { - tetGeometry(position, tetIs); + const triFace triIs(tetIs.faceTriIs(this->mesh_)); dataCell_[tetIs.cell()] += - tetCoordinates_[0]*value + coordinates[0]*value / (0.25*volumeCell_[tetIs.cell()]); for(label i = 0; i < 3; i ++) { - dataDual_[tetVertices_[i]] += - tetCoordinates_[i+1]*value - / (0.25*volumeDual_[tetVertices_[i]]); + dataDual_[triIs[i]] += + coordinates[i+1]*value + / (0.25*volumeDual_[triIs[i]]); } } @@ -164,17 +144,17 @@ void Foam::AveragingMethods::Dual::add template Type Foam::AveragingMethods::Dual::interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { - tetGeometry(position, tetIs); + const triFace triIs(tetIs.faceTriIs(this->mesh_)); return - tetCoordinates_[0]*dataCell_[tetIs.cell()] - + tetCoordinates_[1]*dataDual_[tetVertices_[0]] - + tetCoordinates_[2]*dataDual_[tetVertices_[1]] - + tetCoordinates_[3]*dataDual_[tetVertices_[2]]; + coordinates[0]*dataCell_[tetIs.cell()] + + coordinates[1]*dataDual_[triIs[0]] + + coordinates[2]*dataDual_[triIs[1]] + + coordinates[3]*dataDual_[triIs[2]]; } @@ -182,11 +162,11 @@ template typename Foam::AveragingMethods::Dual::TypeGrad Foam::AveragingMethods::Dual::interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { - tetGeometry(position, tetIs); + const triFace triIs(tetIs.faceTriIs(this->mesh_)); const label celli(tetIs.cell()); @@ -196,9 +176,9 @@ Foam::AveragingMethods::Dual::interpolateGrad ( tensor ( - this->mesh_.points()[tetVertices_[0]] - this->mesh_.C()[celli], - this->mesh_.points()[tetVertices_[1]] - this->mesh_.C()[celli], - this->mesh_.points()[tetVertices_[2]] - this->mesh_.C()[celli] + this->mesh_.points()[triIs[0]] - this->mesh_.C()[celli], + this->mesh_.points()[triIs[1]] - this->mesh_.C()[celli], + this->mesh_.points()[triIs[2]] - this->mesh_.C()[celli] ) ) ); @@ -207,9 +187,9 @@ Foam::AveragingMethods::Dual::interpolateGrad const TypeGrad S ( - dataDual_[tetVertices_[0]], - dataDual_[tetVertices_[1]], - dataDual_[tetVertices_[2]] + dataDual_[triIs[0]], + dataDual_[triIs[1]], + dataDual_[triIs[2]] ); const Type s(dataCell_[celli]); diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H index 94c2a3325..903c2a742 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -90,12 +90,6 @@ private: //- Data on the points Field& dataDual_; - //- Tet vertex labels - mutable FixedList tetVertices_; - - //- Tet barycentric coordinates - mutable FixedList tetCoordinates_; - //- Private static member functions @@ -105,13 +99,6 @@ private: //- Private member functions - //- Calculate indices and barycentric coordinates within a tetrahedron - void tetGeometry - ( - const point position, - const tetIndices& tetIs - ) const; - //- Sync point data over processor boundaries void syncDualData(); @@ -154,7 +141,7 @@ public: //- Add point value to interpolation void add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ); @@ -162,14 +149,14 @@ public: //- Interpolate Type interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; //- Interpolate gradient TypeGrad interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C index 40f048d31..94ef60aa8 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C @@ -128,21 +128,22 @@ void Foam::AveragingMethods::Moment::updateGrad() template void Foam::AveragingMethods::Moment::add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ) { const label celli = tetIs.cell(); + const triFace triIs = tetIs.faceTriIs(this->mesh_); + + const point delta = + (coordinates[0] - 1)*this->mesh_.C()[celli] + + coordinates[1]*this->mesh_.points()[triIs[0]] + + coordinates[2]*this->mesh_.points()[triIs[1]] + + coordinates[3]*this->mesh_.points()[triIs[2]]; const Type v = value/this->mesh_.V()[celli]; - const TypeGrad dv = - transform_[celli] - & ( - v - * (position - this->mesh_.C()[celli]) - / scale_[celli] - ); + const TypeGrad dv = transform_[celli] & (v*delta/scale_[celli]); data_[celli] += v; dataX_[celli] += v + dv.x(); @@ -154,11 +155,18 @@ void Foam::AveragingMethods::Moment::add template Type Foam::AveragingMethods::Moment::interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { const label celli = tetIs.cell(); + const triFace triIs = tetIs.faceTriIs(this->mesh_); + + const point delta = + (coordinates[0] - 1)*this->mesh_.C()[celli] + + coordinates[1]*this->mesh_.points()[triIs[0]] + + coordinates[2]*this->mesh_.points()[triIs[1]] + + coordinates[3]*this->mesh_.points()[triIs[2]]; return data_[celli] @@ -169,8 +177,7 @@ Type Foam::AveragingMethods::Moment::interpolate dataY_[celli] - data_[celli], dataZ_[celli] - data_[celli] ) - & (position - this->mesh_.C()[celli]) - / scale_[celli] + & delta/scale_[celli] ); } @@ -179,7 +186,7 @@ template typename Foam::AveragingMethods::Moment::TypeGrad Foam::AveragingMethods::Moment::interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const { diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H index 8bda35212..6f142ef08 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H +++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -138,7 +138,7 @@ public: //- Add point value to interpolation void add ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs, const Type& value ); @@ -146,14 +146,14 @@ public: //- Interpolate Type interpolate ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; //- Interpolate gradient TypeGrad interpolateGrad ( - const point position, + const barycentric& coordinates, const tetIndices& tetIs ) const; diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C index 811c918ed..a414d5d56 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C @@ -142,9 +142,9 @@ Foam::vector Foam::DampingModels::Relaxation::velocityCorrection const tetIndices tetIs(p.currentTetIndices()); const scalar x = - deltaT*oneByTimeScaleAverage_->interpolate(p.position(), tetIs); + deltaT*oneByTimeScaleAverage_->interpolate(p.coordinates(), tetIs); - const vector u = uAverage_->interpolate(p.position(), tetIs); + const vector u = uAverage_->interpolate(p.coordinates(), tetIs); return (u - p.U())*x/(x + 2.0); } diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C index 0b047b294..9cf935e22 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C @@ -168,15 +168,15 @@ void Foam::IsotropyModels::Stochastic::calculate() typename CloudType::parcelType& p = iter(); const tetIndices tetIs(p.currentTetIndices()); - const scalar x = exponentAverage.interpolate(p.position(), tetIs); + const scalar x = exponentAverage.interpolate(p.coordinates(), tetIs); if (x < rndGen.sample01()) { const vector r(sampleGauss(), sampleGauss(), sampleGauss()); - const vector u = uAverage.interpolate(p.position(), tetIs); + const vector u = uAverage.interpolate(p.coordinates(), tetIs); const scalar uRms = - sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0)); + sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0)); p.U() = u + r*uRms*oneBySqrtThree; } @@ -202,7 +202,7 @@ void Foam::IsotropyModels::Stochastic::calculate() { typename CloudType::parcelType& p = iter(); const tetIndices tetIs(p.currentTetIndices()); - uTildeAverage.add(p.position(), tetIs, p.nParticle()*p.mass()*p.U()); + uTildeAverage.add(p.coordinates(), tetIs, p.nParticle()*p.mass()*p.U()); } uTildeAverage.average(massAverage); @@ -225,10 +225,10 @@ void Foam::IsotropyModels::Stochastic::calculate() { typename CloudType::parcelType& p = iter(); const tetIndices tetIs(p.currentTetIndices()); - const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs); + const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs); uTildeSqrAverage.add ( - p.position(), + p.coordinates(), tetIs, p.nParticle()*p.mass()*magSqr(p.U() - uTilde) ); @@ -241,13 +241,16 @@ void Foam::IsotropyModels::Stochastic::calculate() typename CloudType::parcelType& p = iter(); const tetIndices tetIs(p.currentTetIndices()); - const vector u = uAverage.interpolate(p.position(), tetIs); + const vector u = uAverage.interpolate(p.coordinates(), tetIs); const scalar uRms = - sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0)); + sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0)); - const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs); + const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs); const scalar uTildeRms = - sqrt(max(uTildeSqrAverage.interpolate(p.position(), tetIs), 0.0)); + sqrt + ( + max(uTildeSqrAverage.interpolate(p.coordinates(), tetIs), 0.0) + ); p.U() = u + (p.U() - uTilde)*uRms/max(uTildeRms, SMALL); } diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C index fda103663..9667fcb03 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C @@ -147,15 +147,15 @@ Foam::vector Foam::PackingModels::Explicit::velocityCorrection // interpolated quantities const scalar alpha = - this->volumeAverage_->interpolate(p.position(), tetIs); + this->volumeAverage_->interpolate(p.coordinates(), tetIs); const vector alphaGrad = - this->volumeAverage_->interpolateGrad(p.position(), tetIs); + this->volumeAverage_->interpolateGrad(p.coordinates(), tetIs); const vector uMean = - this->uAverage_->interpolate(p.position(), tetIs); + this->uAverage_->interpolate(p.coordinates(), tetIs); // stress gradient const vector tauGrad = - stressAverage_->interpolateGrad(p.position(), tetIs); + stressAverage_->interpolateGrad(p.coordinates(), tetIs); // parcel relative velocity const vector uRelative = p.U() - uMean; diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C index 9f82d5318..0f65a817d 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C @@ -339,9 +339,6 @@ Foam::vector Foam::PackingModels::Implicit::velocityCorrection // containing tetrahedron and parcel coordinates within const label celli = p.cell(); const label facei = p.tetFace(); - const tetIndices tetIs(p.currentTetIndices()); - FixedList tetCoordinates; - tetIs.tet(mesh).barycentric(p.position(), tetCoordinates); // cell velocity const vector U = uCorrect_()[celli]; @@ -368,7 +365,7 @@ Foam::vector Foam::PackingModels::Implicit::velocityCorrection } // interpolant equal to 1 at the cell centre and 0 at the face - const scalar t = tetCoordinates[0]; + const scalar t = p.coordinates()[0]; // the normal component of the velocity correction is interpolated linearly // the tangential component is equal to that at the cell centre