diff --git a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C index dd05bf88aa..30a5f6fb1a 100644 --- a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C +++ b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C @@ -50,7 +50,7 @@ int main(int argc, char *argv[]) Info<< "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C index e2b1d536d2..1a5f871d96 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C @@ -124,6 +124,41 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap } +Foam::tmp +Foam::solidWallHeatFluxTemperatureFvPatchScalarField::K() const +{ + const fvMesh& mesh = patch().boundaryMesh().mesh(); + + if (mesh.objectRegistry::foundObject(KName_)) + { + return patch().lookupPatchField(KName_); + } + else if (mesh.objectRegistry::foundObject(KName_)) + { + const symmTensorField& KWall = + patch().lookupPatchField(KName_); + + vectorField n = patch().nf(); + + return n & KWall & n; + } + else + { + FatalErrorIn + ( + "solidWallHeatFluxTemperatureFvPatchScalarField::K()" + " const" + ) << "Did not find field " << KName_ + << " on mesh " << mesh.name() << " patch " << patch().name() + << endl + << "Please set 'K' to a valid volScalarField" + << " or a valid volSymmTensorField." << exit(FatalError); + + return scalarField(0); + } +} + + void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs() { if (updated()) @@ -131,12 +166,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs() return; } - const scalarField& Kw = patch().lookupPatchField - ( - KName_ - ); - - gradient() = q_/Kw; + gradient() = q_/K(); fixedGradientFvPatchScalarField::updateCoeffs(); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H index ffbcf52ed5..0544b2bffa 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H @@ -31,9 +31,10 @@ Description myWallPatch { type solidWallHeatFluxTemperature; - K K; // Name of K field - q uniform 1000; // Heat flux / [W/m2] - value 300.0; // Initial temperature / [K] + K K; // Name of K field + q uniform 1000; // Heat flux / [W/m2] + value uniform 300.0; // Initial temperature / [K] + gradient uniform 0.0; // Initial gradient / [K/m] } @@ -140,6 +141,11 @@ public: // Member functions + // Helper + + //- Get K field on this patch + tmp K() const; + // Evaluation functions //- Update the coefficients associated with the patch field diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C index 97992e8192..a18298f826 100644 --- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C +++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C @@ -78,7 +78,7 @@ int main(int argc, char *argv[]) Info<< "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C index 36bccb3911..35525f65f4 100644 --- a/src/OpenFOAM/db/Time/Time.C +++ b/src/OpenFOAM/db/Time/Time.C @@ -508,19 +508,11 @@ bool Foam::Time::run() const } } - return running; -} - - -bool Foam::Time::loop() -{ - bool running = run(); - if (running) { if (!subCycling_) { - readModifiedObjects(); + const_cast(*this).readModifiedObjects(); if (timeIndex_ == startTimeIndex_) { @@ -532,14 +524,22 @@ bool Foam::Time::loop() } } - // Check update the "running" status following the "++" operation - // to take into account possible side-effects from functionObjects - running = run(); + // Update the "running" status following the + // possible side-effects from functionObjects + running = value() < (endTime_ - 0.5*deltaT_); + } - if (running) - { - operator++(); - } + return running; +} + + +bool Foam::Time::loop() +{ + bool running = run(); + + if (running) + { + operator++(); } return running; diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C index 705505a716..7fa8d4b8bc 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C @@ -72,7 +72,8 @@ GeometricBoundaryField ( const BoundaryMesh& bmesh, const DimensionedField& field, - const wordList& patchFieldTypes + const wordList& patchFieldTypes, + const wordList& constraintTypes ) : FieldField(bmesh.size()), @@ -83,18 +84,22 @@ GeometricBoundaryField Info<< "GeometricField::" "GeometricBoundaryField::" "GeometricBoundaryField(const BoundaryMesh&, " - "const Field&, const wordList&)" + "const Field&, const wordList&, const wordList&)" << endl; } - if (patchFieldTypes.size() != this->size()) + if + ( + patchFieldTypes.size() != this->size() + || (constraintTypes.size() && (constraintTypes.size() != this->size())) + ) { FatalErrorIn ( "GeometricField::" "GeometricBoundaryField::" "GeometricBoundaryField(const BoundaryMesh&, " - "const Field&, const wordList&)" + "const Field&, const wordList&, const wordList&)" ) << "Incorrect number of patch type specifications given" << nl << " Number of patches in mesh = " << bmesh.size() << " number of patch type specifications = " @@ -102,18 +107,38 @@ GeometricBoundaryField << abort(FatalError); } - forAll(bmesh_, patchi) + if (constraintTypes.size()) { - set - ( - patchi, - PatchField::New + forAll(bmesh_, patchi) + { + set ( - patchFieldTypes[patchi], - bmesh_[patchi], - field - ) - ); + patchi, + PatchField::New + ( + patchFieldTypes[patchi], + constraintTypes[patchi], + bmesh_[patchi], + field + ) + ); + } + } + else + { + forAll(bmesh_, patchi) + { + set + ( + patchi, + PatchField::New + ( + patchFieldTypes[patchi], + bmesh_[patchi], + field + ) + ); + } } } diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C index 8c17c3886b..b4960c341d 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C @@ -231,14 +231,15 @@ Foam::GeometricField::GeometricField const IOobject& io, const Mesh& mesh, const dimensionSet& ds, - const wordList& patchFieldTypes + const wordList& patchFieldTypes, + const wordList& actualPatchTypes ) : DimensionedField(io, mesh, ds), timeIndex_(this->time().timeIndex()), field0Ptr_(NULL), fieldPrevIterPtr_(NULL), - boundaryField_(mesh.boundary(), *this, patchFieldTypes) + boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes) { if (debug) { @@ -287,14 +288,15 @@ Foam::GeometricField::GeometricField const IOobject& io, const Mesh& mesh, const dimensioned& dt, - const wordList& patchFieldTypes + const wordList& patchFieldTypes, + const wordList& actualPatchTypes ) : DimensionedField(io, mesh, dt), timeIndex_(this->time().timeIndex()), field0Ptr_(NULL), fieldPrevIterPtr_(NULL), - boundaryField_(mesh.boundary(), *this, patchFieldTypes) + boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes) { if (debug) { @@ -653,14 +655,22 @@ Foam::GeometricField::GeometricField ( const IOobject& io, const GeometricField& gf, - const wordList& patchFieldTypes + const wordList& patchFieldTypes, + const wordList& actualPatchTypes + ) : DimensionedField(io, gf), timeIndex_(gf.timeIndex()), field0Ptr_(NULL), fieldPrevIterPtr_(NULL), - boundaryField_(this->mesh().boundary(), *this, patchFieldTypes) + boundaryField_ + ( + this->mesh().boundary(), + *this, + patchFieldTypes, + actualPatchTypes + ) { if (debug) { diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H index 7d529ab40f..fa8ec179e4 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H @@ -128,12 +128,14 @@ public: //- Construct from a BoundaryMesh, // reference to the internal field - // and a wordList of patch types + // and a wordList of patch types and optional the actual patch + // types (to override constraint patches) GeometricBoundaryField ( const BoundaryMesh&, const DimensionedInternalField&, - const wordList& + const wordList& wantedPatchTypes, + const wordList& actualPatchTypes = wordList() ); //- Construct from a BoundaryMesh, @@ -283,7 +285,8 @@ public: const IOobject&, const Mesh&, const dimensionSet&, - const wordList& patchFieldTypes + const wordList& wantedPatchTypes, + const wordList& actualPatchTypes = wordList() ); //- Constructor given IOobject, mesh, dimensioned and patch type. @@ -301,7 +304,8 @@ public: const IOobject&, const Mesh&, const dimensioned&, - const wordList& patchFieldTypes + const wordList& wantedPatchTypes, + const wordList& actualPatchTypes = wordList() ); //- Constructor from components @@ -387,7 +391,8 @@ public: ( const IOobject&, const GeometricField&, - const wordList& patchFieldTypes + const wordList& patchFieldTypes, + const wordList& actualPatchTypes = wordList() ); diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H index 98604bb4bf..869c536365 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H @@ -35,7 +35,7 @@ Description SourceFiles pointPatchField.C - newpointPatchField.C + newPointPatchField.C \*---------------------------------------------------------------------------*/ @@ -199,6 +199,18 @@ public: const DimensionedField& ); + //- Return a pointer to a new patchField created on freestore given + // patch and internal field + // (does not set the patch field values). + // Allows override of constraint type + static autoPtr > New + ( + const word&, + const word& actualPatchType, + const pointPatch&, + const DimensionedField& + ); + //- Return a pointer to a new patchField created on freestore from // a given pointPatchField mapped onto a new patch static autoPtr > New diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldNew.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldNew.C index b14fbeac8f..749511b2cf 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldNew.C +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldNew.C @@ -29,6 +29,7 @@ template Foam::autoPtr > Foam::pointPatchField::New ( const word& patchFieldType, + const word& actualPatchType, const pointPatch& p, const DimensionedField& iF ) @@ -36,7 +37,8 @@ Foam::autoPtr > Foam::pointPatchField::New if (debug) { Info<< "PointPatchField::" - "New(const word&, const pointPatch&, const Field&) : " + "New(const word&, const word&" + ", const pointPatch&, const Field&) : " "constructing pointPatchField" << endl; } @@ -49,7 +51,7 @@ Foam::autoPtr > Foam::pointPatchField::New FatalErrorIn ( "PointPatchField::New" - "(const word&, const pointPatch&, const Field&)" + "(const word&, const word&, const pointPatch&, const Field&)" ) << "Unknown patchFieldType type " << patchFieldType << endl << endl @@ -60,31 +62,48 @@ Foam::autoPtr > Foam::pointPatchField::New autoPtr > pfPtr(cstrIter()(p, iF)); - if (pfPtr().constraintType() == p.constraintType()) + if + ( + actualPatchType == word::null + || actualPatchType != p.type() + ) { - // Compatible (constraint-wise) with the patch type - return pfPtr; - } - else - { - // Use default constraint type - typename pointPatchConstructorTable::iterator patchTypeCstrIter = - pointPatchConstructorTablePtr_->find(p.type()); - - if (patchTypeCstrIter == pointPatchConstructorTablePtr_->end()) + if (pfPtr().constraintType() != p.constraintType()) { - FatalErrorIn - ( - "PointPatchField::New" - "(const word&, const pointPatch&, const Field&)" - ) << "inconsistent patch and patchField types for \n" - << " patch type " << p.type() - << " and patchField type " << patchFieldType - << exit(FatalError); - } + // Use default constraint type + typename pointPatchConstructorTable::iterator patchTypeCstrIter = + pointPatchConstructorTablePtr_->find(p.type()); - return patchTypeCstrIter()(p, iF); + if (patchTypeCstrIter == pointPatchConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "PointPatchField::New" + "(const word&, const word&" + ", const pointPatch&, const Field&)" + ) << "inconsistent patch and patchField types for \n" + << " patch type " << p.type() + << " and patchField type " << patchFieldType + << exit(FatalError); + } + + return patchTypeCstrIter()(p, iF); + } } + + return pfPtr; +} + + +template +Foam::autoPtr > Foam::pointPatchField::New +( + const word& patchFieldType, + const pointPatch& p, + const DimensionedField& iF +) +{ + return New(patchFieldType, word::null, p, iF); } diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H index e43423b1b1..69ee851768 100644 --- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H @@ -226,6 +226,18 @@ public: const DimensionedField& ); + //- Return a pointer to a new patchField created on freestore given + // patch and internal field + // (does not set the patch field values). + // Allows override of constraint type + static tmp > New + ( + const word&, + const word& actualPatchType, + const fvPatch&, + const DimensionedField& + ); + //- Return a pointer to a new patchField created on freestore from // a given fvPatchField mapped onto a new patch static tmp > New diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchFieldNew.C b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchFieldNew.C index 8f89f202d9..86884666a9 100644 --- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchFieldNew.C +++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchFieldNew.C @@ -29,14 +29,16 @@ template Foam::tmp > Foam::fvPatchField::New ( const word& patchFieldType, + const word& actualPatchType, const fvPatch& p, const DimensionedField& iF ) { if (debug) { - Info<< "fvPatchField::New(const word&, const fvPatch&, " - "const DimensionedField&) : patchFieldType=" + Info<< "fvPatchField::New(const word&, const word&, " + "const fvPatch&, const DimensionedField&) :" + " patchFieldType=" << patchFieldType << endl; } @@ -48,7 +50,7 @@ Foam::tmp > Foam::fvPatchField::New { FatalErrorIn ( - "fvPatchField::New(const word&, const fvPatch&, " + "fvPatchField::New(const word&, const word&, const fvPatch&, " "const DimensionedField&)" ) << "Unknown patchTypefield type " << patchFieldType << endl << endl @@ -57,12 +59,23 @@ Foam::tmp > Foam::fvPatchField::New << exit(FatalError); } - typename patchConstructorTable::iterator patchTypeCstrIter = - patchConstructorTablePtr_->find(p.type()); - - if (patchTypeCstrIter != patchConstructorTablePtr_->end()) + if + ( + actualPatchType == word::null + || actualPatchType != p.type() + ) { - return patchTypeCstrIter()(p, iF); + typename patchConstructorTable::iterator patchTypeCstrIter = + patchConstructorTablePtr_->find(p.type()); + + if (patchTypeCstrIter != patchConstructorTablePtr_->end()) + { + return patchTypeCstrIter()(p, iF); + } + else + { + return cstrIter()(p, iF); + } } else { @@ -71,6 +84,18 @@ Foam::tmp > Foam::fvPatchField::New } +template +Foam::tmp > Foam::fvPatchField::New +( + const word& patchFieldType, + const fvPatch& p, + const DimensionedField& iF +) +{ + return New(patchFieldType, word::null, p, iF); +} + + template Foam::tmp > Foam::fvPatchField::New ( diff --git a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H index 85d139ac74..b5011eefcd 100644 --- a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H +++ b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H @@ -216,6 +216,18 @@ public: const DimensionedField& ); + //- Return a pointer to a new patchField created on freestore given + // patch and internal field + // (does not set the patch field values) + // Allows override of constraint type + static tmp > New + ( + const word&, + const word& actualPatchType, + const fvPatch&, + const DimensionedField& + ); + //- Return a pointer to a new patchField created on freestore from // a given fvsPatchField mapped onto a new patch static tmp > New diff --git a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchFieldNew.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchFieldNew.C index b33ee35043..b8ad302fdd 100644 --- a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchFieldNew.C +++ b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchFieldNew.C @@ -34,14 +34,15 @@ template tmp > fvsPatchField::New ( const word& patchFieldType, + const word& actualPatchType, const fvPatch& p, const DimensionedField& iF ) { if (debug) { - Info<< "fvsPatchField::New(const word&, const fvPatch&, " - "const Field&) : " + Info<< "fvsPatchField::New(const word&, const word&" + ", const fvPatch&, const Field&) : " "constructing fvsPatchField" << endl; } @@ -53,8 +54,8 @@ tmp > fvsPatchField::New { FatalErrorIn ( - "fvsPatchField::New(const word&, const fvPatch&, " - "const Field&)" + "fvsPatchField::New(const word&, const word&, const fvPatch&" + ", const Field&)" ) << "Unknown patchTypefield type " << patchFieldType << endl << endl << "Valid patchField types are :" << endl @@ -62,12 +63,23 @@ tmp > fvsPatchField::New << exit(FatalError); } - typename patchConstructorTable::iterator patchTypeCstrIter = - patchConstructorTablePtr_->find(p.type()); - - if (patchTypeCstrIter != patchConstructorTablePtr_->end()) + if + ( + actualPatchType == word::null + || actualPatchType != p.type() + ) { - return patchTypeCstrIter()(p, iF); + typename patchConstructorTable::iterator patchTypeCstrIter = + patchConstructorTablePtr_->find(p.type()); + + if (patchTypeCstrIter != patchConstructorTablePtr_->end()) + { + return patchTypeCstrIter()(p, iF); + } + else + { + return cstrIter()(p, iF); + } } else { @@ -76,6 +88,18 @@ tmp > fvsPatchField::New } +template +tmp > fvsPatchField::New +( + const word& patchFieldType, + const fvPatch& p, + const DimensionedField& iF +) +{ + return New(patchFieldType, word::null, p, iF); +} + + template tmp > fvsPatchField::New ( diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index a381c9d431..e17e7138d6 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -160,10 +160,27 @@ const Foam::vector Foam::KinematicParcel::calcVelocity const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL; // Momentum source due to particle forces - const vector FCoupled = - mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U); - const vector FNonCoupled = - mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U); + const vector FCoupled = mass*td.cloud().forces().calcCoupled + ( + cellI, + dt, + rhoc_, + rho, + Uc_, + U, + d + ); + + const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled + ( + cellI, + dt, + rhoc_, + rho, + Uc_, + U, + d + ); // New particle velocity diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C index 265f9b59ae..7778736683 100644 --- a/src/lagrangian/intermediate/particleForces/particleForces.C +++ b/src/lagrangian/intermediate/particleForces/particleForces.C @@ -27,6 +27,8 @@ License #include "fvMesh.H" #include "volFields.H" #include "fvcGrad.H" +#include "mathematicalConstants.H" +#include "electromagneticConstants.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -45,12 +47,20 @@ Foam::particleForces::particleForces virtualMass_(dict_.lookup("virtualMass")), Cvm_(0.0), pressureGradient_(dict_.lookup("pressureGradient")), - UName_(dict_.lookupOrDefault("U", "U")) + paramagnetic_(dict_.lookup("paramagnetic")), + magneticSusceptibility_(0.0), + UName_(dict_.lookupOrDefault("U", "U")), + HdotGradHName_(dict_.lookupOrDefault("HdotGradH", "HdotGradH")) { if (virtualMass_) { dict_.lookup("Cvm") >> Cvm_; } + + if (paramagnetic_) + { + dict_.lookup("magneticSusceptibility") >> magneticSusceptibility_; + } } @@ -64,7 +74,10 @@ Foam::particleForces::particleForces(const particleForces& f) virtualMass_(f.virtualMass_), Cvm_(f.Cvm_), pressureGradient_(f.pressureGradient_), - UName_(f.UName_) + paramagnetic_(f.paramagnetic_), + magneticSusceptibility_(f.magneticSusceptibility_), + UName_(f.UName_), + HdotGradHName_(f.HdotGradHName_) {} @@ -102,18 +115,42 @@ Foam::Switch Foam::particleForces::virtualMass() const } +Foam::scalar Foam::particleForces::Cvm() const +{ + return Cvm_; +} + + Foam::Switch Foam::particleForces::pressureGradient() const { return pressureGradient_; } +Foam::Switch Foam::particleForces::paramagnetic() const +{ + return paramagnetic_; +} + + +Foam::scalar Foam::particleForces::magneticSusceptibility() const +{ + return magneticSusceptibility_; +} + + const Foam::word& Foam::particleForces::UName() const { return UName_; } +const Foam::word& Foam::particleForces::HdotGradHName() const +{ + return HdotGradHName_; +} + + void Foam::particleForces::cacheFields(const bool store) { if (store && pressureGradient_) @@ -139,10 +176,11 @@ Foam::vector Foam::particleForces::calcCoupled const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const { - vector Ftot = vector::zero; + vector accelTot = vector::zero; // Virtual mass force if (virtualMass_) @@ -151,17 +189,17 @@ Foam::vector Foam::particleForces::calcCoupled ( "Foam::particleForces::calcCoupled(...) - virtual mass force" ); -// Ftot += Cvm_*rhoc/rho*d(Uc - U)/dt; +// accelTot += Cvm_*rhoc/rho*d(Uc - U)/dt; } // Pressure gradient force if (pressureGradient_) { const volTensorField& gradU = *gradUPtr_; - Ftot += rhoc/rho*(U & gradU[cellI]); + accelTot += rhoc/rho*(U & gradU[cellI]); } - return Ftot; + return accelTot; } @@ -172,20 +210,47 @@ Foam::vector Foam::particleForces::calcNonCoupled const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const { - vector Ftot = vector::zero; + vector accelTot = vector::zero; // Gravity force if (gravity_) { - Ftot += g_*(1.0 - rhoc/rho); + accelTot += g_*(1.0 - rhoc/rho); } - return Ftot; + // Magnetic field force + + if (paramagnetic_) + { + const volVectorField& HdotGradH = mesh_.lookupObject + ( + HdotGradHName_ + ); + + accelTot += + 3.0*constant::electromagnetic::mu0.value()/rho + *magneticSusceptibility_/(magneticSusceptibility_ + 3) + *HdotGradH[cellI]; + + // force is: + + // 4.0 + // *constant::mathematical::pi + // *constant::electromagnetic::mu0.value() + // *pow3(d/2) + // *magneticSusceptibility_/(magneticSusceptibility_ + 3) + // *HdotGradH[cellI]; + + // which is divided by mass ((4/3)*pi*r^3*rho) to produce + // acceleration + } + + return accelTot; } // ************************************************************************* // - diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H index 87abf7a7d0..cc47366533 100644 --- a/src/lagrangian/intermediate/particleForces/particleForces.H +++ b/src/lagrangian/intermediate/particleForces/particleForces.H @@ -84,12 +84,22 @@ class particleForces //- Pressure gradient Switch pressureGradient_; + //- Paramagnetic force + Switch paramagnetic_; + + //- Magnetic susceptibility of particle + scalar magneticSusceptibility_; + // Additional info - //- Name of velucity field - default = "U" + //- Name of velocity field - default = "U" const word UName_; + //- Name of paramagnetic field strength field - default = + // "HdotGradH" + const word HdotGradHName_; + public: @@ -128,7 +138,7 @@ public: Switch virtualMass() const; //- Return virtual mass force coefficient - Switch Cvm() const; + scalar Cvm() const; //- Return pressure gradient force activate switch Switch pressureGradient() const; @@ -136,6 +146,15 @@ public: //- Return name of velocity field const word& UName() const; + //- Return paramagnetic force activate switch + Switch paramagnetic() const; + + //- Return magnetic susceptibility + scalar magneticSusceptibility() const; + + //- Return name of paramagnetic field strength field + const word& HdotGradHName() const; + // Evaluation @@ -150,7 +169,8 @@ public: const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const; //- Calculate external forces applied to the particles @@ -161,7 +181,8 @@ public: const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const; }; diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C index f25f0ef448..844272e98b 100644 --- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C +++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C @@ -215,6 +215,43 @@ void Foam::decompositionMethod::calcCellCells } +// Return the minimum face between two cells. Only relevant for +// cells with multiple faces inbetween. +Foam::label Foam::decompositionMethod::masterFace +( + const polyMesh& mesh, + const label own, + const label nei +) +{ + label minFaceI = labelMax; + + // Count multiple faces between own and nei only once + const cell& ownFaces = mesh.cells()[own]; + forAll(ownFaces, i) + { + label otherFaceI = ownFaces[i]; + + if (mesh.isInternalFace(otherFaceI)) + { + label nbrCellI = + ( + mesh.faceNeighbour()[otherFaceI] != own + ? mesh.faceNeighbour()[otherFaceI] + : mesh.faceOwner()[otherFaceI] + ); + + if (nbrCellI == nei) + { + minFaceI = min(minFaceI, otherFaceI); + } + } + } + + return minFaceI; +} + + void Foam::decompositionMethod::calcCSR ( const polyMesh& mesh, @@ -222,61 +259,15 @@ void Foam::decompositionMethod::calcCSR List& xadj ) { + const polyBoundaryMesh& pbm = mesh.boundaryMesh(); + // Make Metis CSR (Compressed Storage Format) storage // adjncy : contains neighbours (= edges in graph) // xadj(celli) : start of information in adjncy for celli - xadj.setSize(mesh.nCells()+1); - // Initialise the number of internal faces of the cells to twice the - // number of internal faces - label nInternalFaces = 2*mesh.nInternalFaces(); - - // Check the boundary for coupled patches and add to the number of - // internal faces - const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - - forAll(pbm, patchi) - { - if (isA(pbm[patchi])) - { - nInternalFaces += pbm[patchi].size(); - } - } - - // Create the adjncy array the size of the total number of internal and - // coupled faces - adjncy.setSize(nInternalFaces); - - // Fill in xadj - // ~~~~~~~~~~~~ - label freeAdj = 0; - - for (label cellI = 0; cellI < mesh.nCells(); cellI++) - { - xadj[cellI] = freeAdj; - - const labelList& cFaces = mesh.cells()[cellI]; - - forAll(cFaces, i) - { - label faceI = cFaces[i]; - - if - ( - mesh.isInternalFace(faceI) - || isA(pbm[pbm.whichPatch(faceI)]) - ) - { - freeAdj++; - } - } - } - xadj[mesh.nCells()] = freeAdj; - - - // Fill in adjncy - // ~~~~~~~~~~~~~~ + // Count unique faces between cells + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ labelList nFacesPerCell(mesh.nCells(), 0); @@ -286,26 +277,95 @@ void Foam::decompositionMethod::calcCSR label own = mesh.faceOwner()[faceI]; label nei = mesh.faceNeighbour()[faceI]; - adjncy[xadj[own] + nFacesPerCell[own]++] = nei; - adjncy[xadj[nei] + nFacesPerCell[nei]++] = own; + if (faceI == masterFace(mesh, own, nei)) + { + nFacesPerCell[own]++; + nFacesPerCell[nei]++; + } } // Coupled faces. Only cyclics done. - forAll(pbm, patchi) + HashSet > cellPair(mesh.nFaces()-mesh.nInternalFaces()); + + forAll(pbm, patchI) { - if (isA(pbm[patchi])) + if (isA(pbm[patchI])) { - const unallocLabelList& faceCells = pbm[patchi].faceCells(); + const unallocLabelList& faceCells = pbm[patchI].faceCells(); label sizeby2 = faceCells.size()/2; - for (label facei=0; facei(pbm[patchI])) + { + const unallocLabelList& faceCells = pbm[patchI].faceCells(); + + label sizeby2 = faceCells.size()/2; + + for (label faceI=0; faceI& xadj ) { + labelHashSet nbrCells; + // Count number of internal faces label nConnections = 0; forAll(cellCells, coarseI) { - nConnections += cellCells[coarseI].size(); + nbrCells.clear(); + + const labelList& cCells = cellCells[coarseI]; + + forAll(cCells, i) + { + if (nbrCells.insert(cCells[i])) + { + nConnections++; + } + } } // Create the adjncy array as twice the size of the total number of @@ -343,11 +415,16 @@ void Foam::decompositionMethod::calcCSR { xadj[coarseI] = freeAdj; + nbrCells.clear(); + const labelList& cCells = cellCells[coarseI]; forAll(cCells, i) { - adjncy[freeAdj++] = cCells[i]; + if (nbrCells.insert(cCells[i])) + { + adjncy[freeAdj++] = cCells[i]; + } } } xadj[cellCells.size()] = freeAdj; @@ -413,15 +490,21 @@ void Foam::decompositionMethod::calcDistributedCSR // Number of faces per cell List nFacesPerCell(mesh.nCells(), 0); - // Number of coupled faces - label nCoupledFaces = 0; - for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) { - nFacesPerCell[faceOwner[faceI]]++; - nFacesPerCell[faceNeighbour[faceI]]++; + label own = faceOwner[faceI]; + label nei = faceNeighbour[faceI]; + + if (faceI == masterFace(mesh, own, nei)) + { + nFacesPerCell[own]++; + nFacesPerCell[nei]++; + } } + // Handle coupled faces + HashSet > cellPair(mesh.nFaces()-mesh.nInternalFaces()); + forAll(patches, patchI) { const polyPatch& pp = patches[patchI]; @@ -429,11 +512,18 @@ void Foam::decompositionMethod::calcDistributedCSR if (pp.coupled()) { label faceI = pp.start(); + label bFaceI = pp.start()-mesh.nInternalFaces(); forAll(pp, i) { - nCoupledFaces++; - nFacesPerCell[faceOwner[faceI++]]++; + label own = faceOwner[faceI]; + label globalNei = globalNeighbour[bFaceI]; + if (cellPair.insert(edge(own, globalNei))) + { + nFacesPerCell[own]++; + } + faceI++; + bFaceI++; } } } @@ -459,7 +549,7 @@ void Foam::decompositionMethod::calcDistributedCSR // Fill in adjncy // ~~~~~~~~~~~~~~ - adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces); + adjncy.setSize(freeAdj); nFacesPerCell = 0; @@ -469,10 +559,17 @@ void Foam::decompositionMethod::calcDistributedCSR label own = faceOwner[faceI]; label nei = faceNeighbour[faceI]; - adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei); - adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own); + if (faceI == masterFace(mesh, own, nei)) + { + adjncy[xadj[own] + nFacesPerCell[own]++] = + globalCells.toGlobal(nei); + adjncy[xadj[nei] + nFacesPerCell[nei]++] = + globalCells.toGlobal(own); + } } + // For boundary faces is offsetted coupled neighbour + cellPair.clear(); forAll(patches, patchI) { const polyPatch& pp = patches[patchI]; @@ -485,9 +582,11 @@ void Foam::decompositionMethod::calcDistributedCSR forAll(pp, i) { label own = faceOwner[faceI]; - adjncy[xadj[own] + nFacesPerCell[own]++] = - globalNeighbour[bFaceI]; - + label globalNei = globalNeighbour[bFaceI]; + if (cellPair.insert(edge(own, globalNei))) + { + adjncy[xadj[own] + nFacesPerCell[own]++] = globalNei; + } faceI++; bFaceI++; } diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H index d729fe9bd3..dabca144c8 100644 --- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H +++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H @@ -65,6 +65,10 @@ protected: labelListList& cellCells ); + //- Calculate the minimum face between two neighbouring cells + // (usually there is only one) + static label masterFace(const polyMesh&, const label, const label); + // From mesh to compact row storage format // (like CompactListList) static void calcCSR diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files index bfed281015..df93e62853 100644 --- a/src/postProcessing/functionObjects/field/Make/files +++ b/src/postProcessing/functionObjects/field/Make/files @@ -19,6 +19,6 @@ streamLine/streamLine.C streamLine/streamLineParticle.C streamLine/streamLineParticleCloud.C streamLine/streamLineFunctionObject.C - +streamLine/vectorIOFieldField.C LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C index 9d65cefc57..97cd9393bd 100644 --- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C +++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C @@ -80,10 +80,11 @@ void Foam::fieldValues::cellSource::setCellZoneCells() } cellId_.setSize(count); + nCells_ = returnReduce(cellId_.size(), sumOp