From def1e7b59f03dac669671c6d4cd35fdee86faa76 Mon Sep 17 00:00:00 2001 From: graham Date: Tue, 8 Sep 2009 21:07:22 +0100 Subject: [PATCH] Adding leapfrog integration parts to the move method. --- .../ReferredCellList/ReferredCellList.C | 2 +- .../InteractingKinematicCloud.C | 26 +++--- .../InteractingKinematicParcel.C | 83 ++++++++++++------- .../InteractingKinematicParcel.H | 28 ++++++- .../InteractingKinematicParcelI.H | 24 +++++- .../PairCollision/PairCollision.C | 2 +- 6 files changed, 117 insertions(+), 48 deletions(-) diff --git a/src/lagrangian/basic/InteractionLists/ReferredCellList/ReferredCellList.C b/src/lagrangian/basic/InteractionLists/ReferredCellList/ReferredCellList.C index ec88f0c4d0..1f5cd9e194 100644 --- a/src/lagrangian/basic/InteractionLists/ReferredCellList/ReferredCellList.C +++ b/src/lagrangian/basic/InteractionLists/ReferredCellList/ReferredCellList.C @@ -1620,7 +1620,7 @@ void Foam::ReferredCellList::referParticles const List >& cellOccupancy ) { - Info<< "Refer particles" << endl; + Info<< " Refer particles" << endl; // Clear all existing referred particles diff --git a/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C index c05b0fe4d3..98cdd756d6 100644 --- a/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C @@ -254,23 +254,25 @@ void Foam::InteractingKinematicCloud::evolve() resetSourceTerms(); } - const scalar deltaT = mesh_.time().deltaT().value(); - - forAllIter(typename InteractingKinematicCloud, *this, iter) - { - ParcelType& p = iter(); - p.U() += 0.5*deltaT*p.f()/p.mass(); - } + // Sympletic leapfrog integration of particle forces: + // + apply half deltaV with stored force + // + move positions with new velocity + // + calculate forces in new position + // + apply half deltaV with new force + td.part() = ParcelType::trackData::LEAPFROG_VELOCITY_STEP; Cloud::move(td); + td.part() = ParcelType::trackData::LINEAR_TRACK; + Cloud::move(td); + + // td.part() = ParcelType::trackData::ROTATIONAL_TRACK; + // Cloud::move(td); + this->collision().collide(); - forAllIter(typename InteractingKinematicCloud, *this, iter) - { - ParcelType& p = iter(); - p.U() += 0.5*deltaT*p.f()/p.mass(); - } + td.part() = ParcelType::trackData::LEAPFROG_VELOCITY_STEP; + Cloud::move(td); postEvolve(); } diff --git a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C index 48214ecd01..0892e64db7 100644 --- a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C @@ -229,48 +229,73 @@ bool Foam::InteractingKinematicParcel::move(TrackData& td) const polyBoundaryMesh& pbMesh = mesh.boundaryMesh(); const scalar deltaT = mesh.time().deltaT().value(); - scalar tEnd = (1.0 - p.stepFraction())*deltaT; - const scalar dtMax = tEnd; - while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL) + if (td.part() == TrackData::LEAPFROG_VELOCITY_STEP) { - // Apply correction to position for reduced-D cases - meshTools::constrainToMeshCentre(mesh, p.position()); + // First and last leapfrog velocity adjust part, required + // before and after tracking and force calculation - // Set the Lagrangian time-step - scalar dt = min(dtMax, tEnd); + p.U() += 0.5*deltaT*p.f()/p.mass(); + } + else if (td.part() == TrackData::LINEAR_TRACK) + { + scalar tEnd = (1.0 - p.stepFraction())*deltaT; + const scalar dtMax = tEnd; - // Remember which cell the Parcel is in since this will change if a - // face is hit - label cellI = p.cell(); - - dt *= p.trackToFace(p.position() + dt*U_, td); - - tEnd -= dt; - p.stepFraction() = 1.0 - tEnd/deltaT; - - // Avoid problems with extremely small timesteps - if (dt > ROOTVSMALL) + while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL) { - // Update cell based properties - p.setCellValues(td, dt, cellI); + // Apply correction to position for reduced-D cases + meshTools::constrainToMeshCentre(mesh, p.position()); - if (td.cloud().cellValueSourceCorrection()) + // Set the Lagrangian time-step + scalar dt = min(dtMax, tEnd); + + // Remember which cell the Parcel is in since this will change if a + // face is hit + label cellI = p.cell(); + + dt *= p.trackToFace(p.position() + dt*U_, td); + + tEnd -= dt; + p.stepFraction() = 1.0 - tEnd/deltaT; + + // Avoid problems with extremely small timesteps + if (dt > ROOTVSMALL) { - p.cellValueSourceCorrection(td, dt, cellI); + // Update cell based properties + p.setCellValues(td, dt, cellI); + + if (td.cloud().cellValueSourceCorrection()) + { + p.cellValueSourceCorrection(td, dt, cellI); + } + + p.calc(td, dt, cellI); } - p.calc(td, dt, cellI); - } - - if (p.onBoundary() && td.keepParticle) - { - if (isType(pbMesh[p.patch(p.face())])) + if (p.onBoundary() && td.keepParticle) { - td.switchProcessor = true; + if (isType(pbMesh[p.patch(p.face())])) + { + td.switchProcessor = true; + } } } } + else if (td.part() == TrackData::ROTATIONAL_TRACK) + { + Info<< "No rotational tracking implementation" << endl; + } + else + { + FatalErrorIn + ( + "InteractingKinematicParcel::move(TrackData& td)" + ) + << td.part() + << " is an invalid part of the tracking method." + << abort(FatalError); + } return td.keepParticle; } diff --git a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.H index 03a207a4d6..41a132df9e 100644 --- a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.H @@ -125,6 +125,18 @@ public: : public Particle::trackData { + public: + + enum trackPart + { + LEAPFROG_VELOCITY_STEP, + LINEAR_TRACK, + ROTATIONAL_TRACK + }; + + + private: + // Private data //- Reference to the cloud containing this particle @@ -133,7 +145,6 @@ public: //- Particle constant properties const constantProperties& constProps_; - // Interpolators for continuous phase fields //- Density interpolator @@ -148,6 +159,10 @@ public: //- Local gravitational or other body-force acceleration const vector& g_; + // label specifying which part of the integration + // algorithm is taking place + trackPart part_; + public: @@ -161,7 +176,8 @@ public: const interpolation& rhoInterp, const interpolation& UInterp, const interpolation& muInterp, - const vector& g + const vector& g, + trackPart part = LINEAR_TRACK ); @@ -185,8 +201,14 @@ public: // phase dynamic viscosity field inline const interpolation& muInterp() const; - // Return const access to the gravitational acceleration vector + //- Return const access to the gravitational acceleration vector inline const vector& g() const; + + //- Return the part of the tracking operation taking place + inline trackPart part() const; + + //- Return access to the part of the tracking operation taking place + inline trackPart& part(); }; diff --git a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcelI.H index 12dc8ff8bc..96276c09a7 100644 --- a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcelI.H +++ b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcelI.H @@ -55,7 +55,8 @@ inline Foam::InteractingKinematicParcel::trackData::trackData const interpolation& rhoInterp, const interpolation& UInterp, const interpolation& muInterp, - const vector& g + const vector& g, + trackPart part ) : Particle::trackData(cloud), @@ -64,7 +65,8 @@ inline Foam::InteractingKinematicParcel::trackData::trackData rhoInterp_(rhoInterp), UInterp_(UInterp), muInterp_(muInterp), - g_(g) + g_(g), + part_(part) {} @@ -205,6 +207,24 @@ Foam::InteractingKinematicParcel::trackData::g() const } +template +inline +typename Foam::InteractingKinematicParcel::trackData::trackPart +Foam::InteractingKinematicParcel::trackData::part() const +{ + return part_; +} + + +template +inline +typename Foam::InteractingKinematicParcel::trackData::trackPart& +Foam::InteractingKinematicParcel::trackData::part() +{ + return part_; +} + + // * * * * * * * InteractingKinematicParcel Member Functions * * * * * * * // template diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C index 54795112bb..5109072f0d 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C @@ -32,7 +32,7 @@ License template void Foam::PairCollision::buildCellOccupancy() { - Info<< "Build cell occupancy" << endl; + Info<< " Build cell occupancy" << endl; forAll(cellOccupancy_, cO) {