From 1c8fd73a790b7259b2094cb0e707c54cdbd84632 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 2 Aug 2010 15:19:37 +0100 Subject: [PATCH 1/5] ENH: Updates to ODEChemistryModel.[CH] --- .../ODEChemistryModel/ODEChemistryModel.C | 190 +++++++++--------- .../ODEChemistryModel/ODEChemistryModel.H | 2 +- 2 files changed, 91 insertions(+), 101 deletions(-) diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C index a3da083706..e9b07da6b1 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -93,7 +93,8 @@ Foam::ODEChemistryModel::~ODEChemistryModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::scalarField Foam::ODEChemistryModel::omega +Foam::tmp +Foam::ODEChemistryModel::omega ( const scalarField& c, const scalar T, @@ -103,7 +104,8 @@ Foam::scalarField Foam::ODEChemistryModel::omega scalar pf, cf, pr, cr; label lRef, rRef; - scalarField om(nEqns(), 0.0); + tmp tom(new scalarField(nEqns(), 0.0)); + scalarField& om = tom(); forAll(reactions_, i) { @@ -116,20 +118,20 @@ Foam::scalarField Foam::ODEChemistryModel::omega forAll(R.lhs(), s) { - label si = R.lhs()[s].index; - scalar sl = R.lhs()[s].stoichCoeff; + const label si = R.lhs()[s].index; + const scalar sl = R.lhs()[s].stoichCoeff; om[si] -= sl*omegai; } forAll(R.rhs(), s) { - label si = R.rhs()[s].index; - scalar sr = R.rhs()[s].stoichCoeff; + const label si = R.rhs()[s].index; + const scalar sr = R.rhs()[s].stoichCoeff; om[si] += sr*omegai; } } - return om; + return tom; } @@ -154,8 +156,8 @@ Foam::scalar Foam::ODEChemistryModel::omega c2[i] = max(0.0, c[i]); } - scalar kf = R.kf(T, p, c2); - scalar kr = R.kr(kf, T, p, c2); + const scalar kf = R.kf(T, p, c2); + const scalar kr = R.kr(kf, T, p, c2); pf = 1.0; pr = 1.0; @@ -169,26 +171,26 @@ Foam::scalar Foam::ODEChemistryModel::omega pf = kf; for (label s=1; s SMALL) { @@ -212,25 +214,25 @@ Foam::scalar Foam::ODEChemistryModel::omega pr = kr; for (label s=1; sSMALL) { @@ -259,8 +261,8 @@ void Foam::ODEChemistryModel::derivatives scalarField& dcdt ) const { - scalar T = c[nSpecie_]; - scalar p = c[nSpecie_ + 1]; + const scalar T = c[nSpecie_]; + const scalar p = c[nSpecie_ + 1]; dcdt = omega(c, T, p); @@ -270,16 +272,16 @@ void Foam::ODEChemistryModel::derivatives scalar cSum = 0.0; for (label i=0; i::derivatives scalar dT = 0.0; for (label i=0; i::jacobian scalarSquareMatrix& dfdc ) const { - scalar T = c[nSpecie_]; - scalar p = c[nSpecie_ + 1]; + const scalar T = c[nSpecie_]; + const scalar p = c[nSpecie_ + 1]; scalarField c2(nSpecie_, 0.0); - for (label i=0; i::jacobian { const Reaction& R = reactions_[ri]; - scalar kf0 = R.kf(T, p, c2); - scalar kr0 = R.kr(T, p, c2); + const scalar kf0 = R.kf(T, p, c2); + const scalar kr0 = R.kr(T, p, c2); forAll(R.lhs(), j) { - label sj = R.lhs()[j].index; + const label sj = R.lhs()[j].index; scalar kf = kf0; forAll(R.lhs(), i) { - label si = R.lhs()[i].index; - scalar el = R.lhs()[i].exponent; + const label si = R.lhs()[i].index; + const scalar el = R.lhs()[i].exponent; if (i == j) { if (el < 1.0) { - if (c2[si]>SMALL) + if (c2[si] > SMALL) { kf *= el*pow(c2[si] + VSMALL, el - 1.0); } @@ -372,31 +374,31 @@ void Foam::ODEChemistryModel::jacobian forAll(R.lhs(), i) { - label si = R.lhs()[i].index; - scalar sl = R.lhs()[i].stoichCoeff; + const label si = R.lhs()[i].index; + const scalar sl = R.lhs()[i].stoichCoeff; dfdc[si][sj] -= sl*kf; } forAll(R.rhs(), i) { - label si = R.rhs()[i].index; - scalar sr = R.rhs()[i].stoichCoeff; + const label si = R.rhs()[i].index; + const scalar sr = R.rhs()[i].stoichCoeff; dfdc[si][sj] += sr*kf; } } forAll(R.rhs(), j) { - label sj = R.rhs()[j].index; + const label sj = R.rhs()[j].index; scalar kr = kr0; forAll(R.rhs(), i) { - label si = R.rhs()[i].index; - scalar er = R.rhs()[i].exponent; - if (i==j) + const label si = R.rhs()[i].index; + const scalar er = R.rhs()[i].exponent; + if (i == j) { - if (er<1.0) + if (er < 1.0) { - if (c2[si]>SMALL) + if (c2[si] > SMALL) { kr *= er*pow(c2[si] + VSMALL, er - 1.0); } @@ -418,25 +420,25 @@ void Foam::ODEChemistryModel::jacobian forAll(R.lhs(), i) { - label si = R.lhs()[i].index; - scalar sl = R.lhs()[i].stoichCoeff; + const label si = R.lhs()[i].index; + const scalar sl = R.lhs()[i].stoichCoeff; dfdc[si][sj] += sl*kr; } forAll(R.rhs(), i) { - label si = R.rhs()[i].index; - scalar sr = R.rhs()[i].stoichCoeff; + const label si = R.rhs()[i].index; + const scalar sr = R.rhs()[i].stoichCoeff; dfdc[si][sj] -= sr*kr; } } } // calculate the dcdT elements numerically - scalar delta = 1.0e-8; - scalarField dcdT0 = omega(c2, T - delta, p); - scalarField dcdT1 = omega(c2, T + delta, p); + const scalar delta = 1.0e-8; + const scalarField dcdT0 = omega(c2, T - delta, p); + const scalarField dcdT1 = omega(c2, T + delta, p); - for (label i=0; i::tc() const scalarField& tc = ttc(); - label nReaction = reactions_.size(); + const label nReaction = reactions_.size(); if (this->chemistry_) { @@ -626,34 +628,31 @@ void Foam::ODEChemistryModel::calculate() this->thermo().rho() ); - for (label i=0; imesh().changing()) { - RR_[i].setSize(rho.size()); + for (label i=0; ichemistry_) { forAll(rho, celli) { + const scalar rhoi = rho[celli]; + const scalar Ti = this->thermo().T()[celli]; + const scalar pi = this->thermo().p()[celli]; + + scalarField c(nSpecie_, 0.0); for (label i=0; ithermo().T()[celli]; - scalar pi = this->thermo().p()[celli]; - - scalarField c(nSpecie_); - scalarField dcdt(nEqns(), 0.0); - - for (label i=0; i::solve const scalar deltaT ) { + scalar deltaTMin = GREAT; + const volScalarField rho ( IOobject @@ -685,35 +686,33 @@ Foam::scalar Foam::ODEChemistryModel::solve this->thermo().rho() ); - for (label i=0; imesh().changing()) { - RR_[i].setSize(rho.size()); + for (label i = 0; i < nSpecie_; i++) + { + RR_[i].setSize(this->mesh().nCells()); + RR_[i] = 0.0; + } } if (!this->chemistry_) { - return GREAT; + return deltaTMin; } - scalar deltaTMin = GREAT; tmp thc = this->thermo().hc(); const scalarField& hc = thc(); forAll(rho, celli) { - for (label i=0; ithermo().hs()[celli] + hc[celli]; + const scalar pi = this->thermo().p()[celli]; scalar Ti = this->thermo().T()[celli]; - scalar hi = this->thermo().hs()[celli] + hc[celli]; - scalar pi = this->thermo().p()[celli]; - scalarField c(nSpecie_); - scalarField c0(nSpecie_); + scalarField c(nSpecie_, 0.0); + scalarField c0(nSpecie_, 0.0); scalarField dc(nSpecie_, 0.0); for (label i=0; i::solve } c0 = c; + // initialise timing parameters scalar t = t0; scalar tauC = this->deltaTChem_[celli]; scalar dt = min(deltaT, tauC); scalar timeLeft = deltaT; // calculate the chemical source terms - scalar cTot = 0.0; - while (timeLeft > SMALL) { tauC = solver().solve(c, Ti, pi, t, dt); t += dt; // update the temperature - cTot = sum(c); + scalar cTot = sum(c); ThermoType mixture(0.0*specieThermo_[0]); for (label i=0; i::solve timeLeft -= dt; this->deltaTChem_[celli] = tauC; - dt = min(timeLeft, tauC); - dt = max(dt, SMALL); + dt = max(SMALL, min(timeLeft, tauC)); } deltaTMin = min(tauC, deltaTMin); dc = c - c0; - scalar WTot = 0.0; - for (label i=0; i& solver() const; //- dc/dt = omega, rate of change in concentration, for each species - virtual scalarField omega + virtual tmp omega ( const scalarField& c, const scalar T, From 195dd9fa543d0a677d0b03dd689c0f675f9d725b Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 6 Aug 2010 12:07:52 +0100 Subject: [PATCH 2/5] ENH: Removed unused 'scale' member data from ode --- .../chemistryModel/chemistrySolver/ode/ode.C | 5 ++--- .../chemistryModel/chemistrySolver/ode/ode.H | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C index e1b4dfad50..35de0dc98c 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C @@ -39,8 +39,7 @@ Foam::ode::ode coeffsDict_(model.subDict(modelName + "Coeffs")), solverName_(coeffsDict_.lookup("ODESolver")), odeSolver_(ODESolver::New(solverName_, model)), - eps_(readScalar(coeffsDict_.lookup("eps"))), - scale_(readScalar(coeffsDict_.lookup("scale"))) + eps_(readScalar(coeffsDict_.lookup("eps"))) {} @@ -67,7 +66,7 @@ Foam::scalar Foam::ode::solve scalarField c1(this->model_.nEqns(), 0.0); // copy the concentration, T and P to the total solve-vector - for (label i=0; i Date: Fri, 6 Aug 2010 12:08:24 +0100 Subject: [PATCH 3/5] ENH: Added 'noChemistrySolver' option/class --- .../chemistrySolver/makeChemistrySolvers.C | 26 +++++ .../noChemistrySolver/noChemistrySolver.C | 66 +++++++++++ .../noChemistrySolver/noChemistrySolver.H | 110 ++++++++++++++++++ 3 files changed, 202 insertions(+) create mode 100644 src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C create mode 100644 src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C index 5f4698fd51..786e90c4b3 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C @@ -29,6 +29,8 @@ License #include "psiChemistryModel.H" #include "rhoChemistryModel.H" +#include "noChemistrySolver.H" + #include "EulerImplicit.H" #include "ode.H" #include "sequential.H" @@ -38,12 +40,24 @@ License namespace Foam { makeChemistrySolver(psiChemistryModel, gasThermoPhysics) + makeChemistrySolverType + ( + noChemistrySolver, + psiChemistryModel, + gasThermoPhysics + ) makeChemistrySolverType(EulerImplicit, psiChemistryModel, gasThermoPhysics) makeChemistrySolverType(ode, psiChemistryModel, gasThermoPhysics) makeChemistrySolverType(sequential, psiChemistryModel, gasThermoPhysics) makeChemistrySolver(psiChemistryModel, icoPoly8ThermoPhysics) makeChemistrySolverType + ( + noChemistrySolver, + psiChemistryModel, + icoPoly8ThermoPhysics + ) + makeChemistrySolverType ( EulerImplicit, psiChemistryModel, @@ -58,12 +72,24 @@ namespace Foam ) makeChemistrySolver(rhoChemistryModel, gasThermoPhysics) + makeChemistrySolverType + ( + noChemistrySolver, + rhoChemistryModel, + gasThermoPhysics + ) makeChemistrySolverType(EulerImplicit, rhoChemistryModel, gasThermoPhysics) makeChemistrySolverType(ode, rhoChemistryModel, gasThermoPhysics) makeChemistrySolverType(sequential, rhoChemistryModel, gasThermoPhysics) makeChemistrySolver(rhoChemistryModel, icoPoly8ThermoPhysics) makeChemistrySolverType + ( + noChemistrySolver, + rhoChemistryModel, + icoPoly8ThermoPhysics + ) + makeChemistrySolverType ( EulerImplicit, rhoChemistryModel, diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C new file mode 100644 index 0000000000..895ca0db2a --- /dev/null +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C @@ -0,0 +1,66 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "noChemistrySolver.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::noChemistrySolver::noChemistrySolver +( + ODEChemistryModel& model, + const word& modelName +) +: + chemistrySolver(model, modelName) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::noChemistrySolver::~noChemistrySolver() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::scalar Foam::noChemistrySolver::solve +( + scalarField&, + const scalar, + const scalar, + const scalar, + const scalar +) const +{ + return GREAT; +} + + +// ************************************************************************* // diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H new file mode 100644 index 0000000000..a844033bc0 --- /dev/null +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H @@ -0,0 +1,110 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::noChemistry + +Description + Dummy chemistry solver for 'none' option + +SourceFiles + noChemistrySolver.H + noChemistrySolver.C + +\*---------------------------------------------------------------------------*/ + +#ifndef noChemistySolver_H +#define noChemistySolver_H + +#include "chemistrySolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +template +class noChemistrySolver; + +/*---------------------------------------------------------------------------*\ + Class noChemistrySolver Declaration +\*---------------------------------------------------------------------------*/ + +template +class noChemistrySolver +: + public chemistrySolver +{ + +public: + + //- Runtime type information + TypeName("none"); + + + // Constructors + + + //- Construct from components + noChemistrySolver + ( + ODEChemistryModel& model, + const word& modelName + ); + + + //- Destructor + virtual ~noChemistrySolver(); + + + // Member Functions + + //- Update the concentrations and return the chemical time + scalar solve + ( + scalarField &c, + const scalar T, + const scalar p, + const scalar t0, + const scalar dt + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "noChemistrySolver.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 02cd5c05da7630431f27c7f850e4677309e10384 Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 6 Aug 2010 12:09:06 +0100 Subject: [PATCH 4/5] STYLE: code clean-up/restructure --- .../ODEChemistryModel/ODEChemistryModel.C | 15 +++---- .../EulerImplicit/EulerImplicit.C | 40 +++++++++---------- .../EulerImplicit/EulerImplicit.H | 7 +++- .../chemistrySolver/sequential/sequential.C | 40 +++++++++---------- .../chemistrySolver/sequential/sequential.H | 8 +++- 5 files changed, 58 insertions(+), 52 deletions(-) diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C index e9b07da6b1..1480fc228d 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C @@ -151,7 +151,7 @@ Foam::scalar Foam::ODEChemistryModel::omega ) const { scalarField c2(nSpecie_, 0.0); - for (label i=0; i::omega lRef = R.lhs()[slRef].index; pf = kf; - for (label s=1; s::omega // find the matrix element and element position for the rhs pr = kr; - for (label s=1; s::derivatives // dT/dt = ... scalar rho = 0.0; scalar cSum = 0.0; - for (label i=0; i::derivatives cp /= mw; scalar dT = 0.0; - for (label i=0; i::Sh() const ) ); + if (this->chemistry_) { scalarField& Sh = tSh(); @@ -560,7 +561,7 @@ Foam::ODEChemistryModel::Sh() const { forAll(Sh, cellI) { - scalar hi = specieThermo_[i].Hc(); + const scalar hi = specieThermo_[i].Hc(); Sh[cellI] -= hi*RR_[i][cellI]; } } @@ -734,7 +735,7 @@ Foam::scalar Foam::ODEChemistryModel::solve t += dt; // update the temperature - scalar cTot = sum(c); + const scalar cTot = sum(c); ThermoType mixture(0.0*specieThermo_[0]); for (label i=0; i::EulerImplicit chemistrySolver(model, modelName), coeffsDict_(model.subDict(modelName + "Coeffs")), cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), - equil_(coeffsDict_.lookup("equilibriumRateLimiter")) + eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")) {} @@ -65,12 +65,12 @@ Foam::scalar Foam::EulerImplicit::solve const label nSpecie = this->model_.nSpecie(); simpleMatrix RR(nSpecie, 0, 0); - for (label i=0; i::solve ); scalar corr = 1.0; - if (equil_) + if (eqRateLimiter_) { - if (omegai<0.0) + if (omegai < 0.0) { corr = 1.0/(1.0 + pr*dt); } @@ -100,31 +100,31 @@ Foam::scalar Foam::EulerImplicit::solve } } - forAll(R.lhs(), s) + forAll(R.lhs(), specieI) { - label si = R.lhs()[s].index; - scalar sl = R.lhs()[s].stoichCoeff; - RR[si][rRef] -= sl*pr*corr; - RR[si][lRef] += sl*pf*corr; + const label id = R.lhs()[specieI].index; + const scalar sc = R.lhs()[specieI].stoichCoeff; + RR[id][rRef] -= sc*pr*corr; + RR[id][lRef] += sc*pf*corr; } - forAll(R.rhs(), s) + forAll(R.rhs(), specieI) { - label si = R.rhs()[s].index; - scalar sr = R.rhs()[s].stoichCoeff; - RR[si][lRef] -= sr*pf*corr; - RR[si][rRef] += sr*pr*corr; + const label id = R.rhs()[specieI].index; + const scalar sc = R.rhs()[specieI].stoichCoeff; + RR[id][lRef] -= sc*pf*corr; + RR[id][rRef] += sc*pr*corr; } } - for (label i=0; i::solve const label nEqns = this->model_.nEqns(); scalarField c1(nEqns, 0.0); - for (label i=0; i::solve scalarField dcdt(nEqns, 0.0); this->model_.derivatives(0.0, c1, dcdt); - scalar sumC = sum(c); + const scalar sumC = sum(c); - for (label i=0; i::sequential chemistrySolver(model, modelName), coeffsDict_(model.subDict(modelName + "Coeffs")), cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), - equil_(coeffsDict_.lookup("equilibriumRateLimiter")) + eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")) {} @@ -70,45 +70,41 @@ Foam::scalar Foam::sequential::solve { const Reaction& R = this->model_.reactions()[i]; - scalar om0 = this->model_.omega + scalar omega = this->model_.omega ( R, c, T, p, pf, cf, lRef, pb, cb, rRef ); - scalar omeg = 0.0; - if (!equil_) + if (eqRateLimiter_) { - omeg = om0; - } - else - { - if (om0<0.0) + if (omega < 0.0) { - omeg = om0/(1.0 + pb*dt); + omega /= 1.0 + pb*dt; } else { - omeg = om0/(1.0 + pf*dt); + omega /= 1.0 + pf*dt; } } - tChemInv = max(tChemInv, mag(omeg)); + + tChemInv = max(tChemInv, mag(omega)); // update species - forAll(R.lhs(), s) + forAll(R.lhs(), specieI) { - label si = R.lhs()[s].index; - scalar sl = R.lhs()[s].stoichCoeff; - c[si] -= dt*sl*omeg; - c[si] = max(0.0, c[si]); + const label id = R.lhs()[specieI].index; + const scalar sc = R.lhs()[specieI].stoichCoeff; + c[id] -= dt*sc*omega; + c[id] = max(0.0, c[id]); } - forAll(R.rhs(), s) + forAll(R.rhs(), specieI) { - label si = R.rhs()[s].index; - scalar sr = R.rhs()[s].stoichCoeff; - c[si] += dt*sr*omeg; - c[si] = max(0.0, c[si]); + const label id = R.rhs()[specieI].index; + const scalar sc = R.rhs()[specieI].stoichCoeff; + c[id] += dt*sc*omega; + c[id] = max(0.0, c[id]); } } diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H index 02e50025c6..32ef0e62bf 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H @@ -59,12 +59,17 @@ class sequential { // Private data + //- Coefficients dictionary dictionary coeffsDict_; + // Model constants + //- Chemistry time scale scalar cTauChem_; - Switch equil_; + + //- Equilibrium rate limiter flag (on/off) + Switch eqRateLimiter_; public: @@ -75,7 +80,6 @@ public: // Constructors - //- Construct from components sequential ( From c91b71b5217d62a12ccf6d727c6e6bc2aebfe0e8 Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 6 Aug 2010 13:20:46 +0100 Subject: [PATCH 5/5] ENH: Added support for DimensionedField to reconstructPar --- .../reconstructPar/reconstructPar.C | 71 +++++------ .../reconstruct/fvFieldReconstructor.C | 3 +- .../reconstruct/fvFieldReconstructor.H | 30 +++-- .../fvFieldReconstructorReconstructFields.C | 111 +++++++++++++++++- .../reconstruct/pointFieldReconstructor.C | 3 +- .../reconstruct/pointFieldReconstructor.H | 9 ++ 6 files changed, 182 insertions(+), 45 deletions(-) diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C index 3734296f51..559bb4a715 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C +++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C @@ -179,23 +179,8 @@ int main(int argc, char *argv[]) // Get list of objects from processor0 database IOobjectList objects(procMeshes.meshes()[0], databases[0].timeName()); - - // If there are any FV fields, reconstruct them - - if - ( - objects.lookupClass(volScalarField::typeName).size() - || objects.lookupClass(volVectorField::typeName).size() - || objects.lookupClass(volSphericalTensorField::typeName).size() - || objects.lookupClass(volSymmTensorField::typeName).size() - || objects.lookupClass(volTensorField::typeName).size() - || objects.lookupClass(surfaceScalarField::typeName).size() - || objects.lookupClass(surfaceVectorField::typeName).size() - || objects.lookupClass(surfaceSphericalTensorField::typeName).size() - || objects.lookupClass(surfaceSymmTensorField::typeName).size() - || objects.lookupClass(surfaceTensorField::typeName).size() - ) { + // If there are any FV fields, reconstruct them Info<< "Reconstructing FV fields" << nl << endl; fvFieldReconstructor fvReconstructor @@ -207,6 +192,32 @@ int main(int argc, char *argv[]) procMeshes.boundaryProcAddressing() ); + fvReconstructor.reconstructFvVolumeInternalFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeInternalFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeInternalFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeInternalFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeInternalFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields ( objects, @@ -258,22 +269,13 @@ int main(int argc, char *argv[]) objects, selectedFields ); - } - else - { - Info<< "No FV fields" << nl << endl; + + if (fvReconstructor.nReconstructed() == 0) + { + Info<< "No FV fields" << nl << endl; + } } - - // If there are any point fields, reconstruct them - if - ( - objects.lookupClass(pointScalarField::typeName).size() - || objects.lookupClass(pointVectorField::typeName).size() - || objects.lookupClass(pointSphericalTensorField::typeName).size() - || objects.lookupClass(pointSymmTensorField::typeName).size() - || objects.lookupClass(pointTensorField::typeName).size() - ) { Info<< "Reconstructing point fields" << nl << endl; @@ -318,10 +320,11 @@ int main(int argc, char *argv[]) objects, selectedFields ); - } - else - { - Info<< "No point fields" << nl << endl; + + if (pointReconstructor.nReconstructed() == 0) + { + Info<< "No point fields" << nl << endl; + } } diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C index 9a8138f5c8..ad92a345aa 100644 --- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C +++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C @@ -40,7 +40,8 @@ Foam::fvFieldReconstructor::fvFieldReconstructor procMeshes_(procMeshes), faceProcAddressing_(faceProcAddressing), cellProcAddressing_(cellProcAddressing), - boundaryProcAddressing_(boundaryProcAddressing) + boundaryProcAddressing_(boundaryProcAddressing), + nReconstructed_(0) {} diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H index f959738fd8..81d40b634d 100644 --- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H +++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H @@ -71,6 +71,9 @@ class fvFieldReconstructor //- List of processor boundary addressing lists const PtrList& boundaryProcAddressing_; + //- Number of fields reconstructed + label nReconstructed_; + // Private Member Functions @@ -134,20 +137,33 @@ public: // Member Functions + //- Return number of fields reconstructed + label nReconstructed() const + { + return nReconstructed_; + } + + //- Reconstruct volume internal field + template + tmp > + reconstructFvVolumeInternalField(const IOobject& fieldIoObject); + //- Reconstruct volume field template tmp > - reconstructFvVolumeField - ( - const IOobject& fieldIoObject - ); + reconstructFvVolumeField(const IOobject& fieldIoObject); //- Reconstruct surface field template tmp > - reconstructFvSurfaceField + reconstructFvSurfaceField(const IOobject& fieldIoObject); + + //- Reconstruct and write all/selected volume internal fields + template + void reconstructFvVolumeInternalFields ( - const IOobject& fieldIoObject + const IOobjectList& objects, + const HashSet& selectedFields ); //- Reconstruct and write all/selected volume fields @@ -158,7 +174,7 @@ public: const HashSet& selectedFields ); - //- Reconstruct and write all/selected volume fields + //- Reconstruct and write all/selected surface fields template void reconstructFvSurfaceFields ( diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C index e01840cd7e..def5b3b1cd 100644 --- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C +++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C @@ -33,6 +33,75 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +Foam::tmp > +Foam::fvFieldReconstructor::reconstructFvVolumeInternalField +( + const IOobject& fieldIoObject +) +{ + // Read the field for all the processors + PtrList > procFields + ( + procMeshes_.size() + ); + + forAll(procMeshes_, procI) + { + procFields.set + ( + procI, + new DimensionedField + ( + IOobject + ( + fieldIoObject.name(), + procMeshes_[procI].time().timeName(), + procMeshes_[procI], + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + procMeshes_[procI] + ) + ); + } + + + // Create the internalField + Field internalField(mesh_.nCells()); + + forAll(procMeshes_, procI) + { + const DimensionedField& procField = procFields[procI]; + + // Set the cell values in the reconstructed field + internalField.rmap + ( + procField.field(), + cellProcAddressing_[procI] + ); + } + + return tmp > + ( + new DimensionedField + ( + IOobject + ( + fieldIoObject.name(), + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + procFields[0].dimensions(), + internalField + ) + ); +} + + template Foam::tmp > Foam::fvFieldReconstructor::reconstructFvVolumeField @@ -451,7 +520,41 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField } -// Reconstruct and write all/selected volume fields +template +void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields +( + const IOobjectList& objects, + const HashSet& selectedFields +) +{ + const word& fieldClassName = DimensionedField::typeName; + + IOobjectList fields = objects.lookupClass(fieldClassName); + + if (fields.size()) + { + Info<< " Reconstructing " << fieldClassName << "s\n" << endl; + + forAllConstIter(IOobjectList, fields, fieldIter) + { + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; + + reconstructFvVolumeInternalField(*fieldIter())().write(); + + nReconstructed_++; + } + } + Info<< endl; + } +} + + template void Foam::fvFieldReconstructor::reconstructFvVolumeFields ( @@ -479,13 +582,15 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields Info<< " " << fieldIter()->name() << endl; reconstructFvVolumeField(*fieldIter())().write(); + + nReconstructed_++; } } Info<< endl; } } -// Reconstruct and write all/selected surface fields + template void Foam::fvFieldReconstructor::reconstructFvSurfaceFields ( @@ -513,6 +618,8 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields Info<< " " << fieldIter()->name() << endl; reconstructFvSurfaceField(*fieldIter())().write(); + + nReconstructed_++; } } Info<< endl; diff --git a/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.C b/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.C index 7de433c8a9..71ef2f2204 100644 --- a/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.C +++ b/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.C @@ -39,7 +39,8 @@ Foam::pointFieldReconstructor::pointFieldReconstructor procMeshes_(procMeshes), pointProcAddressing_(pointProcAddressing), boundaryProcAddressing_(boundaryProcAddressing), - patchPointAddressing_(procMeshes.size()) + patchPointAddressing_(procMeshes.size()), + nReconstructed_(0) { // Inverse-addressing of the patch point labels. labelList pointMap(mesh_.size(), -1); diff --git a/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.H b/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.H index f9e703cc93..faafb1d291 100644 --- a/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.H +++ b/src/parallel/reconstruct/reconstruct/pointFieldReconstructor.H @@ -68,6 +68,9 @@ class pointFieldReconstructor //- Point patch addressing labelListListList patchPointAddressing_; + //- Number of fields reconstructed + label nReconstructed_; + // Private Member Functions @@ -130,6 +133,12 @@ public: // Member Functions + //- Return number of fields reconstructed + label nReconstructed() const + { + return nReconstructed_; + } + //- Reconstruct field template tmp >