diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/files b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/files new file mode 100644 index 0000000000..677ada2ca2 --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/files @@ -0,0 +1,3 @@ +buoyantBoussinesqSimpleFoam.C + +EXE = $(FOAM_APPBIN)/buoyantBoussinesqSimpleFoam diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options new file mode 100644 index 0000000000..a53061a93c --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options @@ -0,0 +1,12 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/turbulenceModels \ + -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lincompressibleRASModels \ + -lincompressibleTransportModels diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H new file mode 100644 index 0000000000..d00f0877ea --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H @@ -0,0 +1,19 @@ +{ + volScalarField kappaEff + ( + "kappaEff", + turbulence->nu() + turbulence->nut()/Prt + ); + + fvScalarMatrix TEqn + ( + fvm::div(phi, T) + - fvm::Sp(fvc::div(phi), T) + - fvm::laplacian(kappaEff, T) + ); + + TEqn.relax(); + + eqnResidual = TEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); +} diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H new file mode 100644 index 0000000000..cf5c4a1978 --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H @@ -0,0 +1,25 @@ + // Solve the momentum equation + + tmp UEqn + ( + fvm::div(phi, U) + - fvm::Sp(fvc::div(phi), U) + + turbulence->divDevReff(U) + ); + + UEqn().relax(); + + eqnResidual = solve + ( + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + - betaghf*fvc::snGrad(T) + ) * mesh.magSf() + ) + ).initialResidual(); + + maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C new file mode 100644 index 0000000000..5e9deff764 --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 + +Application + buoyantBoussinesqSimpleFoam + +Description + Steady-state solver for buoyant, turbulent flow of incompressible fluids + + Uses the Boussinesq approximation: + \f[ + rho_{eff} = 1 - beta(T - T_{ref}) + \f] + + where: + \f$ rho_{eff} \f$ = the effective (driving) density + beta = thermal expansion coefficient [1/K] + T = temperature [K] + \f$ T_{ref} \f$ = reference temperature [K] + + Valid when: + \f[ + rho_{eff} << 1 + \f] + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "readEnvironmentalProperties.H" +# include "createFields.H" +# include "initContinuityErrs.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + for (runTime++; !runTime.end(); runTime++) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + +# include "readSIMPLEControls.H" +# include "initConvergenceCheck.H" + + pd.storePrevIter(); + + // Pressure-velocity SIMPLE corrector + { +# include "UEqn.H" +# include "TEqn.H" +# include "pdEqn.H" + } + + turbulence->correct(); + + if (runTime.write()) + { +# include "writeAdditionalFields.H" + } + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + +# include "convergenceCheck.H" + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H new file mode 100644 index 0000000000..8958063193 --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H @@ -0,0 +1,9 @@ +// check convergence + +if (maxResidual < convergenceCriterion) +{ + Info<< "reached convergence criterion: " << convergenceCriterion << endl; + runTime.writeAndEnd(); + Info<< "latestTime = " << runTime.timeName() << endl; +} + diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H new file mode 100644 index 0000000000..5f3f13626d --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H @@ -0,0 +1,67 @@ + Info<< "Reading thermophysical properties\n" << endl; + + Info<< "Reading field T\n" << endl; + volScalarField T + ( + IOobject + ( + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // kinematic pd + Info<< "Reading field pd\n" << endl; + volScalarField pd + ( + IOobject + ( + "pd", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +# include "createPhi.H" + +# include "readTransportProperties.H" + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + incompressible::RASModel::New(U, phi, laminarTransport) + ); + + Info<< "Calculating field beta*(g.h)\n" << endl; + surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf())); + + label pdRefCell = 0; + scalar pdRefValue = 0.0; + setRefCell + ( + pd, + mesh.solutionDict().subDict("SIMPLE"), + pdRefCell, + pdRefValue + ); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H new file mode 100644 index 0000000000..b56197f22a --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H @@ -0,0 +1,7 @@ +// initialize values for convergence checks + + scalar eqnResidual = 1, maxResidual = 0; + scalar convergenceCriterion = 0; + + simple.readIfPresent("convergence", convergenceCriterion); + diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H new file mode 100644 index 0000000000..e782d26c8d --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H @@ -0,0 +1,49 @@ +{ + volScalarField rUA("rUA", 1.0/UEqn().A()); + surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA)); + + U = rUA*UEqn().H(); + UEqn.clear(); + + phi = fvc::interpolate(U) & mesh.Sf(); + adjustPhi(phi, U, pd); + surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf(); + phi -= buoyancyPhi; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pdEqn + ( + fvm::laplacian(rUAf, pd) == fvc::div(phi) + ); + + pdEqn.setReference(pdRefCell, pdRefValue); + + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pdEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pdEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + // Calculate the conservative fluxes + phi -= pdEqn.flux(); + + // Explicitly relax pressure for momentum corrector + pd.relax(); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rUAf); + U.correctBoundaryConditions(); + } + } + + #include "continuityErrs.H" +} diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H new file mode 100644 index 0000000000..585128dfde --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H @@ -0,0 +1,13 @@ + singlePhaseTransportModel laminarTransport(U, phi); + + // thermal expansion coefficient [1/K] + dimensionedScalar beta(laminarTransport.lookup("beta")); + + // reference temperature [K] + dimensionedScalar TRef(laminarTransport.lookup("TRef")); + + // reference kinematic pressure [m2/s2] + dimensionedScalar pRef(laminarTransport.lookup("pRef")); + + // turbulent Prandtl number + dimensionedScalar Prt(laminarTransport.lookup("Prt")); diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H new file mode 100644 index 0000000000..20f7c6ae1d --- /dev/null +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H @@ -0,0 +1,29 @@ +{ + volScalarField rhoEff + ( + IOobject + ( + "rhoEff", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + 1.0 - beta*(T - TRef) + ); + rhoEff.write(); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + pd + rhoEff*(g & mesh.C()) + pRef + ); + p.write(); +} diff --git a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H index 988268c4ef..24cfbf4f68 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H @@ -9,4 +9,18 @@ UEqn.relax(); - solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh); + if (momentumPredictor) + { + solve + ( + UEqn + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C index cafb677900..75eb401d3d 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C +++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C @@ -41,37 +41,41 @@ Description int main(int argc, char *argv[]) { + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readEnvironmentalProperties.H" + #include "createFields.H" + #include "initContinuityErrs.H" + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "readEnvironmentalProperties.H" -# include "createFields.H" -# include "initContinuityErrs.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { -# include "readPISOControls.H" -# include "compressibleCourantNo.H" -//# include "setDeltaT.H" + #include "readTimeControls.H" + #include "readPISOControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" runTime++; + Info<< "Time = " << runTime.timeName() << nl << endl; -# include "rhoEqn.H" + #include "rhoEqn.H" -# include "UEqn.H" + #include "UEqn.H" // --- PISO loop for (int corr=0; corrcorrect(); diff --git a/applications/solvers/heatTransfer/buoyantFoam/createFields.H b/applications/solvers/heatTransfer/buoyantFoam/createFields.H index 3d604d4daf..a647f0a3ef 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H @@ -52,11 +52,13 @@ ) ); - Info<< "Creating field dpdt\n" << endl; - volScalarField dpdt = fvc::ddt(p); + Info<< "Creating field DpDt\n" << endl; + volScalarField DpDt = + fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantFoam/hEqn.H index 008e3b0f56..f1a87c6a79 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/hEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/hEqn.H @@ -5,9 +5,7 @@ + fvm::div(phi, h) - fvm::laplacian(turbulence->alphaEff(), h) == - dpdt - + fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p)) - - p*fvc::div(phi/fvc::interpolate(rho)) + DpDt ); hEqn.relax(); diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H index 36378f7f6f..91e190b4cd 100644 --- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H @@ -1,60 +1,67 @@ -bool closedVolume = pd.needReference(); - -rho = thermo->rho(); - -volScalarField rUA = 1.0/UEqn.A(); -U = rUA*UEqn.H(); - -phi = - fvc::interpolate(rho) - *( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) - ) - - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn + bool closedVolume = pd.needReference(); + + rho = thermo->rho(); + + volScalarField rUA = 1.0/UEqn.A(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); + + U = rUA*UEqn.H(); + + surfaceScalarField phiU ( - fvm::ddt(psi, pd) - + fvc::ddt(psi)*pRef - + fvc::ddt(psi, rho)*gh - + fvc::div(phi) - - fvm::laplacian(rho*rUA, pd) + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ) ); - if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - pdEqn.solve(mesh.solver(pd.name() + "Final")); - } - else - { - pdEqn.solve(mesh.solver(pd.name())); + fvScalarMatrix pdEqn + ( + fvm::ddt(psi, pd) + + fvc::ddt(psi)*pRef + + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + - fvm::laplacian(rhorUAf, pd) + ); + + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pdEqn.solve(mesh.solver(pd.name() + "Final")); + } + else + { + pdEqn.solve(mesh.solver(pd.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi += pdEqn.flux(); + } } - if (nonOrth == nNonOrthCorr) + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + + p == pd + rho*gh + pRef; + DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); + + #include "rhoEqn.H" + #include "compressibleContinuityErrs.H" + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi += pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); + rho = thermo->rho(); } -} -p == pd + rho*gh + pRef; -dpdt = fvc::ddt(p); - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); - rho = thermo->rho(); } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H index 5bc4c4738c..71a57cd2a9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H @@ -11,7 +11,15 @@ eqnResidual = solve ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) ).initialResidual(); maxResidual = max(eqnResidual, maxResidual); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index 184242be81..d26da5cb32 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -53,6 +53,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 8c3621205b..3c257d249d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,55 +1,65 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); - -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); + volScalarField rUA = 1.0/UEqn().A(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) + U = rUA*UEqn().H(); + UEqn.clear(); + + phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); + bool closedVolume = adjustPhi(phi, U, p); + surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + phi -= buoyancyPhi; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); + fvScalarMatrix pdEqn + ( + fvm::laplacian(rhorUAf, pd) == fvc::div(phi) + ); + + pdEqn.setReference(pdRefCell, pdRefValue); + + // retain the residual from the first iteration + if (nonOrth == 0) + { + eqnResidual = pdEqn.solve().initialResidual(); + maxResidual = max(eqnResidual, maxResidual); + } + else + { + pdEqn.solve(); + } + + if (nonOrth == nNonOrthCorr) + { + // Calculate the conservative fluxes + phi -= pdEqn.flux(); + + // Explicitly relax pressure for momentum corrector + pd.relax(); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf); + U.correctBoundaryConditions(); + } } - if (nonOrth == nNonOrthCorr) + #include "continuityErrs.H" + + p == pd + rho*gh + pRef; + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) { - phi -= pdEqn.flux(); + p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) + /fvc::domainIntegrate(thermo->psi()); } -} -#include "continuityErrs.H" + rho = thermo->rho(); + rho.relax(); + Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); pd == p - (rho*gh + pRef); } - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options index 22d41fdb80..89d7787af6 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I../buoyantSimpleFoam \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ -I$(LIB_SRC)/turbulenceModels \ diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H deleted file mode 100644 index 5bc4c4738c..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H +++ /dev/null @@ -1,18 +0,0 @@ - // Solve the Momentum equation - - tmp UEqn - ( - fvm::div(phi, U) - - fvm::Sp(fvc::div(phi), U) - + turbulence->divDevRhoReff(U) - ); - - UEqn().relax(); - - eqnResidual = solve - ( - UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh - ).initialResidual(); - - maxResidual = max(eqnResidual, maxResidual); - diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H index 76e3805bca..62c06ec38d 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H @@ -54,6 +54,7 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("gh", g & mesh.Cf()); dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef")); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H deleted file mode 100644 index 2a713bac62..0000000000 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H +++ /dev/null @@ -1,54 +0,0 @@ -volScalarField rUA = 1.0/UEqn().A(); -U = rUA*UEqn().H(); -UEqn.clear(); -phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); -bool closedVolume = adjustPhi(phi, U, p); -phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf(); - -for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) -{ - fvScalarMatrix pdEqn - ( - fvm::laplacian(rho*rUA, pd) == fvc::div(phi) - ); - - pdEqn.setReference(pdRefCell, pdRefValue); - // retain the residual from the first iteration - if (nonOrth == 0) - { - eqnResidual = pdEqn.solve().initialResidual(); - maxResidual = max(eqnResidual, maxResidual); - } - else - { - pdEqn.solve(); - } - - if (nonOrth == nNonOrthCorr) - { - phi -= pdEqn.flux(); - } -} - -#include "continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -pd.relax(); - -p = pd + rho*gh + pRef; - -U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh); -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(thermo->psi()*p)) - /fvc::domainIntegrate(thermo->psi()); - pd == p - (rho*gh + pRef); -} - -rho = thermo->rho(); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files index 3ee71c0ea6..9d9152930a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files @@ -7,6 +7,8 @@ derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemper derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C +fluid/compressibleCourantNo.C + chtMultiRegionFoam.C EXE = $(FOAM_APPBIN)/chtMultiRegionFoam diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index dad472930c..84c7c11806 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -36,16 +36,10 @@ Description #include "turbulenceModel.H" #include "fixedGradientFvPatchFields.H" #include "regionProperties.H" +#include "compressibleCourantNo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#include "solveContinuityEquation.C" -#include "solveMomentumEquation.C" -#include "compressibleContinuityErrors.C" -#include "solvePressureDifferenceEquation.C" -#include "solveEnthalpyEquation.C" -#include "compressibleCourantNo.C" - int main(int argc, char *argv[]) { @@ -58,7 +52,6 @@ int main(int argc, char *argv[]) # include "createSolidMeshes.H" # include "createFluidFields.H" - # include "createSolidFields.H" # include "initContinuityErrs.H" @@ -89,6 +82,7 @@ int main(int argc, char *argv[]) { Info<< "\nSolving for fluid region " << fluidRegions[i].name() << endl; +# include "setRegionFluidFields.H" # include "readFluidMultiRegionPISOControls.H" # include "solveFluid.H" } @@ -97,6 +91,7 @@ int main(int argc, char *argv[]) { Info<< "\nSolving for solid region " << solidRegions[i].name() << endl; +# include "setRegionSolidFields.H" # include "readSolidMultiRegionPISOControls.H" # include "solveSolid.H" } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H index ba3d099a9f..e719d92433 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H @@ -1,10 +1,25 @@ - tmp UEqn = solveMomentumEquation + // Solve the Momentum equation + tmp UEqn ( - momentumPredictor, - Uf[i], - rhof[i], - phif[i], - pdf[i], - ghf[i], - turb[i] + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turb.divDevRhoReff(U) ); + + UEqn().relax(); + + if (momentumPredictor) + { + solve + ( + UEqn() + == + -fvc::reconstruct + ( + ( + fvc::snGrad(pd) + + ghf*fvc::snGrad(rho) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C deleted file mode 100644 index 9cf6446aa0..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C +++ /dev/null @@ -1,58 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 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 - -Description - Continuity errors for fluid meshes - -\*---------------------------------------------------------------------------*/ - -void compressibleContinuityErrors -( - scalar& cumulativeContErr, - const volScalarField& rho, - const basicThermo& thermo -) -{ - dimensionedScalar totalMass = fvc::domainIntegrate(rho); - - scalar sumLocalContErr = - ( - fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass - ).value(); - - scalar globalContErr = - ( - fvc::domainIntegrate(rho - thermo.rho())/totalMass - ).value(); - - cumulativeContErr += globalContErr; - - const word& regionName = rho.mesh().name(); - - Info<< "time step continuity errors (" << regionName << ")" - << ": sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H new file mode 100644 index 0000000000..046ca5ec37 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H @@ -0,0 +1,21 @@ +{ + dimensionedScalar totalMass = fvc::domainIntegrate(rho); + + scalar sumLocalContErr = + ( + fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass + ).value(); + + scalar globalContErr = + ( + fvc::domainIntegrate(rho - thermo.rho())/totalMass + ).value(); + + cumulativeContErr[i] += globalContErr; + + Info<< "time step continuity errors (" << mesh.name() << ")" + << ": sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr[i] + << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C index 487fd93825..e671a256c9 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,13 +22,12 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - Calculates and outputs the mean and maximum Courant Numbers for the fluid - regions - \*---------------------------------------------------------------------------*/ -scalar compressibleCourantNo +#include "compressibleCourantNo.H" +#include "fvc.H" + +Foam::scalar Foam::compressibleCourantNo ( const fvMesh& mesh, const Time& runTime, diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H similarity index 65% rename from applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C rename to applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H index 2e9cae95ad..4bdfb84fcb 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,34 +23,27 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Description - Solve enthalpy equation + Calculates and outputs the mean and maximum Courant Numbers for the fluid + regions \*---------------------------------------------------------------------------*/ -void solveEnthalpyEquation -( - const volScalarField& rho, - const volScalarField& DpDt, - const surfaceScalarField& phi, - const compressible::turbulenceModel& turb, - basicThermo& thermo -) +#ifndef compressibleCourantNo_H +#define compressibleCourantNo_H + +#include "fvMesh.H" + +namespace Foam { - volScalarField& h = thermo.h(); - - tmp hEqn + scalar compressibleCourantNo ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::laplacian(turb.alphaEff(), h) - == - DpDt + const fvMesh& mesh, + const Time& runTime, + const volScalarField& rho, + const surfaceScalarField& phi ); - hEqn().relax(); - hEqn().solve(); - - thermo.correct(); - - Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) - << endl; } + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H index 68c3621ec1..3ca2f68581 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H @@ -7,8 +7,8 @@ ( fluidRegions[regionI], runTime, - rhof[regionI], - phif[regionI] + rhoFluid[regionI], + phiFluid[regionI] ), CoNum ); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H index 314b9d028a..1f6e50de0a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H @@ -1,15 +1,16 @@ // Initialise fluid field pointer lists - PtrList thermof(fluidRegions.size()); - PtrList rhof(fluidRegions.size()); - PtrList Kf(fluidRegions.size()); - PtrList Uf(fluidRegions.size()); - PtrList phif(fluidRegions.size()); - PtrList turb(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); - PtrList ghf(fluidRegions.size()); - PtrList pdf(fluidRegions.size()); + PtrList thermoFluid(fluidRegions.size()); + PtrList rhoFluid(fluidRegions.size()); + PtrList KFluid(fluidRegions.size()); + PtrList UFluid(fluidRegions.size()); + PtrList phiFluid(fluidRegions.size()); + PtrList turbulence(fluidRegions.size()); + PtrList DpDtFluid(fluidRegions.size()); + PtrList ghFluid(fluidRegions.size()); + PtrList ghfFluid(fluidRegions.size()); + PtrList pdFluid(fluidRegions.size()); - List initialMassf(fluidRegions.size()); + List initialMassFluid(fluidRegions.size()); dimensionedScalar pRef ( @@ -24,8 +25,8 @@ Info<< "*** Reading fluid mesh thermophysical properties for region " << fluidRegions[i].name() << nl << endl; - Info<< " Adding to pdf\n" << endl; - pdf.set + Info<< " Adding to pdFluid\n" << endl; + pdFluid.set ( i, new volScalarField @@ -42,16 +43,15 @@ ) ); - Info<< " Adding to thermof\n" << endl; - - thermof.set + Info<< " Adding to thermoFluid\n" << endl; + thermoFluid.set ( i, basicThermo::New(fluidRegions[i]).ptr() ); - Info<< " Adding to rhof\n" << endl; - rhof.set + Info<< " Adding to rhoFluid\n" << endl; + rhoFluid.set ( i, new volScalarField @@ -64,12 +64,12 @@ IOobject::NO_READ, IOobject::AUTO_WRITE ), - thermof[i].rho() + thermoFluid[i].rho() ) ); - Info<< " Adding to Kf\n" << endl; - Kf.set + Info<< " Adding to KFluid\n" << endl; + KFluid.set ( i, new volScalarField @@ -82,12 +82,12 @@ IOobject::NO_READ, IOobject::NO_WRITE ), - thermof[i].Cp()*thermof[i].alpha() + thermoFluid[i].Cp()*thermoFluid[i].alpha() ) ); - Info<< " Adding to Uf\n" << endl; - Uf.set + Info<< " Adding to UFluid\n" << endl; + UFluid.set ( i, new volVectorField @@ -104,8 +104,8 @@ ) ); - Info<< " Adding to phif\n" << endl; - phif.set + Info<< " Adding to phiFluid\n" << endl; + phiFluid.set ( i, new surfaceScalarField @@ -118,29 +118,29 @@ IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), - linearInterpolate(rhof[i]*Uf[i]) + linearInterpolate(rhoFluid[i]*UFluid[i]) & fluidRegions[i].Sf() ) ); - Info<< " Adding to turb\n" << endl; - turb.set + Info<< " Adding to turbulence\n" << endl; + turbulence.set ( i, autoPtr ( compressible::turbulenceModel::New ( - rhof[i], - Uf[i], - phif[i], - thermof[i] + rhoFluid[i], + UFluid[i], + phiFluid[i], + thermoFluid[i] ) ).ptr() ); - Info<< " Adding to DpDtf\n" << endl; - DpDtf.set + Info<< " Adding to DpDtFluid\n" << endl; + DpDtFluid.set ( i, new volScalarField @@ -150,9 +150,9 @@ surfaceScalarField ( "phiU", - phif[i]/fvc::interpolate(rhof[i]) + phiFluid[i]/fvc::interpolate(rhoFluid[i]) ), - thermof[i].p() + thermoFluid[i].p() ) ) ); @@ -162,8 +162,8 @@ ("environmentalProperties"); dimensionedVector g(environmentalProperties.lookup("g")); - Info<< " Adding to ghf\n" << endl; - ghf.set + Info<< " Adding to ghFluid\n" << endl; + ghFluid.set ( i, new volScalarField @@ -172,12 +172,21 @@ g & fluidRegions[i].C() ) ); + ghfFluid.set + ( + i, + new surfaceScalarField + ( + "ghf", + g & fluidRegions[i].Cf() + ) + ); Info<< " Updating p from pd\n" << endl; - thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef; - thermof[i].correct(); + thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef; + thermoFluid[i].correct(); - initialMassf[i] = fvc::domainIntegrate(rhof[i]).value(); + initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value(); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H index 301ceddfdb..d421649f13 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H @@ -1,9 +1,17 @@ - solveEnthalpyEquation +{ + tmp hEqn ( - rhof[i], - DpDtf[i], - phif[i], - turb[i], - thermof[i] + fvm::ddt(rho, h) + + fvm::div(phi, h) + - fvm::laplacian(turb.alphaEff(), h) + == + DpDt ); + hEqn().relax(); + hEqn().solve(); + thermo.correct(); + + Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) + << endl; +} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H new file mode 100644 index 0000000000..1a7f5a3262 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H @@ -0,0 +1 @@ +List cumulativeContErr(fluidRegions.size(), 0.0); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 2b1f5fceb1..e297989809 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -1,60 +1,75 @@ { - bool closedVolume = false; + bool closedVolume = pd.needReference(); - rhof[i] = thermof[i].rho(); + rho = thermo.rho(); volScalarField rUA = 1.0/UEqn().A(); - Uf[i] = rUA*UEqn().H(); + surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); - phif[i] = - fvc::interpolate(rhof[i]) + U = rUA*UEqn().H(); + + surfaceScalarField phiU + ( + fvc::interpolate(rho) *( - (fvc::interpolate(Uf[i]) & fluidRegions[i].Sf()) - + fvc::ddtPhiCorr(rUA, rhof[i], Uf[i], phif[i]) + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) ) - - fvc::interpolate(rhof[i]*rUA*ghf[i]) - *fvc::snGrad(rhof[i]) - *fluidRegions[i].magSf(); + ); - // Solve pressure difference -# include "pdEqn.H" + phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf(); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pdEqn + ( + fvm::ddt(psi, pd) + + fvc::ddt(psi)*pRef + + fvc::ddt(psi, rho)*gh + + fvc::div(phi) + - fvm::laplacian(rho*rUA, pd) + ); + + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pdEqn.solve(mesh.solver(pd.name() + "Final")); + } + else + { + pdEqn.solve(mesh.solver(pd.name())); + } + + if (nonOrth == nNonOrthCorr) + { + phi += pdEqn.flux(); + } + } + + // Correct velocity field + U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); + U.correctBoundaryConditions(); + + // Update pressure field (including bc) + p == pd + rho*gh + pRef; + DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); // Solve continuity # include "rhoEqn.H" - // Update pressure field (including bc) - thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef; - DpDtf[i] = fvc::DDt - ( - surfaceScalarField("phiU", phif[i]/fvc::interpolate(rhof[i])), - thermof[i].p() - ); - // Update continuity errors - compressibleContinuityErrors(cumulativeContErr, rhof[i], thermof[i]); - - // Correct velocity field - Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]); - Uf[i].correctBoundaryConditions(); +# include "compressibleContinuityErrors.H" // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity if (closedVolume) { - thermof[i].p() += - ( - dimensionedScalar - ( - "massIni", - dimMass, - initialMassf[i] - ) - - fvc::domainIntegrate(thermof[i].psi()*thermof[i].p()) - )/fvc::domainIntegrate(thermof[i].psi()); - pdf[i] == thermof[i].p() - (rhof[i]*ghf[i] + pRef); - rhof[i] = thermof[i].rho(); + p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi); + rho = thermo.rho(); } // Update thermal conductivity - Kf[i] = thermof[i].Cp()*turb[i].alphaEff(); + K = thermoFluid[i].Cp()*turb.alphaEff(); + + // Update pd (including bc) + pd == p - (rho*gh + pRef); } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H deleted file mode 100644 index a393843256..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H +++ /dev/null @@ -1,14 +0,0 @@ - solvePressureDifferenceEquation - ( - corr, - nCorr, - nNonOrthCorr, - closedVolume, - pdf[i], - pRef, - rhof[i], - thermof[i].psi(), - rUA, - ghf[i], - phif[i] - ); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H index a7e2a98c4a..a6b0ac9fe1 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H @@ -1 +1 @@ - solveContinuityEquation(rhof[i], phif[i]); + solve(fvm::ddt(rho) + fvc::div(phi)); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H deleted file mode 100644 index 161eb742e7..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H +++ /dev/null @@ -1,15 +0,0 @@ - if (adjustTimeStep) - { - if (CoNum > SMALL) - { - runTime.setDeltaT - ( - min - ( - maxCo*runTime.deltaT().value()/CoNum, - maxDeltaT - ) - ); - } - } - diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H new file mode 100644 index 0000000000..bc4590a428 --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H @@ -0,0 +1,18 @@ + const fvMesh& mesh = fluidRegions[i]; + + basicThermo& thermo = thermoFluid[i]; + volScalarField& rho = rhoFluid[i]; + volScalarField& K = KFluid[i]; + volVectorField& U = UFluid[i]; + surfaceScalarField phi = phiFluid[i]; + compressible::turbulenceModel& turb = turbulence[i]; + volScalarField& DpDt = DpDtFluid[i]; + const volScalarField& gh = ghFluid[i]; + const surfaceScalarField& ghf = ghfFluid[i]; + volScalarField& pd = pdFluid[i]; + + volScalarField& p = thermo.p(); + const volScalarField& psi = thermo.psi(); + volScalarField& h = thermo.h(); + + const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H index e24771256f..19ec50cac2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H @@ -12,4 +12,4 @@ # include "pEqn.H" } } - turb[i].correct(); + turb.correct(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C deleted file mode 100644 index 6d8843ab42..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C +++ /dev/null @@ -1,73 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 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 - -Description - Solve pressure difference equation - -\*---------------------------------------------------------------------------*/ - -void solvePressureDifferenceEquation -( - const label corr, - const label nCorr, - const label nNonOrthCorr, - bool& closedVolume, - volScalarField& pd, - const dimensionedScalar& pRef, - const volScalarField& rho, - const volScalarField& psi, - const volScalarField& rUA, - const volScalarField& gh, - surfaceScalarField& phi -) -{ - closedVolume = pd.needReference(); - - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pdEqn - ( - fvm::ddt(psi, pd) - + fvc::ddt(psi)*pRef - + fvc::ddt(psi, rho)*gh - + fvc::div(phi) - - fvm::laplacian(rho*rUA, pd) - ); - - //pdEqn.solve(); - if (corr == nCorr-1 && nonOrth == nNonOrthCorr) - { - pdEqn.solve(pd.mesh().solver(pd.name() + "Final")); - } - else - { - pdEqn.solve(pd.mesh().solver(pd.name())); - } - - if (nonOrth == nNonOrthCorr) - { - phi += pdEqn.flux(); - } - } -} diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H new file mode 100644 index 0000000000..04c90d2c4c --- /dev/null +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H @@ -0,0 +1,6 @@ +// fvMesh& mesh = solidRegions[i]; + + volScalarField& rho = rhos[i]; + volScalarField& cp = cps[i]; + volScalarField& K = Ks[i]; + volScalarField& T = Ts[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H index 790d4ec934..5fa731824f 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H @@ -3,10 +3,9 @@ { solve ( - fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i]) + fvm::ddt(rho*cp, T) - fvm::laplacian(K, T) ); } - Info<< "Min/max T:" << min(Ts[i]) << ' ' << max(Ts[i]) - << endl; + Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl; } diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index 3bf24a845e..0b2fb4ebf8 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -6,7 +6,7 @@ phic = min(interface.cAlpha()*phic, max(phic)); surfaceScalarField phir = phic*interface.nHatf(); - for (int gCorr=0; gCorr phiPtr_; public: @@ -133,12 +133,12 @@ public: const surfaceScalarField& phi() const { - return *phiPtr_; + return phiPtr_(); } surfaceScalarField& phi() { - return *phiPtr_; + return phiPtr_(); } }; diff --git a/applications/test/HashSet/hashSetTest.C b/applications/test/HashSet/hashSetTest.C index 9625b2e0dc..21e1276604 100644 --- a/applications/test/HashSet/hashSetTest.C +++ b/applications/test/HashSet/hashSetTest.C @@ -111,6 +111,33 @@ int main(int argc, char *argv[]) Info<< "setD : " << setD << endl; Info<< "setB ^ setC ^ setD : " << (setB ^ setC ^ setD) << endl; + // test operator[] + + Info<< "setD : " << setD << endl; + if (setD[0]) + { + Info<< "setD has 0" << endl; + } + else + { + Info<< "setD has no 0" << endl; + } + + + if (setD[11]) + { + Info<< "setD has 11" << endl; + } + else + { + Info<< "setD has no 0" << endl; + } + + Info<< "setD : " << setD << endl; + + // this doesn't work (yet?) + // setD[12] = true; + return 0; } diff --git a/applications/test/PackedList/PackedListTest.C b/applications/test/PackedList/PackedListTest.C index b00c703a37..f1f58f3abc 100644 --- a/applications/test/PackedList/PackedListTest.C +++ b/applications/test/PackedList/PackedListTest.C @@ -38,7 +38,6 @@ using namespace Foam; int main(int argc, char *argv[]) { - bool changed; Info<< "PackedList max_bits() = " << PackedList<0>::max_bits() << nl; Info<< "\ntest allocation with value\n"; @@ -46,11 +45,50 @@ int main(int argc, char *argv[]) list1.print(Info); Info<< "\ntest assign uniform value\n"; - list1 = 2; + list1 = 3; + list1.print(Info); + + Info<< "\ntest assign uniform value (with overflow)\n"; + list1 = -1; + list1.print(Info); + + Info<< "\ntest assign between references\n"; + list1[2] = 3; + list1[4] = list1[2]; + list1.print(Info); + + Info<< "\ntest assign between references, with chaining\n"; + list1[4] = list1[2] = 1; + list1.print(Info); + + { + const PackedList<3>& constLst = list1; + Info<< "\ntest operator[] const with out-of-range index\n"; + constLst.print(Info); + if (!constLst[20]) + { + Info<< "[20] is false (expected) list size should be unchanged (const)\n"; + } + constLst.print(Info); + + Info<< "\ntest operator[] non-const with out-of-range index\n"; + if (!list1[20]) + { + Info<< "[20] is false (expected) but list was resized?? (non-const)\n"; + } + list1.print(Info); + } + + + Info<< "\ntest operator[] with out-of-range index\n"; + if (!list1[20]) + { + Info<< "[20] is false, as expected\n"; + } list1.print(Info); Info<< "\ntest resize with value (without reallocation)\n"; - list1.resize(6, 3); + list1.resize(8, list1.max_value()); list1.print(Info); Info<< "\ntest set() function\n"; @@ -96,7 +134,7 @@ int main(int argc, char *argv[]) list1.print(Info); Info<< "\ntest setCapacity() operation\n"; - list1.setCapacity(30); + list1.setCapacity(100); list1.print(Info); Info<< "\ntest operator[] assignment\n"; @@ -108,7 +146,15 @@ int main(int argc, char *argv[]) list1.print(Info); Info<< "\ntest setCapacity smaller\n"; - list1.setCapacity(32); + list1.setCapacity(24); + list1.print(Info); + + Info<< "\ntest resize much smaller\n"; + list1.resize(150); + list1.print(Info); + + Info<< "\ntest trim\n"; + list1.trim(); list1.print(Info); // add in some misc values @@ -118,37 +164,54 @@ int main(int argc, char *argv[]) Info<< "\ntest iterator\n"; PackedList<3>::iterator iter = list1.begin(); - Info<< "iterator:" << iter() << "\n"; + Info<< "begin():"; iter.print(Info) << "\n"; - Info<< "\ntest iterator operator=\n"; - changed = (iter = 5); - Info<< "iterator:" << iter() << "\n"; - Info<< "changed:" << changed << "\n"; - changed = (iter = 5); - Info<< "changed:" << changed << "\n"; + iter() = 5; + iter.print(Info); list1.print(Info); + iter = list1[31]; + Info<< "iterator:" << iter() << "\n"; + iter.print(Info); + + Info<< "\ntest get() method\n"; - Info<< "get(10):" << list1.get(10) - << " and list[10]:" << unsigned(list1[10]) << "\n"; + Info<< "get(10):" << list1.get(10) << " and list[10]:" << list1[10] << "\n"; list1.print(Info); Info<< "\ntest iterator indexing\n"; - Info<< "end() "; - list1.end().print(Info) << "\n"; + Info<< "cend() "; + list1.cend().print(Info) << "\n"; - for (iter = list1[31]; iter != list1.end(); ++iter) { - iter.print(Info); + Info<< "\ntest assignment of iterator\n"; + list1.print(Info); + PackedList<3>::iterator cit = list1[25]; + cit.print(Info); + list1.end().print(Info); } - Info<< "\ntest operator[] auto-vivify\n"; - const unsigned int val = list1[45]; - Info<< "list[45]:" << val << "\n"; - list1.print(Info); + for + ( + PackedList<3>::iterator cit = list1[5]; + cit != list1.end(); + ++cit + ) + { + cit.print(Info); + } + +// Info<< "\ntest operator[] auto-vivify\n"; +// const unsigned int val = list1[45]; +// +// Info<< "list[45]:" << val << "\n"; +// list1[45] = list1.max_value(); +// Info<< "list[45]:" << list1[45] << "\n"; +// list1[49] = list1.max_value(); +// list1.print(Info); Info<< "\ntest copy constructor + append\n"; @@ -161,8 +224,15 @@ int main(int argc, char *argv[]) Info<< "\ntest pattern that fills all bits\n"; PackedList<4> list3(8, 8); - list3[list3.size()-2] = 0; - list3[list3.size()-1] = list3.max_value(); + + label pos = list3.size() - 1; + + list3[pos--] = list3.max_value(); + list3[pos--] = 0; + list3[pos--] = list3.max_value(); + list3.print(Info); + + Info<< "removed final value: " << list3.remove() << endl; list3.print(Info); Info<< "\n\nDone.\n"; diff --git a/applications/test/PackedList2/Make/files b/applications/test/PackedList2/Make/files new file mode 100644 index 0000000000..bb1f3085fa --- /dev/null +++ b/applications/test/PackedList2/Make/files @@ -0,0 +1,3 @@ +PackedListTest2.C + +EXE = $(FOAM_USER_APPBIN)/PackedListTest2 diff --git a/applications/test/PackedList2/Make/options b/applications/test/PackedList2/Make/options new file mode 100644 index 0000000000..e69de29bb2 diff --git a/applications/test/PackedList2/PackedListTest2.C b/applications/test/PackedList2/PackedListTest2.C new file mode 100644 index 0000000000..22282058ac --- /dev/null +++ b/applications/test/PackedList2/PackedListTest2.C @@ -0,0 +1,348 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 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 + +Application + +Description + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "boolList.H" +#include "PackedBoolList.H" +#include "HashSet.H" +#include "cpuTime.H" +#include + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +// Main program: + +int main(int argc, char *argv[]) +{ + const label n = 1000000; + const label nIters = 1000; + + unsigned int sum = 0; + + PackedBoolList packed(n, 1); + boolList unpacked(n, true); + std::vector stlVector(n, true); + + labelHashSet emptyHash; + labelHashSet fullHash(1000); + for(label i = 0; i < n; i++) + { + fullHash.insert(i); + } + + cpuTime timer; + + for (label iter = 0; iter < nIters; ++iter) + { + packed.resize(40); + packed.shrink(); + packed.resize(n, 1); + } + Info<< "resize/shrink/resize:" << timer.cpuTimeIncrement() << " s\n\n"; + + // set every other bit on: + Info<< "set every other bit on and count\n"; + packed.storage() = 0xAAAAAAAAu; + + // Count packed + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(packed, i) + { + sum += packed[i]; + } + } + Info<< "Counting brute-force:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Count packed + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + sum += packed.count(); + } + Info<< "Counting via count():" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Dummy addition + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(unpacked, i) + { + sum += i + 1; + } + } + Info<< "Dummy loop:" << timer.cpuTimeIncrement() << " s" << endl; + Info<< " sum " << sum << endl; + + // + // Read + // + + // Read stl + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + for(unsigned int i = 0; i < stlVector.size(); i++) + { + sum += stlVector[i]; + } + } + Info<< "Reading stl:" << timer.cpuTimeIncrement() << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read unpacked + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(unpacked, i) + { + sum += unpacked[i]; + } + } + Info<< "Reading unpacked:" << timer.cpuTimeIncrement() << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read packed + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(packed, i) + { + sum += packed.get(i); + } + } + Info<< "Reading packed using get:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read packed + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(packed, i) + { + sum += packed[i]; + } + } + Info<< "Reading packed using reference:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read via iterator + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + for + ( + PackedBoolList::iterator it = packed.begin(); + it != packed.end(); + ++it + ) + { + sum += it; + } + } + Info<< "Reading packed using iterator:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read via iterator + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + for + ( + PackedBoolList::const_iterator cit = packed.cbegin(); + cit != packed.cend(); + ++cit + ) + { + sum += cit(); + } + } + Info<< "Reading packed using const_iterator():" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read empty hash + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(unpacked, i) + { + sum += emptyHash.found(i); + } + } + Info<< "Reading empty labelHashSet:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // Read full hash + sum = 0; + for (label iter = 0; iter < nIters; ++iter) + { + forAll(unpacked, i) + { + sum += fullHash.found(i); + } + } + Info<< "Reading full labelHashSet:" << timer.cpuTimeIncrement() + << " s" << endl; + Info<< " sum " << sum << endl; + + + // + // Write + // + + // Write stl + for (label iter = 0; iter < nIters; ++iter) + { + for (unsigned int i = 0; i < stlVector.size(); i++) + { + stlVector[i] = true; + } + } + Info<< "Writing stl:" << timer.cpuTimeIncrement() << " s" << endl; + + // Write unpacked + for (label iter = 0; iter < nIters; ++iter) + { + forAll(unpacked, i) + { + unpacked[i] = true; + } + } + Info<< "Writing unpacked:" << timer.cpuTimeIncrement() << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + forAll(packed, i) + { + packed[i] = 1; + } + } + Info<< "Writing packed using reference:" << timer.cpuTimeIncrement() + << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + forAll(packed, i) + { + packed.set(i, 1); + } + } + Info<< "Writing packed using set:" << timer.cpuTimeIncrement() + << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + for + ( + PackedBoolList::iterator it = packed.begin(); + it != packed.end(); + ++it + ) + { + it() = 1; + } + } + Info<< "Writing packed using iterator:" << timer.cpuTimeIncrement() + << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + packed = 0; + } + Info<< "Writing packed uniform 0:" << timer.cpuTimeIncrement() + << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + packed = 1; + } + Info<< "Writing packed uniform 1:" << timer.cpuTimeIncrement() + << " s" << endl; + + + PackedList<3> oddPacked(n, 3); + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + packed = 0; + } + Info<< "Writing packed<3> uniform 0:" << timer.cpuTimeIncrement() + << " s" << endl; + + + // Write packed + for (label iter = 0; iter < nIters; ++iter) + { + packed = 1; + } + Info<< "Writing packed<3> uniform 1:" << timer.cpuTimeIncrement() + << " s" << endl; + + + Info << "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/fileName/fileNameTest.C b/applications/test/fileName/fileNameTest.C index c48e4ba2fc..f5a4a93542 100644 --- a/applications/test/fileName/fileNameTest.C +++ b/applications/test/fileName/fileNameTest.C @@ -31,6 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fileName.H" +#include "SubList.H" #include "IOstreams.H" #include "OSspecific.H" @@ -50,19 +51,55 @@ int main() fileName pathName(wrdList); - Info<< "pathName = " << pathName << endl; - Info<< "pathName.name() = " << pathName.name() << endl; - Info<< "pathName.path() = " << pathName.path() << endl; - Info<< "pathName.ext() = " << pathName.ext() << endl; + Info<< "pathName = " << pathName << nl + << "pathName.name() = " << pathName.name() << nl + << "pathName.path() = " << pathName.path() << nl + << "pathName.ext() = " << pathName.ext() << endl; - Info<< "pathName.components() = " << pathName.components() << endl; - Info<< "pathName.component(2) = " << pathName.component(2) << endl; + Info<< "pathName.components() = " << pathName.components() << nl + << "pathName.component(2) = " << pathName.component(2) << nl + << endl; + // try with different combination + for (label start = 0; start < wrdList.size(); ++start) + { + fileName instance, local; + word name; + + fileName path(SubList(wrdList, wrdList.size()-start, start)); + fileName path2 = "." / path; + + path.IOobjectComponents + ( + instance, + local, + name + ); + + Info<< "IOobjectComponents for " << path << nl + << " instance = " << instance << nl + << " local = " << local << nl + << " name = " << name << endl; + + path2.IOobjectComponents + ( + instance, + local, + name + ); + + Info<< "IOobjectComponents for " << path2 << nl + << " instance = " << instance << nl + << " local = " << local << nl + << " name = " << name << endl; + + } // test findEtcFile Info<< "\n\nfindEtcFile tests:" << nl << " controlDict => " << findEtcFile("controlDict") << nl << " badName => " << findEtcFile("badName") << endl; + Info<< "This should emit a fatal error:" << endl; Info<< " badName(die) => " << findEtcFile("badName", true) << nl << endl; diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index c92418dc80..a780d8aecb 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -464,6 +464,7 @@ int main(int argc, char *argv[]) # include "createTime.H" runTime.functionObjects().off(); # include "createPolyMesh.H" + const word oldInstance = mesh.pointsInstance(); scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])())); scalar angle(readScalar(IStringStream(args.additionalArgs()[1])())); @@ -587,8 +588,12 @@ int main(int argc, char *argv[]) { runTime++; } + else + { + mesh.setInstance(oldInstance); + } - Info << "Writing collapsed mesh to time " << runTime.value() << endl; + Info<< "Writing collapsed mesh to time " << runTime.timeName() << endl; mesh.write(); } diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C index c4e47890d5..17092d5a6e 100644 --- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C +++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C @@ -441,6 +441,7 @@ int main(int argc, char *argv[]) # include "createTime.H" runTime.functionObjects().off(); # include "createPolyMesh.H" + const word oldInstance = mesh.pointsInstance(); scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])())); @@ -502,6 +503,11 @@ int main(int argc, char *argv[]) if (nChanged > 0) { + if (overwrite) + { + mesh.setInstance(oldInstance); + } + Info<< "Writing morphed mesh to time " << runTime.timeName() << endl; mesh.write(); diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C index 67fdc0807e..0011d750e5 100644 --- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C +++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C @@ -334,6 +334,7 @@ int main(int argc, char *argv[]) # include "createTime.H" runTime.functionObjects().off(); # include "createPolyMesh.H" + const word oldInstance = mesh.pointsInstance(); bool overwrite = args.options().found("overwrite"); @@ -553,9 +554,13 @@ int main(int argc, char *argv[]) { runTime++; } + else + { + mesh.setInstance(oldInstance); + } // Write resulting mesh - Info << "Writing modified mesh to time " << runTime.value() << endl; + Info << "Writing modified mesh to time " << runTime.timeName() << endl; mesh.write(); } else if (edgeToPos.size()) @@ -602,9 +607,13 @@ int main(int argc, char *argv[]) { runTime++; } + else + { + mesh.setInstance(oldInstance); + } // Write resulting mesh - Info << "Writing modified mesh to time " << runTime.value() << endl; + Info << "Writing modified mesh to time " << runTime.timeName() << endl; mesh.write(); } else @@ -641,9 +650,13 @@ int main(int argc, char *argv[]) { runTime++; } + else + { + mesh.setInstance(oldInstance); + } // Write resulting mesh - Info << "Writing modified mesh to time " << runTime.value() << endl; + Info << "Writing modified mesh to time " << runTime.timeName() << endl; mesh.write(); } diff --git a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C index 4a45f214f6..364aed2e16 100644 --- a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C +++ b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C @@ -58,6 +58,8 @@ int main(int argc, char *argv[]) # include "createTime.H" runTime.functionObjects().off(); # include "createMesh.H" + const word oldInstance = mesh.pointsInstance(); + pointMesh pMesh(mesh); word cellSetName(args.args()[1]); @@ -177,6 +179,10 @@ int main(int argc, char *argv[]) Pout<< "Refined from " << returnReduce(map().nOldCells(), sumOp