From c1cb13307b68e357e3a75c78cdca2caaa0b1998d Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 6 Jun 2017 16:07:56 +0100 Subject: [PATCH] interpolation: Optimise by using particle local coordinates This change changes the point-tetIndices-face interpolation function method to take barycentric-tetIndices-face arguments instead. This function is, at present, only used for interpolating Eulerian data to Lagrangian particles. This change prevents an inefficiency in cellPointInterpolation whereby the position of the particle is calculated from it's barycentric coordinates, before immediately being converted back to barycentric coordinates to perform the interpolation. --- .../interpolation/interpolation.H | 14 ++++++++++---- .../interpolationCell/interpolationCell.H | 15 ++++++++++++++- .../interpolationCellPatchConstrained.H | 15 ++++++++++++++- .../interpolationCellPoint.H | 6 +++--- .../interpolationCellPointI.H | 19 ++++++------------- .../interpolationCellPointWallModified.H | 6 +++--- .../interpolationCellPointWallModifiedI.H | 4 ++-- .../KinematicParcel/KinematicParcel.C | 6 +++--- .../Templates/ReactingParcel/ReactingParcel.C | 4 ++-- .../Templates/ThermoParcel/ThermoParcel.C | 12 ++++++------ .../ParticleForces/Lift/LiftForce/LiftForce.C | 4 ++-- .../Paramagnetic/ParamagneticForce.C | 4 ++-- .../PressureGradient/PressureGradientForce.C | 4 ++-- src/lagrangian/solidParticle/solidParticle.C | 9 ++++----- 14 files changed, 73 insertions(+), 49 deletions(-) diff --git a/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H b/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H index f47e52f9c..ca60aff4f 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H +++ b/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -135,18 +135,24 @@ public: const label facei = -1 ) const = 0; - //- Interpolate field to the given point in the tetrahedron + //- Interpolate field to the given coordinates in the tetrahedron // defined by the given indices. Calls interpolate function // above here execpt where overridden by derived // interpolation types. virtual Type interpolate ( - const vector& position, + const barycentric& coordinates, const tetIndices& tetIs, const label facei = -1 ) const { - return interpolate(position, tetIs.cell(), facei); + return + interpolate + ( + tetIs.tet(pMesh_).barycentricToPoint(coordinates), + tetIs.cell(), + facei + ); } }; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H b/src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H index ce6c0df88..d2e61b8e4 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -75,6 +75,19 @@ public: const label celli, const label facei = -1 ) const; + + //- Interpolate field to the given coordinates in the tetrahedron + // defined by the given indices. This is an optimisation which skips + // calculating the position, as cell interpolation doesn't need it. + inline Type interpolate + ( + const barycentric& coordinates, + const tetIndices& tetIs, + const label facei = -1 + ) const + { + return interpolate(vector::zero, tetIs.cell(), facei); + } }; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H index d0b11f8ff..d95ad7dcb 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,6 +77,19 @@ public: const label celli, const label facei = -1 ) const; + + //- Interpolate field to the given coordinates in the tetrahedron + // defined by the given indices. This is an optimisation which skips + // calculating the position, as cell interpolation doesn't need it. + inline Type interpolate + ( + const barycentric& coordinates, + const tetIndices& tetIs, + const label facei = -1 + ) const + { + return interpolate(vector::zero, tetIs.cell(), facei); + } }; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H index 7a4302f7a..428054ed1 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -86,11 +86,11 @@ public: const label facei = -1 ) const; - //- Interpolate field to the given point in the tetrahedron + //- Interpolate field to the given coordinates in the tetrahedron // defined by the given indices. inline Type interpolate ( - const vector& position, + const barycentric& coordinates, const tetIndices& tetIs, const label facei = -1 ) const; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H index 13942628b..e46e76b91 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H @@ -58,7 +58,7 @@ inline Type Foam::interpolationCellPoint::interpolate template inline Type Foam::interpolationCellPoint::interpolate ( - const vector& position, + const barycentric& coordinates, const tetIndices& tetIs, const label facei ) const @@ -79,20 +79,13 @@ inline Type Foam::interpolationCellPoint::interpolate } } - const barycentric weights = - tetIs.tet(this->pMesh_).pointToBarycentric(position); - - // Order of weights is the same as that of the vertices of the tet, i.e. - // cellCentre, faceBasePt, facePtA, facePtB. - - Type t = this->psi_[tetIs.cell()]*weights[0]; - const triFace triIs = tetIs.faceTriIs(this->pMesh_); - t += psip_[triIs[0]]*weights[1]; - t += psip_[triIs[1]]*weights[2]; - t += psip_[triIs[2]]*weights[3]; - return t; + return + this->psi_[tetIs.cell()]*coordinates[0] + + psip_[triIs[0]]*coordinates[1] + + psip_[triIs[1]]*coordinates[2] + + psip_[triIs[2]]*coordinates[3]; } diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H index a2bd9c751..2cfa869fc 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,11 +78,11 @@ public: const label facei = -1 ) const; - //- Interpolate field to the given point in the tetrahedron + //- Interpolate field to the given coordinates in the tetrahedron // defined by the given indices. inline Type interpolate ( - const vector& position, + const barycentric& coordinates, const tetIndices& tetIs, const label facei = -1 ) const; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H index ad72ef453..e55a68f4c 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H @@ -67,7 +67,7 @@ inline Type Foam::interpolationCellPointWallModified::interpolate template inline Type Foam::interpolationCellPointWallModified::interpolate ( - const vector& position, + const barycentric& coordinates, const tetIndices& tetIs, const label facei ) const @@ -101,7 +101,7 @@ inline Type Foam::interpolationCellPointWallModified::interpolate return interpolationCellPoint::interpolate ( - position, + coordinates, tetIs, facei ); diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 26ae482e7..49bb96921 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -47,7 +47,7 @@ void Foam::KinematicParcel::setCellValues { tetIndices tetIs = this->currentTetIndices(); - rhoc_ = td.rhoInterp().interpolate(this->position(), tetIs); + rhoc_ = td.rhoInterp().interpolate(this->coordinates(), tetIs); if (rhoc_ < td.cloud().constProps().rhoMin()) { @@ -61,9 +61,9 @@ void Foam::KinematicParcel::setCellValues rhoc_ = td.cloud().constProps().rhoMin(); } - Uc_ = td.UInterp().interpolate(this->position(), tetIs); + Uc_ = td.UInterp().interpolate(this->coordinates(), tetIs); - muc_ = td.muInterp().interpolate(this->position(), tetIs); + muc_ = td.muInterp().interpolate(this->coordinates(), tetIs); // Apply dispersion components to carrier phase velocity Uc_ = td.cloud().dispersion().update diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index a861982cd..bad44e1bc 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -213,7 +213,7 @@ void Foam::ReactingParcel::setCellValues pc_ = td.pInterp().interpolate ( - this->position(), + this->coordinates(), this->currentTetIndices() ); diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C index 103c97c34..46454ec7a 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,9 +43,9 @@ void Foam::ThermoParcel::setCellValues tetIndices tetIs = this->currentTetIndices(); - Cpc_ = td.CpInterp().interpolate(this->position(), tetIs); + Cpc_ = td.CpInterp().interpolate(this->coordinates(), tetIs); - Tc_ = td.TInterp().interpolate(this->position(), tetIs); + Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs); if (Tc_ < td.cloud().constProps().TMin()) { @@ -124,8 +124,8 @@ void Foam::ThermoParcel::calcSurfaceValues rhos = this->rhoc_*TRatio; tetIndices tetIs = this->currentTetIndices(); - mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio; - kappas = td.kappaInterp().interpolate(this->position(), tetIs)/TRatio; + mus = td.muInterp().interpolate(this->coordinates(), tetIs)/TRatio; + kappas = td.kappaInterp().interpolate(this->coordinates(), tetIs)/TRatio; Pr = Cpc_*mus/kappas; Pr = max(ROOTVSMALL, Pr); @@ -286,7 +286,7 @@ Foam::scalar Foam::ThermoParcel::calcHeatTransfer if (td.cloud().radiation()) { tetIndices tetIs = this->currentTetIndices(); - const scalar Gc = td.GInterp().interpolate(this->position(), tetIs); + const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs); const scalar sigma = physicoChemical::sigma.value(); const scalar epsilon = td.cloud().constProps().epsilon0(); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C index c9c7fbc0b..1268179b7 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -137,7 +137,7 @@ Foam::forceSuSp Foam::LiftForce::calcCoupled forceSuSp value(Zero, 0.0); vector curlUc = - curlUcInterp().interpolate(p.position(), p.currentTetIndices()); + curlUcInterp().interpolate(p.coordinates(), p.currentTetIndices()); scalar Cl = this->Cl(p, curlUc, Re, muc); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C index 8426d3714..a268a07b0 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -110,7 +110,7 @@ Foam::forceSuSp Foam::ParamagneticForce::calcNonCoupled value.Su()= mass*3.0*constant::electromagnetic::mu0.value()/p.rho() *magneticSusceptibility_/(magneticSusceptibility_ + 3) - *HdotGradHInterp.interpolate(p.position(), p.currentTetIndices()); + *HdotGradHInterp.interpolate(p.coordinates(), p.currentTetIndices()); return value; } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C index e0d6518fb..0059a6627 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -128,7 +128,7 @@ Foam::forceSuSp Foam::PressureGradientForce::calcCoupled forceSuSp value(Zero, 0.0); vector DUcDt = - DUcDtInterp().interpolate(p.position(), p.currentTetIndices()); + DUcDtInterp().interpolate(p.coordinates(), p.currentTetIndices()); value.Su() = mass*p.rhoc()/p.rho()*DUcDt; diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C index ec189b4a1..d21041a4e 100644 --- a/src/lagrangian/solidParticle/solidParticle.C +++ b/src/lagrangian/solidParticle/solidParticle.C @@ -55,7 +55,6 @@ bool Foam::solidParticle::move } - const label celli = cell(); const scalar sfrac = stepFraction(); const scalar f = 1 - stepFraction(); @@ -63,10 +62,10 @@ bool Foam::solidParticle::move const scalar dt = (stepFraction() - sfrac)*trackTime; - cellPointWeight cpw(mesh(), position(), celli, face()); - scalar rhoc = td.rhoInterp().interpolate(cpw); - vector Uc = td.UInterp().interpolate(cpw); - scalar nuc = td.nuInterp().interpolate(cpw); + const tetIndices tetIs = this->currentTetIndices(); + scalar rhoc = td.rhoInterp().interpolate(this->coordinates(), tetIs); + vector Uc = td.UInterp().interpolate(this->coordinates(), tetIs); + scalar nuc = td.nuInterp().interpolate(this->coordinates(), tetIs); scalar rhop = td.cloud().rhop(); scalar magUr = mag(Uc - U_);