Commit 7a02a507 authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

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.
parent e9fb8b85
......@@ -180,10 +180,10 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::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<ParcelType>::TrackingData<CloudType>::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<ParcelType>::TrackingData<CloudType>::updateAverages
weightAverage.add
(
p.position(),
p.coordinates(),
tetIs,
p.nParticle()*pow(p.volume(), 2.0/3.0)
);
......@@ -232,15 +232,15 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::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);
}
......
......@@ -194,6 +194,16 @@ bool Foam::AveragingMethod<Type>::write() const
dimensioned<TypeGrad>("zero", dimless, Zero)
);
// Barycentric coordinates of the tet vertices
const FixedList<barycentric, 4>
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)
{
......@@ -206,18 +216,16 @@ bool Foam::AveragingMethod<Type>::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);
}
}
}
......
......@@ -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
......@@ -41,6 +41,7 @@ SourceFiles
#include "IOdictionary.H"
#include "autoPtr.H"
#include "barycentric.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -137,7 +138,7 @@ public:
//- Add point value to interpolation
virtual void add
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs,
const Type& value
) = 0;
......@@ -145,14 +146,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;
......
......@@ -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<Type>::updateGrad()
template<class Type>
void Foam::AveragingMethods::Basic<Type>::add
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs,
const Type& value
)
......@@ -103,7 +103,7 @@ void Foam::AveragingMethods::Basic<Type>::add
template<class Type>
Type Foam::AveragingMethods::Basic<Type>::interpolate
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs
) const
{
......@@ -115,7 +115,7 @@ template<class Type>
typename Foam::AveragingMethods::Basic<Type>::TypeGrad
Foam::AveragingMethods::Basic<Type>::interpolateGrad
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs
) const
{
......
......@@ -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
......@@ -128,7 +128,7 @@ public:
//- Add point value to interpolation
void add
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs,
const Type& value
);
......@@ -136,14 +136,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;
......
......@@ -91,9 +91,7 @@ Foam::AveragingMethods::Dual<Type>::Dual
volumeCell_(am.volumeCell_),
volumeDual_(am.volumeDual_),
dataCell_(FieldField<Field, Type>::operator[](0)),
dataDual_(FieldField<Field, Type>::operator[](1)),
tetVertices_(am.tetVertices_),
tetCoordinates_(am.tetCoordinates_)
dataDual_(FieldField<Field, Type>::operator[](1))
{}
......@@ -106,24 +104,6 @@ Foam::AveragingMethods::Dual<Type>::~Dual()
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
void Foam::AveragingMethods::Dual<Type>::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<class Type>
void Foam::AveragingMethods::Dual<Type>::syncDualData()
{
......@@ -141,22 +121,22 @@ void Foam::AveragingMethods::Dual<Type>::syncDualData()
template<class Type>
void Foam::AveragingMethods::Dual<Type>::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<Type>::add
template<class Type>
Type Foam::AveragingMethods::Dual<Type>::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<class Type>
typename Foam::AveragingMethods::Dual<Type>::TypeGrad
Foam::AveragingMethods::Dual<Type>::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<Type>::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<Type>::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]);
......
......@@ -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
......@@ -93,12 +93,6 @@ private:
//- Data on the points
Field<Type>& dataDual_;
//- Tet vertex labels
mutable FixedList<label, 3> tetVertices_;
//- Tet barycentric coordinates
mutable FixedList<scalar, 4> tetCoordinates_;
//- Private static member functions
......@@ -108,13 +102,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();
......@@ -157,7 +144,7 @@ public:
//- Add point value to interpolation
void add
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs,
const Type& value
);
......@@ -165,14 +152,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;
......
......@@ -128,21 +128,22 @@ void Foam::AveragingMethods::Moment<Type>::updateGrad()
template<class Type>
void Foam::AveragingMethods::Moment<Type>::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<Type>::add
template<class Type>
Type Foam::AveragingMethods::Moment<Type>::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<Type>::interpolate
dataY_[celli] - data_[celli],
dataZ_[celli] - data_[celli]
)
& (position - this->mesh_.C()[celli])
/ scale_[celli]
& delta/scale_[celli]
);
}
......@@ -179,7 +186,7 @@ template<class Type>
typename Foam::AveragingMethods::Moment<Type>::TypeGrad
Foam::AveragingMethods::Moment<Type>::interpolateGrad
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs
) const
{
......
......@@ -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
......@@ -141,7 +141,7 @@ public:
//- Add point value to interpolation
void add
(
const point position,
const barycentric& coordinates,
const tetIndices& tetIs,
const Type& value
);
......@@ -149,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;
......
......@@ -142,9 +142,9 @@ Foam::vector Foam::DampingModels::Relaxation<CloudType>::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);
}
......
......@@ -168,15 +168,15 @@ void Foam::IsotropyModels::Stochastic<CloudType>::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<scalar>())
{
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<CloudType>::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<CloudType>::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<CloudType>::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)