From 2cc2770530d4f11b33bb2bfb0d2ddffa83933b7e Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 23 Aug 2013 16:20:42 +0100 Subject: [PATCH 1/9] multiphaseEulerFoam: Avoid cyclic linkage dependency --- .../multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options index ce7e7a2d85..68bcd9ef06 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options @@ -10,5 +10,4 @@ EXE_INC = \ LIB_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ - -lcompressibleMultiphaseEulerianInterfacialModels \ -lfiniteVolume From 94c4d37f63c8f3a5efd16957d028e66af5ccbfeb Mon Sep 17 00:00:00 2001 From: laurence Date: Wed, 28 Aug 2013 12:55:22 +0100 Subject: [PATCH 2/9] ENH: foamyHexMesh: check for intersection in inside/outside test --- .../conformationSurfaces.C | 64 ++++++++----------- 1 file changed, 26 insertions(+), 38 deletions(-) diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C index 20a8b278eb..0360cf0d83 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C @@ -670,6 +670,32 @@ Foam::Field Foam::conformationSurfaces::wellInOutSide continue; } + const searchableSurface& surface(allGeometry_[surfaces_[s]]); + + if (!surface.hasVolumeType()) + { + pointField sample(1, samplePts[i]); + scalarField nearestDistSqr(1, GREAT); + List info; + + surface.findNearest(sample, nearestDistSqr, info); + + vector hitDir = info[0].rawPoint() - samplePts[i]; + hitDir /= mag(hitDir) + SMALL; + + if + ( + findSurfaceAnyIntersection + ( + samplePts[i], + info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir) + ) + ) + { + continue; + } + } + if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE) { if @@ -694,44 +720,6 @@ Foam::Field Foam::conformationSurfaces::wellInOutSide break; } } -// else -// { -// // Surface volume type is unknown -// Info<< "UNKNOWN" << endl; -// // Get nearest face normal -// -// pointField sample(1, samplePts[i]); -// scalarField nearestDistSqr(1, GREAT); -// List info; -// vectorField norms(1); -// -// surface.findNearest(sample, nearestDistSqr, info); -// surface.getNormal(info, norms); -// -// vector fN = norms[0]; -// fN /= mag(fN); -// -// vector hitDir = info[0].rawPoint() - samplePts[i]; -// hitDir /= mag(hitDir); -// -// if ((fN & hitDir) < 0) -// { -// // Point is OUTSIDE -// -// if -// ( -// normalVolumeTypes_[regionI] -// == extendedFeatureEdgeMesh::OUTSIDE -// ) -// { -// } -// else -// { -// insidePoint[i] = false; -// break; -// } -// } -// } } } From d185f8b9fe56dbf96c4cc560e0b6bb3ed40e02b9 Mon Sep 17 00:00:00 2001 From: Sergio Ferraris Date: Thu, 29 Aug 2013 11:05:27 +0100 Subject: [PATCH 3/9] ENH: Adding useChemistrySolver in reactingOneDim solid region to avoid using chemistry solvers in reacting solids Adding non const access function to chemistryModel.H for combustion models that calculates RR externally Removing EulerImplicit from options in solidReactions. BUG: Correcting function used by sequential solver in pyrolysisChemistryModel.C --- .../reactingOneDim/reactingOneDim.C | 26 ++++++++++---- .../reactingOneDim/reactingOneDim.H | 3 ++ .../basicChemistryModel/basicChemistryModel.H | 8 ++++- .../chemistryModel/chemistryModel.H | 6 ++++ .../chemistryModel/chemistryModelI.H | 12 ++++++- .../basicSolidChemistryModel.C | 31 +++++++++++++++++ .../basicSolidChemistryModel.H | 9 +++++ .../pyrolysisChemistryModel.C | 34 +++++++------------ .../pyrolysisChemistryModel.H | 13 ++++--- .../solidChemistryModel/solidChemistryModel.H | 22 ++++++------ .../solidChemistryModelI.H | 14 +------- .../makeSolidChemistrySolverType.H | 14 ++------ .../solidThermo/solidThermo/makeSolidThermo.H | 2 +- 13 files changed, 122 insertions(+), 72 deletions(-) diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C index 10c8a4b4ac..6fa334158a 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C +++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C @@ -61,6 +61,9 @@ void reactingOneDim::readReactingOneDimControls() coeffs().lookup("gasHSource") >> gasHSource_; coeffs().lookup("QrHSource") >> QrHSource_; + useChemistrySolvers_ = + coeffs().lookupOrDefault("useChemistrySolvers", true); + } @@ -462,7 +465,8 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh) totalGasMassFlux_(0.0), totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)), gasHSource_(false), - QrHSource_(false) + QrHSource_(false), + useChemistrySolvers_(true) { if (active_) { @@ -560,7 +564,8 @@ reactingOneDim::reactingOneDim totalGasMassFlux_(0.0), totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)), gasHSource_(false), - QrHSource_(false) + QrHSource_(false), + useChemistrySolvers_(true) { if (active_) { @@ -681,11 +686,18 @@ void reactingOneDim::evolveRegion() { Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl; - solidChemistry_->solve - ( - time().value() - time().deltaTValue(), - time().deltaTValue() - ); + if (useChemistrySolvers_) + { + solidChemistry_->solve + ( + time().value() - time().deltaTValue(), + time().deltaTValue() + ); + } + else + { + solidChemistry_->calculate(); + } solveContinuity(); diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H index 4a92fd6967..a5950a9804 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H +++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H @@ -154,6 +154,9 @@ protected: //- Add in depth radiation source term bool QrHSource_; + //- Use chemistry solvers (ode or sequential) + bool useChemistrySolvers_; + // Protected member functions diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H index 67704bb9df..75ea14eb67 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -140,6 +140,12 @@ public: const label i ) const = 0; + //- Return access to chemical source terms [kg/m3/s] + virtual DimensionedField& RR + ( + const label i + ) = 0; + // Chemistry solution diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H index 7b3d1c663a..475fc7c8fd 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H @@ -211,6 +211,12 @@ public: const label i ) const; + //- Return non const access to chemical source terms [kg/m3/s] + virtual DimensionedField& RR + ( + const label i + ); + //- Solve the reaction system for the given start time and time // step and return the characteristic time virtual scalar solve(const scalar t0, const scalar deltaT); diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H index d128443fe9..c009a11e46 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,5 +78,15 @@ Foam::chemistryModel::RR return RR_[i]; } +template +Foam::DimensionedField& +Foam::chemistryModel::RR +( + const label i +) +{ + return RR_[i]; +} + // ************************************************************************* // diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C index cfff71b9ed..f94db5946b 100644 --- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C +++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C @@ -53,4 +53,35 @@ Foam::basicSolidChemistryModel::~basicSolidChemistryModel() {} +const Foam::DimensionedField& +Foam::basicSolidChemistryModel::RR(const label i) const +{ + notImplemented + ( + "const Foam::DimensionedField&" + "basicSolidChemistryModel::RR(const label)" + ); + return (DimensionedField::null()); +} + + +Foam::DimensionedField& +Foam::basicSolidChemistryModel::RR(const label i) +{ + notImplemented + ( + "Foam::DimensionedField&" + "basicSolidChemistryModel::RR(const label)" + ); + + return dynamic_cast&> + ( + const_cast& > + ( + DimensionedField::null() + ) + ); +} + + // ************************************************************************* // diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H index 2134bffca1..50d6a2e96b 100644 --- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H @@ -149,6 +149,15 @@ public: //- Calculates the reaction rates virtual void calculate() = 0; + + //- Return const access to the total source terms + virtual const DimensionedField& RR + ( + const label i + ) const; + + //- Return non-const access to the total source terms + virtual DimensionedField& RR(const label i); }; diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C index e301a1b85e..4962bf716b 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C @@ -535,14 +535,21 @@ updateConcsInReactionI c[si] = max(0.0, c[si]); } + scalar sr = 0.0; forAll(R.rhs(), s) { label si = R.rhs()[s].index; const scalar rhoR = this->solidThermo_[si].rho(p, T); - const scalar sr = rhoR/rhoL; + sr = rhoR/rhoL; c[si] += dt*sr*omeg; c[si] = max(0.0, c[si]); } + + forAll(R.grhs(), g) + { + label gi = R.grhs()[g].index; + c[gi + this->nSolids_] += (1.0 - sr)*omeg*dt; + } } @@ -561,24 +568,11 @@ updateRRInReactionI simpleMatrix& RR ) const { - const Reaction& R = this->reactions_[index]; - scalar rhoL = 0.0; - forAll(R.lhs(), s) - { - label si = R.lhs()[s].index; - rhoL = this->solidThermo_[si].rho(p, T); - RR[si][rRef] -= pr*corr; - RR[si][lRef] += pf*corr; - } - - forAll(R.rhs(), s) - { - label si = R.rhs()[s].index; - const scalar rhoR = this->solidThermo_[si].rho(p, T); - const scalar sr = rhoR/rhoL; - RR[si][lRef] -= sr*pf*corr; - RR[si][rRef] += sr*pr*corr; - } + notImplemented + ( + "void Foam::pyrolysisChemistryModel::" + "updateRRInReactionI" + ); } @@ -666,7 +660,6 @@ Foam::pyrolysisChemistryModel::solve scalar newCp = 0.0; scalar newhi = 0.0; - scalar invRho = 0.0; scalarList dcdt = (c - c0)/dt; for (label i=0; inSolids_; i++) @@ -675,7 +668,6 @@ Foam::pyrolysisChemistryModel::solve scalar Yi = c[i]/cTot; newCp += Yi*this->solidThermo_[i].Cp(pi, Ti); newhi -= dYi*this->solidThermo_[i].Hc(); - invRho += Yi/this->solidThermo_[i].rho(pi, Ti); } scalar dTi = (newhi/newCp)*dt; diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H index 2d4a11cc90..7ae76c5233 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H @@ -137,8 +137,9 @@ public: const bool updateC0 = false ) const; - //- Return the reaction rate for reaction r and the reference - // species and charateristic times + //- Return the reaction rate for reaction r + // NOTE: Currently does not calculate reference specie and + // characteristic times (pf, cf,l Ref, etc.) virtual scalar omega ( const Reaction& r, @@ -153,8 +154,10 @@ public: label& rRef ) const; - //- Return the reaction rate for iReaction and the reference - // species and charateristic times + + //- Return the reaction rate for iReaction + // NOTE: Currently does not calculate reference specie and + // characteristic times (pf, cf,l Ref, etc.) virtual scalar omegaI ( label iReaction, @@ -169,6 +172,7 @@ public: label& rRef ) const; + //- Calculates the reaction rates virtual void calculate(); @@ -186,6 +190,7 @@ public: //- Update matrix RR for reaction i. Used by EulerImplicit + // (Not implemented) virtual void updateRRInReactionI ( const label i, diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H index 03076d79ce..bf67f1ca4d 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -65,6 +65,9 @@ class solidChemistryModel { // Private Member Functions + //- Disallow copy constructor + solidChemistryModel(const solidChemistryModel&); + //- Disallow default bitwise assignment void operator=(const solidChemistryModel&); @@ -151,6 +154,7 @@ public: label& rRef ) const = 0; + //- Return the reaction rate for iReaction and the reference // species and charateristic times virtual scalar omegaI @@ -167,6 +171,10 @@ public: label& rRef ) const = 0; + //- Calculates the reaction rates + virtual void calculate() = 0; + + //- Update concentrations in reaction i given dt and reaction rate // omega used by sequential solver virtual void updateConcsInReactionI @@ -194,11 +202,8 @@ public: simpleMatrix& RR ) const = 0; - //- Calculates the reaction rates - virtual void calculate() = 0; - - // Chemistry model functions + // Solid Chemistry model functions //- Return const access to the chemical source terms for solids inline const DimensionedField& RRs @@ -209,13 +214,6 @@ public: //- Return total solid source term inline tmp > RRs() const; - //- Return const access to the total source terms - inline const DimensionedField& RR - ( - const label i - ) const; - - //- Solve the reaction system for the given start time and time // step and return the characteristic time virtual scalar solve(const scalar t0, const scalar deltaT) = 0; diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H index 92b461e11e..f2d0dad79a 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -95,16 +95,4 @@ Foam::solidChemistryModel::RRs() const } -template -inline const Foam::DimensionedField& -Foam::solidChemistryModel::RR -( - const label i -) const -{ - notImplemented("solidChemistryModel::RR(const label)"); - return (DimensionedField::null()); -} - - // ************************************************************************* // diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H index 12f842a054..447020e2f7 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H @@ -52,9 +52,8 @@ namespace Foam defineTemplateTypeNameAndDebugWithName \ ( \ SS##Schem##Comp##SThermo##GThermo, \ - (#SS"<" + word(Schem::typeName_()) \ - + "<"#Comp"," + SThermo::typeName() \ - + "," + GThermo::typeName() + ">>").c_str(), \ + (#SS"<"#Schem"<"#Comp"," + SThermo::typeName() + "," \ + + GThermo::typeName() + ">>").c_str(), \ 0 \ ); \ \ @@ -77,14 +76,6 @@ namespace Foam GThermo \ ); \ \ - makeSolidChemistrySolverType \ - ( \ - EulerImplicit, \ - SolidChem, \ - Comp, \ - SThermo, \ - GThermo \ - ); \ \ makeSolidChemistrySolverType \ ( \ @@ -104,7 +95,6 @@ namespace Foam GThermo \ ); - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H index 5193a4b9de..07b6ffb050 100644 --- a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H +++ b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H @@ -30,7 +30,7 @@ Description \*---------------------------------------------------------------------------*/ #ifndef makeSolidThermo_H -#define makesolidThermo_H +#define makeSolidThermo_H #include "addToRunTimeSelectionTable.H" From 7531ccba9be86308a9e69a84b784bad8737591be Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 29 Aug 2013 15:32:43 +0100 Subject: [PATCH 4/9] ENH: foamyHexMesh: add rayShooting initial points method --- .../conformalVoronoiMesh/Make/files | 1 + .../conformalVoronoiMesh.C | 2 +- .../rayShooting/rayShooting.C | 218 ++++++++++++++++++ .../rayShooting/rayShooting.H | 105 +++++++++ 4 files changed, 325 insertions(+), 1 deletion(-) create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files index e42c52451f..bc4a363517 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files @@ -64,6 +64,7 @@ initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C initialPointsMethod/faceCentredCubic/faceCentredCubic.C initialPointsMethod/pointFile/pointFile.C initialPointsMethod/autoDensity/autoDensity.C +initialPointsMethod/rayShooting/rayShooting.C relaxationModel/relaxationModel/relaxationModel.C relaxationModel/adaptiveLinear/adaptiveLinear.C diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index 6f137249c1..e580febdd8 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -1249,7 +1249,7 @@ void Foam::conformalVoronoiMesh::move() if ( ( - (vA->internalPoint() && vB->internalPoint()) + (vA->internalPoint() || vB->internalPoint()) && (!vA->referred() || !vB->referred()) // || // ( diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C new file mode 100644 index 0000000000..c4beb20b8e --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C @@ -0,0 +1,218 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "rayShooting.H" +#include "addToRunTimeSelectionTable.H" +#include "triSurfaceMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(rayShooting, 0); +addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +rayShooting::rayShooting +( + const dictionary& initialPointsDict, + const Time& runTime, + Random& rndGen, + const conformationSurfaces& geometryToConformTo, + const cellShapeControl& cellShapeControls, + const autoPtr& decomposition +) +: + initialPointsMethod + ( + typeName, + initialPointsDict, + runTime, + rndGen, + geometryToConformTo, + cellShapeControls, + decomposition + ), + maxRayLength_(readScalar(detailsDict().lookup("maxRayLength"))), + randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")), + randomPerturbationCoeff_ + ( + readScalar(detailsDict().lookup("randomPerturbationCoeff")) + ) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +List rayShooting::initialPoints() const +{ + // Loop over surface faces + const searchableSurfaces& surfaces = geometryToConformTo().geometry(); + const labelList& surfacesToConformTo = geometryToConformTo().surfaces(); + + // Initialise points list + label initialPointsSize = 0; + forAll(surfaces, surfI) + { + initialPointsSize += surfaces[surfI].size(); + } + + DynamicList initialPoints(initialPointsSize); + + forAll(surfacesToConformTo, surfI) + { + const searchableSurface& s = surfaces[surfacesToConformTo[surfI]]; + + tmp faceCentresTmp(s.coordinates()); + const pointField& faceCentres = faceCentresTmp(); + + Info<< " Shoot rays from " << s.name() << nl + << " nRays = " << faceCentres.size() << endl; + + + forAll(faceCentres, fcI) + { + const Foam::point& fC = faceCentres[fcI]; + + if + ( + Pstream::parRun() + && !decomposition().positionOnThisProcessor(fC) + ) + { + continue; + } + + const scalar pert = + randomPerturbationCoeff_ + *cellShapeControls().cellSize(fC); + + pointIndexHit surfHitStart; + label hitSurfaceStart; + + // Face centres should be on the surface so search distance can be + // small + geometryToConformTo().findSurfaceNearest + ( + fC, + sqr(SMALL), + surfHitStart, + hitSurfaceStart + ); + + vectorField normStart(1, vector::min); + geometryToConformTo().getNormal + ( + hitSurfaceStart, + List(1, surfHitStart), + normStart + ); + + pointIndexHit surfHitEnd; + label hitSurfaceEnd; + + geometryToConformTo().findSurfaceNearestIntersection + ( + fC - normStart[0]*SMALL, + fC - normStart[0]*maxRayLength_, + surfHitEnd, + hitSurfaceEnd + ); + + vectorField normEnd(1, vector::min); + geometryToConformTo().getNormal + ( + hitSurfaceEnd, + List(1, surfHitEnd), + normEnd + ); + + if + ( + surfHitEnd.hit() + && ((normStart[0] & normEnd[0]) < 0) + ) + { + line l(fC, surfHitEnd.hitPoint()); + + if (Pstream::parRun()) + { + // Clip the line in parallel + pointIndexHit procIntersection = + decomposition().findLine + ( + l.start() + l.vec()*SMALL, + l.end() - l.vec()*maxRayLength_ + ); + + if (procIntersection.hit()) + { + l = + line + ( + l.start(), + procIntersection.hitPoint() + ); + } + } + + Foam::point midPoint(l.centre()); + + const scalar minDistFromSurfaceSqr = + minimumSurfaceDistanceCoeffSqr_ + *sqr(cellShapeControls().cellSize(midPoint)); + + if + ( + magSqr(midPoint - l.start()) > minDistFromSurfaceSqr + && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr + ) + { + if (randomiseInitialGrid_) + { + midPoint.x() += pert*(rndGen().scalar01() - 0.5); + midPoint.y() += pert*(rndGen().scalar01() - 0.5); + midPoint.z() += pert*(rndGen().scalar01() - 0.5); + } + + initialPoints.append(toPoint(midPoint)); + } + } + } + } + + return initialPoints.shrink(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H new file mode 100644 index 0000000000..c199a31e0a --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::rayShooting + +Description + +SourceFiles + rayShooting.C + +\*---------------------------------------------------------------------------*/ + +#ifndef rayShooting_H +#define rayShooting_H + +#include "initialPointsMethod.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class rayShooting Declaration +\*---------------------------------------------------------------------------*/ + +class rayShooting +: + public initialPointsMethod +{ + +private: + + // Private data + + const scalar maxRayLength_; + + //- Should the initial positions be randomised + Switch randomiseInitialGrid_; + + //- Randomise the initial positions by fraction of the initialCellSize_ + scalar randomPerturbationCoeff_; + + +public: + + //- Runtime type information + TypeName("rayShooting"); + + // Constructors + + //- Construct from components + rayShooting + ( + const dictionary& initialPointsDict, + const Time& runTime, + Random& rndGen, + const conformationSurfaces& geometryToConformTo, + const cellShapeControl& cellShapeControls, + const autoPtr& decomposition + ); + + + //- Destructor + virtual ~rayShooting() + {} + + + // Member Functions + + //- Return the initial points for the conformalVoronoiMesh + virtual List initialPoints() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 8dda57ddbfb9dbccafa9494be3a37d92bbd1e066 Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 29 Aug 2013 16:00:04 +0100 Subject: [PATCH 5/9] ENH: foamyHexMesh: correct rayShooting point addition on intersection --- .../rayShooting/rayShooting.C | 103 +++++++++--------- 1 file changed, 51 insertions(+), 52 deletions(-) diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C index c4beb20b8e..e0be01bcac 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C @@ -121,7 +121,7 @@ List rayShooting::initialPoints() const geometryToConformTo().findSurfaceNearest ( fC, - sqr(SMALL), + sqr(pert), surfHitStart, hitSurfaceStart ); @@ -145,63 +145,62 @@ List rayShooting::initialPoints() const hitSurfaceEnd ); - vectorField normEnd(1, vector::min); - geometryToConformTo().getNormal - ( - hitSurfaceEnd, - List(1, surfHitEnd), - normEnd - ); - - if - ( - surfHitEnd.hit() - && ((normStart[0] & normEnd[0]) < 0) - ) + if (surfHitEnd.hit()) { - line l(fC, surfHitEnd.hitPoint()); - - if (Pstream::parRun()) - { - // Clip the line in parallel - pointIndexHit procIntersection = - decomposition().findLine - ( - l.start() + l.vec()*SMALL, - l.end() - l.vec()*maxRayLength_ - ); - - if (procIntersection.hit()) - { - l = - line - ( - l.start(), - procIntersection.hitPoint() - ); - } - } - - Foam::point midPoint(l.centre()); - - const scalar minDistFromSurfaceSqr = - minimumSurfaceDistanceCoeffSqr_ - *sqr(cellShapeControls().cellSize(midPoint)); - - if + vectorField normEnd(1, vector::min); + geometryToConformTo().getNormal ( - magSqr(midPoint - l.start()) > minDistFromSurfaceSqr - && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr - ) + hitSurfaceEnd, + List(1, surfHitEnd), + normEnd + ); + + if ((normStart[0] & normEnd[0]) < 0) { - if (randomiseInitialGrid_) + line l(fC, surfHitEnd.hitPoint()); + + if (Pstream::parRun()) { - midPoint.x() += pert*(rndGen().scalar01() - 0.5); - midPoint.y() += pert*(rndGen().scalar01() - 0.5); - midPoint.z() += pert*(rndGen().scalar01() - 0.5); + // Clip the line in parallel + pointIndexHit procIntersection = + decomposition().findLine + ( + l.start() + l.vec()*SMALL, + l.end() - l.vec()*maxRayLength_ + ); + + if (procIntersection.hit()) + { + l = + line + ( + l.start(), + procIntersection.hitPoint() + ); + } } - initialPoints.append(toPoint(midPoint)); + Foam::point midPoint(l.centre()); + + const scalar minDistFromSurfaceSqr = + minimumSurfaceDistanceCoeffSqr_ + *sqr(cellShapeControls().cellSize(midPoint)); + + if + ( + magSqr(midPoint - l.start()) > minDistFromSurfaceSqr + && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr + ) + { + if (randomiseInitialGrid_) + { + midPoint.x() += pert*(rndGen().scalar01() - 0.5); + midPoint.y() += pert*(rndGen().scalar01() - 0.5); + midPoint.z() += pert*(rndGen().scalar01() - 0.5); + } + + initialPoints.append(toPoint(midPoint)); + } } } } From 20a8649a0094f4377f767829c1e33a263a205437 Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 29 Aug 2013 16:04:24 +0100 Subject: [PATCH 6/9] ENH: foamyHexMesh: calc maxRayLength from bounds of surface --- .../initialPointsMethod/rayShooting/rayShooting.C | 7 ++++--- .../initialPointsMethod/rayShooting/rayShooting.H | 2 -- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C index e0be01bcac..880871165e 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C @@ -59,7 +59,6 @@ rayShooting::rayShooting cellShapeControls, decomposition ), - maxRayLength_(readScalar(detailsDict().lookup("maxRayLength"))), randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")), randomPerturbationCoeff_ ( @@ -76,6 +75,8 @@ List rayShooting::initialPoints() const const searchableSurfaces& surfaces = geometryToConformTo().geometry(); const labelList& surfacesToConformTo = geometryToConformTo().surfaces(); + const scalar maxRayLength = surfaces.bounds().mag(); + // Initialise points list label initialPointsSize = 0; forAll(surfaces, surfI) @@ -140,7 +141,7 @@ List rayShooting::initialPoints() const geometryToConformTo().findSurfaceNearestIntersection ( fC - normStart[0]*SMALL, - fC - normStart[0]*maxRayLength_, + fC - normStart[0]*maxRayLength, surfHitEnd, hitSurfaceEnd ); @@ -166,7 +167,7 @@ List rayShooting::initialPoints() const decomposition().findLine ( l.start() + l.vec()*SMALL, - l.end() - l.vec()*maxRayLength_ + l.end() - l.vec()*maxRayLength ); if (procIntersection.hit()) diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H index c199a31e0a..e66df62b7c 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H @@ -54,8 +54,6 @@ private: // Private data - const scalar maxRayLength_; - //- Should the initial positions be randomised Switch randomiseInitialGrid_; From 067c80b0d73042f6ec7629c1e561cfbe3bcc6dbd Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 29 Aug 2013 16:34:39 +0100 Subject: [PATCH 7/9] ENH: foamyHexMesh: correct proc boundary intersection in ray shooting --- .../rayShooting/rayShooting.C | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C index 880871165e..ebecad88e1 100644 --- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C +++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C @@ -140,7 +140,7 @@ List rayShooting::initialPoints() const geometryToConformTo().findSurfaceNearestIntersection ( - fC - normStart[0]*SMALL, + fC - normStart[0]*pert, fC - normStart[0]*maxRayLength, surfHitEnd, hitSurfaceEnd @@ -166,8 +166,8 @@ List rayShooting::initialPoints() const pointIndexHit procIntersection = decomposition().findLine ( - l.start() + l.vec()*SMALL, - l.end() - l.vec()*maxRayLength + l.start(), + l.end() ); if (procIntersection.hit()) @@ -187,19 +187,19 @@ List rayShooting::initialPoints() const minimumSurfaceDistanceCoeffSqr_ *sqr(cellShapeControls().cellSize(midPoint)); + if (randomiseInitialGrid_) + { + midPoint.x() += pert*(rndGen().scalar01() - 0.5); + midPoint.y() += pert*(rndGen().scalar01() - 0.5); + midPoint.z() += pert*(rndGen().scalar01() - 0.5); + } + if ( magSqr(midPoint - l.start()) > minDistFromSurfaceSqr && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr ) { - if (randomiseInitialGrid_) - { - midPoint.x() += pert*(rndGen().scalar01() - 0.5); - midPoint.y() += pert*(rndGen().scalar01() - 0.5); - midPoint.z() += pert*(rndGen().scalar01() - 0.5); - } - initialPoints.append(toPoint(midPoint)); } } From f7145eb63ffb26e7c2dda79a2f63c4baac144376 Mon Sep 17 00:00:00 2001 From: Henry Date: Mon, 2 Sep 2013 09:46:25 +0100 Subject: [PATCH 8/9] Removed unnecessary include directory --- .../compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options index 11053f31a9..502938c53c 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options @@ -1,5 +1,4 @@ EXE_INC = \ - -I../rhoPorousMRFPimpleFoam \ -I.. \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ From d730c49e86e6a5aab4e5a2ce3312dc465b311e93 Mon Sep 17 00:00:00 2001 From: Sergio Ferraris Date: Mon, 2 Sep 2013 09:56:39 +0100 Subject: [PATCH 9/9] ENH: Modification of hEqn in reactingOneDim.C introducing laplacian(kappa, T) instead of laplacian(alpha, h). This is necessary when Cp is not uniform. --- .../pyrolysisModels/reactingOneDim/reactingOneDim.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C index 6fa334158a..36ff9f4cad 100644 --- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C +++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C @@ -324,6 +324,8 @@ void reactingOneDim::solveEnergy() ( fvm::ddt(rho_, h_) - fvm::laplacian(alpha, h_) + + fvc::laplacian(alpha, h_) + - fvc::laplacian(kappa(), T()) == chemistrySh_ - fvm::Sp(solidChemistry_->RRg(), h_)