From 6b316ba3a92cacfafe1dd7cda6756ff7a8bd57d2 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 | 20 ++++++------------- .../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(+), 50 deletions(-) diff --git a/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H b/src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H index f47e52f9c9..ca60aff4f4 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 ce6c0df882..d2e61b8e43 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 d0b11f8ff8..d95ad7dcb3 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 7a4302f7ae..428054ed11 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 e8e5bd3fa5..a547a72e60 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,21 +79,13 @@ inline Type Foam::interpolationCellPoint::interpolate } } - FixedList weights; - - tetIs.tet(this->pMesh_).barycentric(position, weights); - - // 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 a2bd9c7518..2cfa869fc3 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 803c7a440a..96f90becdc 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 0a0b755ea3..49888b5140 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -48,7 +48,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()) { @@ -62,9 +62,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 08be00dae3..3300c2c3e0 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 @@ -211,7 +211,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 90f52a9bbc..4a906cb642 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 c9c7fbc0b5..1268179b73 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 8426d37140..a268a07b09 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 e0d6518fba..0059a6627f 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 ec189b4a14..d21041a4ec 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_);