From 19cb3b0454a4d837cbff950dd868ce008b48b71a Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Mon, 21 Nov 2016 09:44:00 +0000 Subject: [PATCH 1/7] STYLE: porosityModel - minor cleanup of bounds output --- .../porosityModel/porosityModel.C | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C index 77d28f589c..eb5550222e 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C +++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C @@ -128,30 +128,30 @@ Foam::porosityModel::porosityModel const pointField& points = mesh_.points(); const cellList& cells = mesh_.cells(); const faceList& faces = mesh_.faces(); - DynamicList localPoints; forAll(cellZoneIDs_, zoneI) { const cellZone& cZone = mesh_.cellZones()[cellZoneIDs_[zoneI]]; - localPoints.setCapacity(10*cells.size()); + point bbMin = point::max; + point bbMax = point::min; forAll(cZone, i) { const label cellI = cZone[i]; - const cell& c = mesh_.cells()[cellI]; + const cell& c = cells[cellI]; const pointField cellPoints(c.points(faces, points)); forAll(cellPoints, pointI) { - const point& pt = cellPoints[pointI]; - localPoints.append(coordSys_.localPosition(pt)); + const point pt = coordSys_.localPosition(cellPoints[pointI]); + bbMin = min(bbMin, pt); + bbMax = max(bbMax, pt); } } - boundBox bb(localPoints, true); + reduce(bbMin, minOp()); + reduce(bbMax, maxOp()); - Info<< " local bounds: " << bb << endl; - - localPoints.clear(); + Info<< " local bounds: " << (bbMax - bbMin) << nl << endl; } } From 5325af46265b4fe3bf33f991402a114506697a49 Mon Sep 17 00:00:00 2001 From: sergio Date: Mon, 21 Nov 2016 11:54:31 -0800 Subject: [PATCH 2/7] BUG: Following issue 284 for particle switching processors flag --- src/lagrangian/basic/Cloud/Cloud.C | 7 ++++++- .../parcels/Templates/KinematicParcel/KinematicParcel.C | 4 ++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C index 0ea3551bb1..a0ebd79048 100644 --- a/src/lagrangian/basic/Cloud/Cloud.C +++ b/src/lagrangian/basic/Cloud/Cloud.C @@ -261,7 +261,12 @@ void Foam::Cloud::move(TrackData& td, const scalar trackTime) { // If we are running in parallel and the particle is on a // boundary face - if (Pstream::parRun() && p.face() >= pMesh().nInternalFaces()) + if + ( + Pstream::parRun() + && td.switchProcessor + && p.face() >= pMesh().nInternalFaces() + ) { label patchI = pbm.whichPatch(p.face()); diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 62f5d50162..7fdf811542 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -345,6 +345,8 @@ bool Foam::KinematicParcel::move p.calc(td, dt, cellI); } + p.age() += dt; + if (p.onBoundary() && td.keepParticle) { if (isA(pbMesh[p.patch(p.face())])) @@ -353,8 +355,6 @@ bool Foam::KinematicParcel::move } } - p.age() += dt; - td.cloud().functions().postMove(p, cellI, dt, start, td.keepParticle); } From 05df4c71b68ad79565de2f84c909d900892a57ee Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Tue, 29 Nov 2016 09:10:37 +0000 Subject: [PATCH 3/7] ENH: Lagrangian models - added headers to enable derived libraries to compile --- .../parcels/Templates/CollidingParcel/CollidingParcel.H | 1 + .../Kinematic/ParticleForces/ParticleForce/ParticleForce.H | 5 +++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H index 338c3186d9..4354d08297 100644 --- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H @@ -44,6 +44,7 @@ SourceFiles #include "CollisionRecordList.H" #include "labelFieldIOField.H" #include "vectorFieldIOField.H" +#include "demandDrivenEntry.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H index ba0ef93db4..872d7f87f6 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H @@ -43,6 +43,7 @@ SourceFiles #include "dictionary.H" #include "forceSuSp.H" #include "fvMesh.H" +#include "runTimeSelectionTables.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -64,7 +65,7 @@ class ParticleForce //- Reference to the mesh database const fvMesh& mesh_; - //- Force coefficients dictaionary + //- Force coefficients dictionary const dictionary coeffs_; @@ -138,7 +139,7 @@ public: //- Return const access to the cloud owner inline const CloudType& owner() const; - //- Return refernce to the cloud owner + //- Return reference to the cloud owner inline CloudType& owner(); //- Return the mesh database From e73d2c5d29c5ebeca85b4b159beb99408c97b6f2 Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Tue, 6 Dec 2016 16:19:01 +0000 Subject: [PATCH 4/7] BUG: surfaceNoise - corrected fftWriteInterval operation --- .../noise/noiseModels/surfaceNoise/surfaceNoise.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C index 92a3df57fd..7d1a6e3e6c 100644 --- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C +++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C @@ -562,7 +562,7 @@ void surfaceNoise::calculate() forAll(surfPrmsf, i) { - label freqI = i*fftWriteInterval_; + label freqI = (i + 1)*fftWriteInterval_ - 1; fOut[i] = freq1[freqI]; const word& fName = inputFileName_.name(true); const word gName = "fft"; From 11479c51d776489cef47cae041a9d9102364a93a Mon Sep 17 00:00:00 2001 From: sergio Date: Wed, 7 Dec 2016 12:57:32 -0800 Subject: [PATCH 5/7] ENH: Changes in handling topological changes in VOF solvers 1) Using divU instead of fvc::absolute(phi,U) in TEqn as the latter uses latest time meshPhi which is inconsistent 2) Adding fvc::interpolate(U) when topo changes 3) in pEq for compressible dgdt is updated using the latest rho1 and rho2 after compressible effects are considered --- .../multiphase/compressibleInterFoam/TEqn.H | 6 +- .../compressibleInterDyMFoam.C | 104 ++++++++++++------ .../compressibleInterDyMFoam/pEqn.H | 20 ++-- .../compressibleInterFoam.C | 1 + .../multiphase/compressibleInterFoam/pEqn.H | 20 ++-- .../interFoam/interDyMFoam/interDyMFoam.C | 7 +- 6 files changed, 98 insertions(+), 60 deletions(-) diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H index 13e1feb546..970e97157f 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H @@ -5,8 +5,8 @@ + fvm::div(rhoPhi, T) - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T) + ( - fvc::div(fvc::absolute(phi, U), p) - + fvc::ddt(rho, K) + fvc::div(rhoPhi, K) + divU*p + + fvc::ddt(rho, K) + fvc::div(rhoPhi, K) ) *( alpha1/mixture.thermo1().Cv() @@ -20,4 +20,4 @@ mixture.correct(); Info<< "min(T) " << min(T).value() << endl; -} +} \ No newline at end of file diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C index 7cbbd26b14..daef88206b 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C @@ -70,6 +70,21 @@ int main(int argc, char *argv[]) turbulence->validate(); + + volScalarField rAU + ( + IOobject + ( + "rAU", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("rAU", dimTime/rho.dimensions(), 1.0) + ); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -77,45 +92,17 @@ int main(int argc, char *argv[]) { #include "readControls.H" - { - // Store divU from the previous mesh so that it can be mapped - // and used in correctPhi to ensure the corrected phi has the - // same divergence - volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U))); + // Store divU from the previous mesh so that it can be mapped + // and used in correctPhi to ensure the corrected phi has the + // same divergence + volScalarField divU("divU0",fvc::div(fvc::absolute(phi, U))); - #include "CourantNo.H" - #include "setDeltaT.H" + #include "CourantNo.H" + #include "setDeltaT.H" - runTime++; + runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; - - scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); - - // Do any mesh changes - mesh.update(); - - if (mesh.changing()) - { - Info<< "Execution time for mesh.update() = " - << runTime.elapsedCpuTime() - timeBeforeMeshUpdate - << " s" << endl; - - gh = (g & mesh.C()) - ghRef; - ghf = (g & mesh.Cf()) - ghRef; - } - - if (mesh.changing() && correctPhi) - { - // Calculate absolute flux from the mapped surface velocity - phi = mesh.Sf() & Uf; - - #include "correctPhi.H" - - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); - } - } + Info<< "Time = " << runTime.timeName() << nl << endl; if (mesh.changing() && checkMeshCourantNo) { @@ -127,6 +114,46 @@ int main(int argc, char *argv[]) // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { + if (pimple.firstIter()) + { + scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); + + mesh.update(); + + if (mesh.changing()) + { + Info<< "Execution time for mesh.update() = " + << runTime.elapsedCpuTime() - timeBeforeMeshUpdate + << " s" << endl; + + gh = (g & mesh.C()) - ghRef; + ghf = (g & mesh.Cf()) - ghRef; + } + + if ((correctPhi && mesh.changing()) || mesh.topoChanging()) + { + // Calculate absolute flux from the mapped surface velocity + // SAF: temporary fix until mapped Uf is assessed + Uf = fvc::interpolate(U); + + phi = mesh.Sf() & Uf; + + #include "correctPhi.H" + + fvc::makeRelative(phi, U); + + mixture.correct(); + + mesh.topoChanging(false); + } + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + } + + #include "alphaEqnsSubCycle.H" // correct interface on first PIMPLE corrector @@ -145,6 +172,11 @@ int main(int argc, char *argv[]) { #include "pEqn.H" } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } } rho = alpha1*rho1 + alpha2*rho2; diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H index 859fc9cc4e..34e39fadde 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H @@ -87,15 +87,6 @@ if (pimple.finalNonOrthogonalIter()) { - p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); - p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; - - dgdt = - ( - pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 - - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 - ); - phi = phiHbyA + p_rghEqnIncomp.flux(); U = HbyA @@ -116,8 +107,19 @@ rho = alpha1*rho1 + alpha2*rho2; + p = max(p_rgh + rho*gh, pMin); + p_rgh = p - rho*gh; + + dgdt = + ( + pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 + - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 + ); + K = 0.5*magSqr(U); Info<< "max(U) " << max(mag(U)).value() << endl; Info<< "min(p_rgh) " << min(p_rgh).value() << endl; } + + diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index 2a113880de..e4d413964d 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C @@ -93,6 +93,7 @@ int main(int argc, char *argv[]) solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" + volScalarField divU(fvc::div(fvc::absolute(phi, U))); #include "TEqn.H" // --- Pressure corrector loop diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index a4a3fb2e29..f24b72eb68 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -84,15 +84,6 @@ if (pimple.finalNonOrthogonalIter()) { - p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); - p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; - - dgdt = - ( - pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 - - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 - ); - phi = phiHbyA + p_rghEqnIncomp.flux(); U = HbyA @@ -101,14 +92,21 @@ } } - // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); - // Update densities from change in p_rgh rho1 += psi1*(p_rgh - p_rgh_0); rho2 += psi2*(p_rgh - p_rgh_0); rho = alpha1*rho1 + alpha2*rho2; + p = max(p_rgh + rho*gh, pMin); + p_rgh = p - rho*gh; + + dgdt = + ( + pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2 + - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1 + ); + K = 0.5*magSqr(U); Info<< "max(U) " << max(mag(U)).value() << endl; diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 340f2c948e..401ece3fde 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -133,9 +133,12 @@ int main(int argc, char *argv[]) ghf = (g & mesh.Cf()) - ghRef; } - if (mesh.changing() && correctPhi) + if ((mesh.changing() && correctPhi) || mesh.topoChanging()) { // Calculate absolute flux from the mapped surface velocity + // SAF: temporary fix until mapped Uf is assessed + Uf = fvc::interpolate(U); + phi = mesh.Sf() & Uf; #include "correctPhi.H" @@ -144,6 +147,8 @@ int main(int argc, char *argv[]) fvc::makeRelative(phi, U); mixture.correct(); + + mesh.topoChanging(false); } if (mesh.changing() && checkMeshCourantNo) From dfb5fae5de8a67e55c769c403ffffc285c2ab852 Mon Sep 17 00:00:00 2001 From: sergio Date: Tue, 13 Dec 2016 09:44:25 -0800 Subject: [PATCH 6/7] BUG: Fixing oldTime Tw and Reference Value of humidityTemperatureCoupledMixed BC --- ...TemperatureCoupledMixedFvPatchScalarField.C | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C b/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C index 4e32aedf59..3d5f1e8e5c 100644 --- a/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -449,7 +449,14 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs() scalarField nbrIntFld(nbrField.patchInternalField()); mpp.distribute(nbrIntFld); + scalarField& Tp = *this; + + const volScalarField& T = + static_cast(dimensionedInternalField()); + + const fvPatchField& TpOld = T.boundaryField()[patch().index()]; + scalarField Tin(patchInternalField()); const scalarField K(this->kappa(*this)); @@ -463,6 +470,7 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs() mpp.distribute(nrbDeltaCoeffs); scalarField KDeltaNbr(nbrK*nrbDeltaCoeffs); + mpp.distribute(KDeltaNbr); myKDelta_ = K*patch().deltaCoeffs(); @@ -655,15 +663,11 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs() } } - scalarField myKDeltaNbr(patch().size(), 0.0); scalarField mpCpTpNbr(patch().size(), 0.0); scalarField dmHfgNbr(patch().size(), 0.0); if (!fluid_) { - myKDeltaNbr = nbrField.myKDelta(); - mpp.distribute(myKDeltaNbr); - mpCpTpNbr = nbrField.mpCpTp(); mpp.distribute(mpCpTpNbr); @@ -690,11 +694,11 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs() const scalarField mpCpdt(mpCpTpNbr + mpCpTp_); // Qr > 0 (heat up the wall) - scalarField alpha(KDeltaNbr + mpCpdt - (Qr + QrNbr + dmHfg)/Tp); + scalarField alpha(KDeltaNbr + mpCpdt - (Qr + QrNbr)/Tp); valueFraction() = alpha/(alpha + myKDelta_); - refValue() = (KDeltaNbr*nbrIntFld + mpCpdt*Tp)/alpha; + refValue() = (KDeltaNbr*nbrIntFld + mpCpdt*TpOld + dmHfg)/alpha; mixedFvPatchScalarField::updateCoeffs(); From 162a0eac4d3cc99e5388553501911c70f293ba72 Mon Sep 17 00:00:00 2001 From: sergio Date: Tue, 13 Dec 2016 13:30:42 -0800 Subject: [PATCH 7/7] Adding autoPtr clone constructors and correcting T.oldTime() --- .../humidityTemperatureCoupledMixedFvPatchScalarField.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C b/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C index 3d5f1e8e5c..c97553d74d 100644 --- a/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/properties/liquidPropertiesFvPatchFields/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C @@ -192,7 +192,7 @@ humidityTemperatureCoupledMixedFvPatchScalarField QrNbrName_(psf.QrNbrName_), QrName_(psf.QrName_), specieName_(psf.specieName_), - liquid_(psf.liquid_), + liquid_(psf.liquid_, false), liquidDict_(psf.liquidDict_), mass_(psf.mass_, mapper), Tvap_(psf.Tvap_), @@ -351,7 +351,7 @@ humidityTemperatureCoupledMixedFvPatchScalarField QrNbrName_(psf.QrNbrName_), QrName_(psf.QrName_), specieName_(psf.specieName_), - liquid_(psf.liquid_), + liquid_(psf.liquid_, false), liquidDict_(psf.liquidDict_), mass_(psf.mass_), Tvap_(psf.Tvap_), @@ -455,7 +455,7 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs() const volScalarField& T = static_cast(dimensionedInternalField()); - const fvPatchField& TpOld = T.boundaryField()[patch().index()]; + const scalarField TpOld(T.oldTime().boundaryField()[patch().index()]); scalarField Tin(patchInternalField());