diff --git a/applications/solvers/compressible/rhoPimpleFoam/Allwmake b/applications/solvers/compressible/rhoPimpleFoam/Allwmake index 241e22eb1b6c4fd85d2b42836dcc7d78cbfa4915..f6f8ad3635d10f11e2557d9e6f0af3258a738bb2 100755 --- a/applications/solvers/compressible/rhoPimpleFoam/Allwmake +++ b/applications/solvers/compressible/rhoPimpleFoam/Allwmake @@ -3,6 +3,7 @@ cd ${0%/*} || exit 1 # run from this directory set -x wmake +wmake rhoPimplecFoam wmake rhoPorousMRFPimpleFoam wmake rhoPorousMRFLTSPimpleFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H index 1e39e983e42a476d62c73a65a0ef531771f43991..b7b4b04db5c477cefc3dd2dfeab98d3ee6557a11 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H @@ -4,13 +4,12 @@ tmp<fvVectorMatrix> UEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) + - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U) + turbulence->divDevRhoReff(U) ); UEqn().relax(); -volScalarField rAU(1.0/UEqn().A()); - if (pimple.momentumPredictor()) { solve(UEqn() == -fvc::grad(p)); diff --git a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H index b1d9e4e8b0237ca0b9f28f6d36c8f5162121c846..57d72887e60ba1f6dda559ff3353d1294ff3346d 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H @@ -3,10 +3,14 @@ ( fvm::ddt(rho, h) + fvm::div(phi, h) + - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h) - fvm::laplacian(turbulence->alphaEff(), h) == dpdt - - (fvc::ddt(rho, K) + fvc::div(phi, K)) + - ( + fvc::ddt(rho, K) + fvc::div(phi, K) + - (fvc::ddt(rho) + fvc::div(phi))*K + ) ); hEqn.relax(); diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index cfe34f2873146dabbcc4880755867b57e7609766..3d324e085b943574aec0f3ef60f8901e13d8cb0c 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -3,6 +3,7 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); +volScalarField rAU(1.0/UEqn().A()); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn().H(); @@ -23,13 +24,15 @@ if (pimple.transonic()) ) ); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); @@ -52,6 +55,8 @@ else ) ); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { // Pressure corrector @@ -59,7 +64,7 @@ else ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C index 26d11a14aa6694211f3c9dfba3e97825bd7a79e1..28787716e0ac365dbb84a8a130e6cf23fdad0c5b 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -66,7 +66,10 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..5bb1f8d0239f1f82a86fbe88e38b7203ca9f8023 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/files @@ -0,0 +1,3 @@ +rhoPimplecFoam.C + +EXE = $(FOAM_APPBIN)/rhoPimplecFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..d16a0ee1d88a2d9a037f976e6f5b9bc7c7074261 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/Make/options @@ -0,0 +1,15 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lbasicThermophysicalModels \ + -lspecie \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H new file mode 100644 index 0000000000000000000000000000000000000000..7ef1b6253319db33d97a6617e6ba464e34eef840 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H @@ -0,0 +1,117 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn().A()); +volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1())); + +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn().H(); + +if (pimple.nCorrPIMPLE() <= 1) +{ + UEqn.clear(); +} + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + surfaceScalarField phic + ( + "phic", + fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() + ); + + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + + fvc::div(phic) + - fvm::laplacian(Dp, p) + ); + + // Relax the pressure equation to maintain diagonal dominance + pEqn.relax(); + + pEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi == phic + pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho) + *( + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ) + ); + + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(Dp, p) + ); + + pEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +U = HbyA - rAtU*fvc::grad(p); +U.correctBoundaryConditions(); +K = 0.5*magSqr(U); + +dpdt = fvc::ddt(p); + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); + +if (!pimple.transonic()) +{ + rho.relax(); +} + +Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C new file mode 100644 index 0000000000000000000000000000000000000000..0b01f51c8217540f654a7d10d57817f9d69b6fc6 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C @@ -0,0 +1,105 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Application + rhoPimplecFoam + +Description + Transient solver for laminar or turbulent flow of compressible fluids + for HVAC and similar applications. + + Uses the flexible PIMPLEC (PISOC-SIMPLEC) solution for time-resolved and + pseudo-transient simulations. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "basicPsiThermo.H" +#include "turbulenceModel.H" +#include "bound.H" +#include "pimpleControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + pimpleControl pimple(mesh); + + #include "createFields.H" + #include "initContinuityErrs.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "hEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C index de32f1f5c633c48aae0434b92799818183829cf5..e305916f2928e2dd7dcb1a66fc2a6d2ad8f27dc6 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,10 @@ int main(int argc, char *argv[]) #include "setrDeltaT.H" - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H index c63f0f8e4b080cdb6c6662788c176982a6c42a59..c24b1f587a59bac6cd7b305b5dc6ccdf55986613 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H @@ -13,8 +13,6 @@ UEqn().relax(); mrfZones.addCoriolis(rho, UEqn()); pZones.addResistance(UEqn()); -volScalarField rAU(1.0/UEqn().A()); - if (pimple.momentumPredictor()) { solve(UEqn() == -fvc::grad(p)); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H index c5c7602a43ff4bbbf5230305157bcdc3e9b03cb3..3c75e87e6461998422e897f490c5edff9668dfe2 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H @@ -25,13 +25,15 @@ if (pimple.transonic()) ); mrfZones.relativeFlux(fvc::interpolate(psi), phid); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); @@ -56,6 +58,8 @@ else mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA); + volScalarField Dp("Dp", rho*rAU); + while (pimple.correctNonOrthogonal()) { // Pressure corrector @@ -63,7 +67,7 @@ else ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C index 068de89952e1f5b940ffac33f4af63c3c82b20ec..e02787c35d6201c2c7c92635b4177521ce356707 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,7 +69,10 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; - #include "rhoEqn.H" + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + } // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H index 879578cf49d9d34b0e30b56e426a59dcd969cc09..92cea341ca69cd8df156c5ebaaa8696611182499 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H @@ -11,7 +11,6 @@ ); hEqn.relax(); - hEqn.solve(); thermo.correct(); diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index cefee48969acf984baf1c8a8b54211290576ab45..6c64056be2e8e2e11a00c5286515b93e277fa085 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -54,7 +54,8 @@ else { fvScalarMatrix pEqn ( - fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA) + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) ); pEqn.setReference(pRefCell, pRefValue); @@ -63,7 +64,7 @@ else if (simple.finalNonOrthogonalIter()) { - phi = phiHbyA - pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } @@ -88,5 +89,10 @@ if (closedVolume) rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); -rho.relax(); + +if (!simple.transonic()) +{ + rho.relax(); +} + Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H deleted file mode 100644 index 46a382864e72740185c60ca49560ff95f2a2bfdb..0000000000000000000000000000000000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/createFields.H +++ /dev/null @@ -1,61 +0,0 @@ - Info<< "Reading thermophysical properties\n" << endl; - - autoPtr<basicPsiThermo> pThermo - ( - basicPsiThermo::New(mesh) - ); - basicPsiThermo& thermo = pThermo(); - - volScalarField rho - ( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - thermo.rho() - ); - - volScalarField& p = thermo.p(); - volScalarField& h = thermo.h(); - const volScalarField& psi = thermo.psi(); - - Info<< "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - #include "compressibleCreatePhi.H" - - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, simple.dict(), pRefCell, pRefValue); - - dimensionedScalar rhoMax(simple.dict().lookup("rhoMax")); - dimensionedScalar rhoMin(simple.dict().lookup("rhoMin")); - - Info<< "Creating turbulence model\n" << endl; - autoPtr<compressible::RASModel> turbulence - ( - compressible::RASModel::New - ( - rho, - U, - phi, - thermo - ) - ); - - dimensionedScalar initialMass = fvc::domainIntegrate(rho); diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H index 379e720c042eb8d08cdae81ed4064549a84aa7f1..8c7405c346ffdebce86d052fe60be505cb52d584 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H @@ -3,13 +3,11 @@ rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); -volScalarField p0(p); - -volScalarField AU(UEqn().A()); -volScalarField AtU(AU - UEqn().H1()); +volScalarField rAU(1.0/UEqn().A()); +volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1())); volVectorField HbyA("HbyA", U); -HbyA = UEqn().H()/AU; +HbyA = rAU*UEqn().H(); UEqn.clear(); @@ -17,28 +15,32 @@ bool closedVolume = false; if (simple.transonic()) { - while (simple.correctNonOrthogonal()) - { - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi*HbyA) & mesh.Sf() - ); + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); - surfaceScalarField phic - ( - "phic", - fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf() - ); + surfaceScalarField phic + ( + "phic", + fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() + ); + + HbyA -= (rAU - rAtU)*fvc::grad(p); + volScalarField Dp("Dp", rho*rAtU); + + while (simple.correctNonOrthogonal()) + { fvScalarMatrix pEqn ( fvm::div(phid, p) + fvc::div(phic) - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); - // Relax the pressure equation to ensure diagonal-dominance + // Relax the pressure equation to maintain diagonal dominance pEqn.relax(); pEqn.setReference(pRefCell, pRefValue); @@ -53,21 +55,25 @@ if (simple.transonic()) } else { - while (simple.correctNonOrthogonal()) - { - surfaceScalarField phiHbyA - ( - "phiHbyA", - fvc::interpolate(rho*HbyA) & mesh.Sf() - ); + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); - closedVolume = adjustPhi(phi, U, p); - phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf(); + closedVolume = adjustPhi(phiHbyA, U, p); + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField Dp("Dp", rho*rAtU); + + while (simple.correctNonOrthogonal()) + { fvScalarMatrix pEqn ( fvc::div(phiHbyA) - - fvm::laplacian(rho/AtU, p) + - fvm::laplacian(Dp, p) ); pEqn.setReference(pRefCell, pRefValue); @@ -81,16 +87,14 @@ else } } -// The incompressibe for of the continuity error check is appropriate for +// The incompressibe form of the continuity error check is appropriate for // steady-state compressible also. #include "incompressible/continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); -U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU); -//U = HbyA - fvc::grad(p)/AU; - +U = HbyA - rAtU*fvc::grad(p); U.correctBoundaryConditions(); // For closed-volume cases adjust the pressure and density levels @@ -101,6 +105,7 @@ if (closedVolume) /fvc::domainIntegrate(psi); } +// Recalculate density from the relaxed pressure rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H index 46cc36216c45da4dde31afbddf045ae740952fe5..aeebee2df51502346e832f342259848a83934e1a 100644 --- a/applications/solvers/compressible/sonicFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -14,13 +14,16 @@ surfaceScalarField phid ) ); + +volScalarField Dp("Dp", rho*rAU); + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H index cad5bb72170900e8d4fae2c7600468f4284f1ce2..21dc48614ec948ebd3d51c6fb22556aefdf85c5b 100644 --- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H +++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H @@ -14,13 +14,15 @@ surfaceScalarField phid ) ); +volScalarField Dp("Dp", rho*rAU); + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C index 6388caea6952dfac09eebb645de7cd25df0080ab..19a530a6aaa52c9d41198f4d3aa6a42346ea9af6 100644 --- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,13 +85,14 @@ int main(int argc, char *argv[]) ); phi = (rhoO/psi)*phid; + volScalarField Dp("Dp", rho*rAU); fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phi) + fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) + - fvm::laplacian(Dp, p) ); pEqn.solve(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files index 908d7b3a95b5ac8370055a4e221543e030942d85..522f90e3efa41e7539958f671d5fac688179bd4c 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files @@ -2,5 +2,4 @@ fluid/compressibleCourantNo.C solid/solidRegionDiffNo.C chtMultiRegionFoam.C -EXE = $(FOAM_USER_APPBIN)/chtMultiRegionFoam - +EXE = $(FOAM_APPBIN)/chtMultiRegionFoam diff --git a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C index b620536abbb13117c0e6154ae1c0bb0bb89e86b0..bf12294638d314c0b7d96d241afaea4c23e0b2ce 100644 --- a/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C +++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C @@ -61,6 +61,7 @@ int main(int argc, char *argv[]) fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(fluid.nu(), U) + - (fvc::grad(U) & fvc::grad(fluid.nu())) ); solve(UEqn == -fvc::grad(p)); diff --git a/applications/test/PisoFoam/PisoFoam.C b/applications/test/PisoFoam/PisoFoam.C index 67f0084e7192f7024446e7d8994b9a3247b5060a..1713e9747c77b09aa357ff547f9bf34743f40ea4 100644 --- a/applications/test/PisoFoam/PisoFoam.C +++ b/applications/test/PisoFoam/PisoFoam.C @@ -102,7 +102,7 @@ int main(int argc, char *argv[]) "{" " /*solver SmoothSolver;*/" " smoother GaussSeidel;" - " solver PBiCICG;" + " solver PBiCCCG;" " preconditioner none;" " tolerance (1e-7 1e-7 1);" " relTol (0 0 0);" diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files index d89ca6e737c6c84ce878bd4e57566fbd070333a4..a15838abe84da087d0c63693617ed5b410e56b22 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/Make/files +++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files @@ -1,3 +1,4 @@ collapseEdges.C +pointEdgeCollapse/pointEdgeCollapse.C EXE = $(FOAM_APPBIN)/collapseEdges diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/options b/applications/utilities/mesh/advanced/collapseEdges/Make/options index 21b17b14c9d7dd8ae8c0c55fdcf698ce54019fcb..d1efa61fd56aef52623d7e352c300f8948503fd2 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/Make/options +++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options @@ -1,6 +1,10 @@ EXE_INC = \ - -I$(LIB_SRC)/dynamicMesh/lnInclude + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -IpointEdgeCollapse EXE_LIBS = \ -ldynamicMesh \ - -lmeshTools + -lmeshTools \ + -lfiniteVolume diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index f1eb28e5947759fe22aee59bf0c28a3ed3a2c0f0..899159b2a29ba4005f6fdd740260f1fee794674b 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,170 +41,634 @@ Description #include "argList.H" #include "Time.H" -#include "edgeCollapser.H" #include "polyTopoChange.H" #include "polyTopoChanger.H" #include "polyMesh.H" +#include "fvMesh.H" #include "mapPolyMesh.H" #include "mathematicalConstants.H" #include "PackedBoolList.H" #include "SortableList.H" #include "unitConversion.H" +#include "globalMeshData.H" +#include "globalIndex.H" +#include "PointEdgeWave.H" +#include "pointEdgeCollapse.H" +#include "motionSmoother.H" + +#include "OFstream.H" +#include "meshTools.H" using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Get faceEdges in order of face points, i.e. faceEdges[0] is between -// f[0] and f[1] -labelList getSortedEdges +label findIndex ( - const edgeList& edges, - const labelList& f, - const labelList& edgeLabels + const labelList& elems, + const label start, + const label size, + const label val ) { - labelList faceEdges(edgeLabels.size(), -1); + for (label i = start; i < size; i++) + { + if (elems[i] == val) + { + return i; + } + } + return -1; +} - // Find starting pos in f for every edgeLabels - forAll(edgeLabels, i) + +void filterFace +( + const label faceI, + face& f, + const List<pointEdgeCollapse>& allPointInfo, + const Map<DynamicList<label> >& collapseStrings +) +{ + label newFp = 0; + + face oldFace = f; + + forAll(f, fp) { - label edgeI = edgeLabels[i]; + label pointI = f[fp]; - const edge& e = edges[edgeI]; + label collapsePoint = allPointInfo[pointI].collapseIndex(); - label fp = findIndex(f, e[0]); - label fp1 = f.fcIndex(fp); + if (collapseStrings.found(collapsePoint)) + { + collapsePoint = collapseStrings[collapsePoint][0]; + } - if (f[fp1] == e[1]) + if (collapsePoint == -1) { - // EdgeI between fp -> fp1 - faceEdges[fp] = edgeI; + WarningIn + ( + "filterFace" + "(const label, face&, const List<pointEdgeCollapse>&)" + ) + << "Point " << pointI << " was not visited by PointEdgeWave" + << endl; + } + else if (collapsePoint == -2) + { + f[newFp++] = pointI; } else { - // EdgeI between fp-1 -> fp - faceEdges[f.rcIndex(fp)] = edgeI; + if (findIndex(f, 0, newFp, collapsePoint) == -1) + { + f[newFp++] = collapsePoint; + } } } - return faceEdges; + + // Check for pinched face. Tries to correct + // - consecutive duplicate vertex. Removes duplicate vertex. + // - duplicate vertex with one other vertex in between (spike). + // Both of these should not really occur! and should be checked before + // collapsing edges. + + const label size = newFp; + + newFp = 2; + + for (label fp = 2; fp < size; fp++) + { + label fp1 = fp-1; + label fp2 = fp-2; + + label pointI = f[fp]; + + // Search for previous occurrence. + label index = findIndex(f, 0, fp, pointI); + + if (index == fp1) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Removing consecutive duplicate vertex in face " + << f << endl; + // Don't store current pointI + } + else if (index == fp2) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Removing non-consecutive duplicate vertex in face " + << f << endl; + // Don't store current pointI and remove previous + newFp--; + } + else if (index != -1) + { + WarningIn + ( + "Foam::edgeCollapser::filterFace(const label faceI, " + "face& f) const" + ) << "Pinched face " << f << endl; + f[newFp++] = pointI; + } + else + { + f[newFp++] = pointI; + } + } + + f.setSize(newFp); } -// Merges edges which are in straight line. I.e. edge split by point. -label mergeEdges +bool setRefinement ( const polyMesh& mesh, - const scalar maxCos, - edgeCollapser& collapser + polyTopoChange& meshMod, + const List<pointEdgeCollapse>& allPointInfo ) { - const pointField& points = mesh.points(); - const edgeList& edges = mesh.edges(); - const labelListList& pointEdges = mesh.pointEdges(); - const labelList& region = collapser.pointRegion(); - const labelList& master = collapser.pointRegionMaster(); + const cellList& cells = mesh.cells(); + const labelList& faceOwner = mesh.faceOwner(); + const labelList& faceNeighbour = mesh.faceNeighbour(); + const labelListList& pointFaces = mesh.pointFaces(); + const pointZoneMesh& pointZones = mesh.pointZones(); + globalIndex globalStrings(mesh.nPoints()); + + boolList removedPoints(mesh.nPoints(), false); + + // Create strings of edges + Map<DynamicList<label> > collapseStrings; + + forAll(allPointInfo, pointI) + { + const label collapseIndex = allPointInfo[pointI].collapseIndex(); + + if (collapseIndex != -1 && collapseIndex != -2) + { + collapseStrings(collapseIndex).append(pointI); + } + } + + bool meshChanged = false; + + // Current faces (is also collapseStatus: f.size() < 3) + faceList newFaces(mesh.faces()); + + // Current cellCollapse status + boolList cellRemoved(mesh.nCells(), false); + + label nUnvisited = 0; + label nUncollapsed = 0; label nCollapsed = 0; - forAll(pointEdges, pointI) + forAll(allPointInfo, pI) { - const labelList& pEdges = pointEdges[pointI]; + const pointEdgeCollapse& pec = allPointInfo[pI]; - if (pEdges.size() == 2) + if (pec.collapseIndex() == -1) + { + nUnvisited++; + } + else if (pec.collapseIndex() == -2) { - const edge& leftE = edges[pEdges[0]]; - const edge& rightE = edges[pEdges[1]]; + nUncollapsed++; + } + else if (pec.collapseIndex() >= 0) + { + nCollapsed++; + } + } + + label nPoints = allPointInfo.size(); + + reduce(nPoints, sumOp<label>()); + reduce(nUnvisited, sumOp<label>()); + reduce(nUncollapsed, sumOp<label>()); + reduce(nCollapsed, sumOp<label>()); - // Get the two vertices on both sides of the point - label leftV = leftE.otherVertex(pointI); - label rightV = rightE.otherVertex(pointI); + Info<< incrIndent; + Info<< indent << "Number of points : " << nPoints << nl + << indent << "Not visited : " << nUnvisited << nl + << indent << "Not collapsed : " << nUncollapsed << nl + << indent << "Collapsed : " << nCollapsed << nl + << endl; + Info<< decrIndent; - // Collapse only if none of the points part of merge network - // or all of networks with different masters. - label midMaster = -1; - if (region[pointI] != -1) + do + { + // Update face collapse from edge collapses + forAll(newFaces, faceI) + { + filterFace(faceI, newFaces[faceI], allPointInfo, collapseStrings); + } + + // Check if faces to be collapsed cause cells to become collapsed. + label nCellCollapsed = 0; + + forAll(cells, cellI) + { + if (!cellRemoved[cellI]) { - midMaster = master[region[pointI]]; + const cell& cFaces = cells[cellI]; + + label nFaces = cFaces.size(); + + forAll(cFaces, i) + { + label faceI = cFaces[i]; + + if (newFaces[faceI].size() < 3) + { + --nFaces; + + if (nFaces < 4) + { + Info<< "Cell:" << cellI + << " uses faces:" << cFaces + << " of which too many are marked for removal:" + << endl + << " "; + forAll(cFaces, j) + { + if (newFaces[cFaces[j]].size() < 3) + { + Info<< ' '<< cFaces[j]; + } + } + Info<< endl; + + cellRemoved[cellI] = true; + + // Collapse all edges of cell to nothing + //collapseEdges(cellEdges[cellI]); + + nCellCollapsed++; + + break; + } + } + } } + } + + if (nCellCollapsed == 0) + { + break; + } + } while (true); + + + // Keep track of faces that have been done already. + boolList doneFace(mesh.nFaces(), false); + + { + // Mark points used. + boolList usedPoint(mesh.nPoints(), false); - label leftMaster = -2; - if (region[leftV] != -1) + forAll(cellRemoved, cellI) + { + if (cellRemoved[cellI]) { - leftMaster = master[region[leftV]]; + meshMod.removeCell(cellI, -1); } + } + + // Remove faces + forAll(newFaces, faceI) + { + const face& f = newFaces[faceI]; + + if (f.size() < 3) + { + meshMod.removeFace(faceI, -1); + meshChanged = true; - label rightMaster = -3; - if (region[rightV] != -1) + // Mark face as been done. + doneFace[faceI] = true; + } + else { - rightMaster = master[region[rightV]]; + // Kept face. Mark vertices + forAll(f, fp) + { + usedPoint[f[fp]] = true; + } } + } - if - ( - midMaster != leftMaster - && midMaster != rightMaster - && leftMaster != rightMaster - ) + // Remove unused vertices that have not been marked for removal already + forAll(usedPoint, pointI) + { + if (!usedPoint[pointI]) { - // Check if the two edge are in line - vector leftVec = points[pointI] - points[leftV]; - leftVec /= mag(leftVec) + VSMALL; + removedPoints[pointI] = true; + meshMod.removePoint(pointI, -1); + meshChanged = true; + } + } + } + + // Modify the point location of the remaining points + forAll(allPointInfo, pointI) + { + const label collapseIndex = allPointInfo[pointI].collapseIndex(); + const point& collapsePoint = allPointInfo[pointI].collapsePoint(); + + if + ( + removedPoints[pointI] == false + && collapseIndex != -1 + && collapseIndex != -2 + ) + { + meshMod.modifyPoint + ( + pointI, + collapsePoint, + pointZones.whichZone(pointI), + false + ); + } + } + + + const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh(); + const faceZoneMesh& faceZones = mesh.faceZones(); + + // Renumber faces that use points + forAll(allPointInfo, pointI) + { + if (removedPoints[pointI] == true) + { + const labelList& changedFaces = pointFaces[pointI]; - vector rightVec = points[rightV] - points[pointI]; - rightVec /= mag(rightVec) + VSMALL; + forAll(changedFaces, changedFaceI) + { + label faceI = changedFaces[changedFaceI]; - if ((leftVec & rightVec) > maxCos) + if (!doneFace[faceI]) { - // Collapse one (left) side of the edge. Make left vertex - // the master. - //if (collapser.unaffectedEdge(pEdges[0])) + doneFace[faceI] = true; + + // Get current zone info + label zoneID = faceZones.whichZone(faceI); + + bool zoneFlip = false; + + if (zoneID >= 0) { - collapser.collapseEdge(pEdges[0], leftV); - nCollapsed++; + const faceZone& fZone = faceZones[zoneID]; + + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + // Get current connectivity + label own = faceOwner[faceI]; + label nei = -1; + label patchID = -1; + + if (mesh.isInternalFace(faceI)) + { + nei = faceNeighbour[faceI]; + } + else + { + patchID = boundaryMesh.whichPatch(faceI); } + + meshMod.modifyFace + ( + newFaces[faceI], // face + faceI, // faceI to change + own, // owner + nei, // neighbour + false, // flipFaceFlux + patchID, // patch + zoneID, + zoneFlip + ); + meshChanged = true; } } } } - return nCollapsed; + // Print regions: +// printRegions(); + + return meshChanged; +} + + +// Get faceEdges in order of face points, i.e. faceEdges[0] is between +// f[0] and f[1] +labelList getSortedEdges +( + const edgeList& edges, + const labelList& f, + const labelList& edgeLabels +) +{ + labelList faceEdges(edgeLabels.size(), -1); + + // Find starting pos in f for every edgeLabels + forAll(edgeLabels, i) + { + label edgeI = edgeLabels[i]; + + const edge& e = edges[edgeI]; + + label fp = findIndex(f, e[0]); + label fp1 = f.fcIndex(fp); + + if (f[fp1] == e[1]) + { + // EdgeI between fp -> fp1 + faceEdges[fp] = edgeI; + } + else + { + // EdgeI between fp-1 -> fp + faceEdges[f.rcIndex(fp)] = edgeI; + } + } + + return faceEdges; } +// Merges edges which are in straight line. I.e. edge split by point. +//label mergeEdges +//( +// const polyMesh& mesh, +// const scalar maxCos, +// List<pointEdgeCollapse>& allPointInfo +//) +//{ +// const pointField& points = mesh.points(); +// const edgeList& edges = mesh.edges(); +// const labelListList& pointEdges = mesh.pointEdges(); +// const labelList& region = collapser.pointRegion(); +// const labelList& master = collapser.pointRegionMaster(); +// +// label nCollapsed = 0; +// +// forAll(pointEdges, pointI) +// { +// const labelList& pEdges = pointEdges[pointI]; +// +// if (pEdges.size() == 2) +// { +// const edge& leftE = edges[pEdges[0]]; +// const edge& rightE = edges[pEdges[1]]; +// +// // Get the two vertices on both sides of the point +// label leftV = leftE.otherVertex(pointI); +// label rightV = rightE.otherVertex(pointI); +// +// // Collapse only if none of the points part of merge network +// // or all of networks with different masters. +// label midMaster = -1; +// if (region[pointI] != -1) +// { +// midMaster = master[region[pointI]]; +// } +// +// label leftMaster = -2; +// if (region[leftV] != -1) +// { +// leftMaster = master[region[leftV]]; +// } +// +// label rightMaster = -3; +// if (region[rightV] != -1) +// { +// rightMaster = master[region[rightV]]; +// } +// +// if +// ( +// midMaster != leftMaster +// && midMaster != rightMaster +// && leftMaster != rightMaster +// ) +// { +// // Check if the two edge are in line +// vector leftVec = points[pointI] - points[leftV]; +// leftVec /= mag(leftVec) + VSMALL; +// +// vector rightVec = points[rightV] - points[pointI]; +// rightVec /= mag(rightVec) + VSMALL; +// +// if ((leftVec & rightVec) > maxCos) +// { +// // Collapse one (left) side of the edge. Make left vertex +// // the master. +// //if (collapser.unaffectedEdge(pEdges[0])) +// const edge& e = mesh.edges()[pEdges[0]]; +// +// if +// ( +// allPointInfo[e[0]].collapseIndex() < 0 +// && allPointInfo[e[1]].collapseIndex() < 0 +// ) +// { +// //pointEdgeCollapse pec(mesh.points()[leftV], leftV); +// +// allPointInfo[e[0]] = pec; +// allPointInfo[e[1]] = pec; +// +// //collapser.collapseEdge(pEdges[0], leftV); +// nCollapsed++; +// } +// } +// } +// } +// } +// +// return nCollapsed; +//} + + // Return master point edge needs to be collapsed to (or -1) -label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e) +label edgeMaster +( + const labelList& boundaryPoint, + const bool flipEdge, + const edge& e +) { label masterPoint = -1; + label e0 = e[0]; + label e1 = e[1]; + + if (flipEdge) + { + e0 = e[1]; + e1 = e[0]; + } + + // Check if one of the points is on a processor +// if +// ( +// boundaryPoint[e0] > 0 +// && boundaryPoint[e1] > 0 +// ) +// { +// if (boundaryPoint[e0] != boundaryPoint[e1]) +// { +// return -1; +// } +// } +// +// if (boundaryPoint[e0] > 0) +// { +// return e0; +// } +// else if (boundaryPoint[e1] > 0) +// { +// return e1; +// } + // Collapse edge to boundary point. - if (boundaryPoint.get(e[0])) + if (boundaryPoint[e0] == 0) { - if (boundaryPoint.get(e[1])) + if (boundaryPoint[e1] == 0) { // Both points on boundary. Choose one to collapse to. // Note: should look at feature edges/points! - masterPoint = e[0]; + masterPoint = e0; } else { - masterPoint = e[0]; + masterPoint = e0; } } else { - if (boundaryPoint.get(e[1])) + if (boundaryPoint[e1] == 0) { - masterPoint = e[1]; + masterPoint = e1; } else { // None on boundary. Choose arbitrary. // Note: should look at geometry? - masterPoint = e[0]; + masterPoint = e0; } } + return masterPoint; } @@ -212,124 +676,266 @@ label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e) label collapseSmallEdges ( const polyMesh& mesh, - const PackedBoolList& boundaryPoint, + const scalarList& freezeEdges, + const labelList& boundaryPoint, const scalar minLen, - edgeCollapser& collapser + List<pointEdgeCollapse>& allPointInfo ) { const pointField& points = mesh.points(); const edgeList& edges = mesh.edges(); - // Collapse all edges that are too small. Choose intelligently which - // point to collapse edge to. - - label nCollapsed = 0; + // Store collapse direction in collapseEdge + // -1 -> Do not collapse + // 0 -> Collapse to start point + // 1 -> Collapse to end point + labelList collapseEdge(edges.size(), -1); forAll(edges, edgeI) { const edge& e = edges[edgeI]; - if (e.mag(points) < minLen) + if (e.mag(points) < minLen*freezeEdges[edgeI]) { - label master = edgeMaster(boundaryPoint, e); - - if (master != -1) // && collapser.unaffectedEdge(edgeI)) - { - collapser.collapseEdge(edgeI, master); - nCollapsed++; - } + collapseEdge[edgeI] = 0; } } - return nCollapsed; -} - - -// Faces which have edges just larger than collapse length but faces which -// are very small. This one tries to collapse them if it can be done with -// edge collapse. For faces where a face gets replace by two edges use -// collapseFaces -label collapseHighAspectFaces -( - const polyMesh& mesh, - const PackedBoolList& boundaryPoint, - const scalar areaFac, - const scalar edgeRatio, - edgeCollapser& collapser -) -{ - const pointField& points = mesh.points(); - const edgeList& edges = mesh.edges(); - const faceList& faces = mesh.faces(); - const labelListList& faceEdges = mesh.faceEdges(); - scalarField magArea(mag(mesh.faceAreas())); + // Check whether edge point order is reversed from mesh to coupledPatch +// const globalMeshData& globalData = mesh.globalData(); +// const mapDistribute& map = globalData.globalEdgeSlavesMap(); +// const labelList& coupledMeshEdges = globalData.coupledPatchMeshEdges(); +// const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch(); +// const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation(); +// PackedBoolList meshToPatchSameOrientation(coupledMeshEdges.size(), true); +// +// forAll(coupledMeshEdges, eI) +// { +// const label meshEdgeIndex = coupledMeshEdges[eI]; +// +// if (collapseEdge[meshEdgeIndex] != -1) +// { +// const edge& meshEdge = edges[meshEdgeIndex]; +// const edge& coupledPatchEdge = coupledPatch.edges()[eI]; +// +// if +// ( +// meshEdge[0] == coupledPatch.meshPoints()[coupledPatchEdge[1]] +// && meshEdge[1] == coupledPatch.meshPoints()[coupledPatchEdge[0]] +// ) +// { +// meshToPatchSameOrientation[eI] = false; +// } +// } +// } +// +// +// labelList cppEdgeData(coupledMeshEdges.size(), -1); +// +// forAll(coupledMeshEdges, eI) +// { +// const label meshEdgeIndex = coupledMeshEdges[eI]; +// +// if (collapseEdge[meshEdgeIndex] != -1) +// { +// if (meshToPatchSameOrientation[eI] == cppOrientation[eI]) +// { +// cppEdgeData[eI] = 0; +// } +// else +// { +// cppEdgeData[eI] = 1; +// } +// } +// } +// +// +// // Synchronise cppEdgeData +// // Use minEqOp reduction, so that edge will only be collapsed on processor +// // boundary if both processors agree to collapse it +// globalData.syncData +// ( +// cppEdgeData, +// globalData.globalEdgeSlaves(), +// globalData.globalEdgeTransformedSlaves(), +// map, +// minEqOp<label>() +// ); +// +// +// forAll(coupledMeshEdges, eI) +// { +// const label meshEdgeIndex = coupledMeshEdges[eI]; +// +// if (collapseEdge[meshEdgeIndex] != -1) +// { +// if (meshToPatchSameOrientation[eI] == cppOrientation[eI]) +// { +// collapseEdge[meshEdgeIndex] = 0; +// } +// else +// { +// collapseEdge[meshEdgeIndex] = 1; +// } +// } +// } - label maxIndex = findMax(magArea); + label nCollapsed = 0; - scalar minArea = areaFac * magArea[maxIndex]; + DynamicList<label> initPoints(mesh.nPoints()); + DynamicList<pointEdgeCollapse> initPointInfo(mesh.nPoints()); + allPointInfo.resize(mesh.nPoints()); - Info<< "Max face area:" << magArea[maxIndex] << endl - << "Collapse area factor:" << areaFac << endl - << "Collapse area:" << minArea << endl; + globalIndex globalStrings(mesh.nPoints()); - label nCollapsed = 0; + List<pointEdgeCollapse> allEdgeInfo(mesh.nEdges()); + forAll(allEdgeInfo, edgeI) + { + allEdgeInfo[edgeI] = pointEdgeCollapse(vector::zero, -1); + } - forAll(faces, faceI) + forAll(edges, edgeI) { - if (magArea[faceI] < minArea) - { - const face& f = faces[faceI]; + const edge& e = edges[edgeI]; - // Get the edges in face point order - labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI])); + if (collapseEdge[edgeI] != -1) + { + const label master = + edgeMaster + ( + boundaryPoint, + collapseEdge[edgeI], + e + ); - SortableList<scalar> lengths(fEdges.size()); - forAll(fEdges, i) +// if (master != -1) { - lengths[i] = edges[fEdges[i]].mag(points); - } - lengths.sort(); - - - label edgeI = -1; + pointEdgeCollapse pec + ( + points[master], + globalStrings.toGlobal(master) + ); - if (f.size() == 4) - { - // Compare second largest to smallest - if (lengths[2] > edgeRatio*lengths[0]) - { - // Collapse smallest only. Triangle should be cleared - // next time around. - edgeI = fEdges[lengths.indices()[0]]; - } - } - else if (f.size() == 3) - { - // Compare second largest to smallest - if (lengths[1] > edgeRatio*lengths[0]) - { - edgeI = fEdges[lengths.indices()[0]]; - } - } + allEdgeInfo[edgeI] = pec; + initPointInfo.append(pec); + initPoints.append(e.start()); - if (edgeI != -1) - { - label master = edgeMaster(boundaryPoint, edges[edgeI]); + initPointInfo.append(pec); + initPoints.append(e.end()); - if (master != -1)// && collapser.unaffectedEdge(edgeI)) - { - collapser.collapseEdge(edgeI, master); - nCollapsed++; - } + nCollapsed++; } } } + PointEdgeWave<pointEdgeCollapse> collapsePropagator + ( + mesh, + initPoints, + initPointInfo, + allPointInfo, + allEdgeInfo, + mesh.globalData().nTotalPoints() // Maximum number of iterations + ); + return nCollapsed; } +// Faces which have edges just larger than collapse length but faces which +// are very small. This one tries to collapse them if it can be done with +// edge collapse. For faces where a face gets replace by two edges use +// collapseFaces +//label collapseHighAspectFaces +//( +// const polyMesh& mesh, +// const PackedBoolList& boundaryPoint, +// const Map<label>& processorPoints, +// const scalar areaFac, +// const scalar edgeRatio, +// edgeCollapser& collapser +//) +//{ +// const pointField& points = mesh.points(); +// const edgeList& edges = mesh.edges(); +// const faceList& faces = mesh.faces(); +// const labelListList& faceEdges = mesh.faceEdges(); +// +// scalarField magArea(mag(mesh.faceAreas())); +// +// label maxIndex = findMax(magArea); +// +// scalar minArea = areaFac * magArea[maxIndex]; +// +// Info<< "Max face area:" << magArea[maxIndex] << endl +// << "Collapse area factor:" << areaFac << endl +// << "Collapse area:" << minArea << endl; +// +// label nCollapsed = 0; +// +// forAll(faces, faceI) +// { +// if (magArea[faceI] < minArea) +// { +// const face& f = faces[faceI]; +// +// // Get the edges in face point order +// labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI])); +// +// SortableList<scalar> lengths(fEdges.size()); +// forAll(fEdges, i) +// { +// lengths[i] = edges[fEdges[i]].mag(points); +// } +// lengths.sort(); +// +// +// label edgeI = -1; +// +// if (f.size() == 4) +// { +// // Compare second largest to smallest +// if (lengths[2] > edgeRatio*lengths[0]) +// { +// // Collapse smallest only. Triangle should be cleared +// // next time around. +// edgeI = fEdges[lengths.indices()[0]]; +// } +// } +// else if (f.size() == 3) +// { +// // Compare second largest to smallest +// if (lengths[1] > edgeRatio*lengths[0]) +// { +// edgeI = fEdges[lengths.indices()[0]]; +// } +// } +// +// +// if (edgeI != -1) +// { +// label master = edgeMaster +// ( +// boundaryPoint, +// processorPoints, +// false, +// edges[edgeI] +// ); +// +// if (master != -1)// && collapser.unaffectedEdge(edgeI)) +// { +// collapser.collapseEdge(edgeI, master); +// nCollapsed++; +// } +// } +// } +// } +// +// return nCollapsed; +//} + + void set(const labelList& elems, const bool val, boolList& status) { forAll(elems, i) @@ -340,122 +946,254 @@ void set(const labelList& elems, const bool val, boolList& status) // Tries to simplify polygons to face of minSize (4=quad, 3=triangle) -label simplifyFaces -( - const polyMesh& mesh, - const PackedBoolList& boundaryPoint, - const label minSize, - const scalar lenGap, - edgeCollapser& collapser -) +//label simplifyFaces +//( +// const polyMesh& mesh, +// const PackedBoolList& boundaryPoint, +// const Map<label>& processorPoints, +// const label minSize, +// const scalar lenGap, +// edgeCollapser& collapser +//) +//{ +// const pointField& points = mesh.points(); +// const edgeList& edges = mesh.edges(); +// const faceList& faces = mesh.faces(); +// const cellList& cells = mesh.cells(); +// const labelListList& faceEdges = mesh.faceEdges(); +// const labelList& faceOwner = mesh.faceOwner(); +// const labelList& faceNeighbour = mesh.faceNeighbour(); +// const labelListList& pointCells = mesh.pointCells(); +// const labelListList& cellEdges = mesh.cellEdges(); +// +// label nCollapsed = 0; +// +// boolList protectedEdge(mesh.nEdges(), false); +// +// forAll(faces, faceI) +// { +// const face& f = faces[faceI]; +// +// if +// ( +// f.size() > minSize +// && cells[faceOwner[faceI]].size() >= 6 +// && ( +// mesh.isInternalFace(faceI) +// && cells[faceNeighbour[faceI]].size() >= 6 +// ) +// ) +// { +// // Get the edges in face point order +// labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI])); +// +// SortableList<scalar> lengths(fEdges.size()); +// forAll(fEdges, i) +// { +// lengths[i] = edges[fEdges[i]].mag(points); +// } +// lengths.sort(); +// +// +// // Now find a gap in length between consecutive elements greater +// // than lenGap. +// +// label gapPos = -1; +// +// for (label i = f.size()-1-minSize; i >= 0; --i) +// { +// if (lengths[i+1] > lenGap*lengths[i]) +// { +// gapPos = i; +// +// break; +// } +// } +// +// if (gapPos != -1) +// { +// //for (label i = gapPos; i >= 0; --i) +// label i = 0; // Hack: collapse smallest edge only. +// { +// label edgeI = fEdges[lengths.indices()[i]]; +// +// if (!protectedEdge[edgeI]) +// { +// const edge& e = edges[edgeI]; +// +// label master +// = edgeMaster(boundaryPoint, processorPoints, false, e); +// +// if (master != -1) +// { +// collapser.collapseEdge(edgeI, master); +// +// // Protect all other edges on all cells using edge +// // points. +// +// const labelList& pCells0 = pointCells[e[0]]; +// +// forAll(pCells0, i) +// { +// set(cellEdges[pCells0[i]], true, protectedEdge); +// } +// const labelList& pCells1 = pointCells[e[1]]; +// +// forAll(pCells1, i) +// { +// set(cellEdges[pCells1[i]], true, protectedEdge); +// } +// +// nCollapsed++; +// } +// } +// } +// } +// } +// } +// +// return nCollapsed; +//} + + +labelHashSet checkMeshQuality(const polyMesh& mesh) { - const pointField& points = mesh.points(); - const edgeList& edges = mesh.edges(); - const faceList& faces = mesh.faces(); - const cellList& cells = mesh.cells(); - const labelListList& faceEdges = mesh.faceEdges(); - const labelList& faceOwner = mesh.faceOwner(); - const labelList& faceNeighbour = mesh.faceNeighbour(); - const labelListList& pointCells = mesh.pointCells(); - const labelListList& cellEdges = mesh.cellEdges(); + //mesh.checkMesh(true); + labelHashSet freezePoints; - label nCollapsed = 0; + IOdictionary meshQualityDict + ( + IOobject + ( + "meshQualityControls", + mesh.time().system(), + mesh.time(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); - boolList protectedEdge(mesh.nEdges(), false); + labelHashSet badFaces(mesh.nFaces()/100); + DynamicList<label> checkFaces(mesh.nFaces()); - forAll(faces, faceI) - { - const face& f = faces[faceI]; + const vectorField& fAreas = mesh.faceAreas(); - if - ( - f.size() > minSize - && cells[faceOwner[faceI]].size() >= 6 - && ( - mesh.isInternalFace(faceI) - && cells[faceNeighbour[faceI]].size() >= 6 - ) - ) + scalar faceAreaLimit = SMALL; + + forAll(fAreas, fI) + { + if (mag(fAreas[fI]) > faceAreaLimit) { - // Get the edges in face point order - labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI])); + checkFaces.append(fI); + } + } - SortableList<scalar> lengths(fEdges.size()); - forAll(fEdges, i) - { - lengths[i] = edges[fEdges[i]].mag(points); - } - lengths.sort(); + motionSmoother::checkMesh + ( + false, + mesh, + meshQualityDict, + checkFaces, + badFaces + ); + label nBadFaces = badFaces.size(); + reduce(nBadFaces, sumOp<label>()); - // Now find a gap in length between consecutive elements greater - // than lenGap. + Info<< nl << "Number of bad faces : " << nBadFaces << endl; - label gapPos = -1; + forAllConstIter(labelHashSet, badFaces, iter) + { + const face& f = mesh.faces()[iter.key()]; - for (label i = f.size()-1-minSize; i >= 0; --i) - { - if (lengths[i+1] > lenGap*lengths[i]) - { - gapPos = i; + forAll(f, pI) + { + freezePoints.insert(f[pI]); + } + } - break; - } - } +// const edgeList& edges = mesh.edges(); +// +// label nFrozenEdges = 0; +// +// OFstream str("frozenEdges.obj"); +// +// label count = 0; +// forAll(edges, eI) +// { +// const edge& e = edges[eI]; +// +// if (freezePoints.found(e[0]) && freezePoints.found(e[1])) +// { +// freezeEdges[eI] = true; +// nFrozenEdges++; +// } +// } + + return freezePoints; +} - if (gapPos != -1) - { - //for (label i = gapPos; i >= 0; --i) - label i = 0; // Hack: collapse smallest edge only. - { - label edgeI = fEdges[lengths.indices()[i]]; - if (!protectedEdge[edgeI]) - { - const edge& e = edges[edgeI]; +// Mark boundary points +// boundaryPoint: +// + -1 : point not on boundary +// + 0 : point on a real boundary +// + >0 : point on a processor patch with that ID +labelList findBoundaryPoints(const polyMesh& mesh) +{ + const faceList& faces = mesh.faces(); + const polyBoundaryMesh& bMesh = mesh.boundaryMesh(); - label master = edgeMaster(boundaryPoint, e); - if (master != -1) - { - collapser.collapseEdge(edgeI, master); + labelList boundaryPoint(mesh.nPoints(), -1); - // Protect all other edges on all cells using edge - // points. + // Get all points on a boundary + label nIntFaces = mesh.nInternalFaces(); + for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++) + { + const face& f = faces[faceI]; - const labelList& pCells0 = pointCells[e[0]]; + forAll(f, fp) + { + boundaryPoint[f[fp]] = 0; + } + } - forAll(pCells0, i) - { - set(cellEdges[pCells0[i]], true, protectedEdge); - } - const labelList& pCells1 = pointCells[e[1]]; + // Get all processor boundary points and the processor patch label + // that they are on. + forAll(bMesh, patchI) + { + const polyPatch& patch = bMesh[patchI]; - forAll(pCells1, i) - { - set(cellEdges[pCells1[i]], true, protectedEdge); - } + if (isA<processorPolyPatch>(patch)) + { + const processorPolyPatch& pPatch = + refCast<const processorPolyPatch>(patch); - nCollapsed++; - } - } + forAll(pPatch, fI) + { + const face& f = pPatch[fI]; + + forAll(f, fp) + { + boundaryPoint[f[fp]] = patchI; } } } } - return nCollapsed; + return boundaryPoint; } // Main program: - int main(int argc, char *argv[]) { # include "addOverwriteOption.H" - argList::noParallel(); + argList::validArgs.append("edge length [m]"); argList::validArgs.append("merge angle (degrees)"); + argList::addOption("minLenFactor", "scalar", "edge length factor"); # include "setRootCase.H" # include "createTime.H" @@ -463,8 +1201,11 @@ int main(int argc, char *argv[]) # include "createPolyMesh.H" const word oldInstance = mesh.pointsInstance(); - const scalar minLen = args.argRead<scalar>(1); + scalar minLen = args.argRead<scalar>(1); const scalar angle = args.argRead<scalar>(2); + const scalar minLenFactor + = args.optionLookupOrDefault<scalar>("minLenFactor", 0.5); + const bool overwrite = args.optionFound("overwrite"); scalar maxCos = Foam::cos(degToRad(angle)); @@ -475,108 +1216,196 @@ int main(int argc, char *argv[]) << " degrees" << nl << endl; + Info<< "If an invalid mesh is generated then the edge length will be " << nl + << "multiplied by a factor of " << minLenFactor << " and collapsing " + << "will be reattempted" << nl << endl; bool meshChanged = false; - while (true) - { - const faceList& faces = mesh.faces(); + checkMeshQuality(mesh); - // Get all points on the boundary - PackedBoolList boundaryPoint(mesh.nPoints()); + autoPtr<fvMesh> fvMeshPtr; - label nIntFaces = mesh.nInternalFaces(); - for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++) - { - const face& f = faces[faceI]; + scalarList freezeEdges(mesh.nEdges(), 1.0); - forAll(f, fp) - { - boundaryPoint.set(f[fp], 1); - } - } + do + { + label nIterations = 0; + label nFrozenEdges = 0; - // Edge collapsing engine - edgeCollapser collapser(mesh); + fvMeshPtr.reset(new fvMesh(mesh)); + fvMesh& fvMeshRef = fvMeshPtr(); + // Contains new point label for original points + labelList pointMap(identity(mesh.nPoints())); - // Collapse all edges that are too small. - label nCollapsed = - collapseSmallEdges - ( - mesh, - boundaryPoint, - minLen, - collapser - ); - Info<< "Collapsing " << nCollapsed << " small edges" << endl; + scalarList tmpFreezeEdges = freezeEdges; + autoPtr<mapPolyMesh> morphMap; - // Remove midpoints on straight edges. - if (nCollapsed == 0) + while (true) { - nCollapsed = mergeEdges(mesh, maxCos, collapser); - Info<< "Collapsing " << nCollapsed << " in line edges" << endl; - } + Info<< "Iteration " << nIterations << incrIndent << endl; + labelList boundaryPoint = findBoundaryPoints(fvMeshRef); - // Remove small sliver faces that can be collapsed to single edge - if (nCollapsed == 0) - { - nCollapsed = - collapseHighAspectFaces + List<pointEdgeCollapse> allPointInfo; + + // Collapse all edges that are too small. + label nSmallCollapsed = + collapseSmallEdges ( - mesh, + fvMeshRef, + tmpFreezeEdges, boundaryPoint, - 1e-9, // factor of largest face area - 5, // factor between smallest and largest edge on - // face - collapser + minLen, + allPointInfo ); - Info<< "Collapsing " << nCollapsed + + reduce(nSmallCollapsed, sumOp<label>()); + + Info<< indent << "Collapsing " << nSmallCollapsed + << " small edges" << endl; + + + + + + label nMerged = 0; + + // Remove midpoints on straight edges. + if (nSmallCollapsed == 0) + { + //nMerged = mergeEdges(fvMeshRef, maxCos, allPointInfo); + } + + reduce(nMerged, sumOp<label>()); + + Info<< indent << "Collapsing " << nMerged << " in line edges" + << endl; + + + + + + + label nSliversCollapsed = 0; + + // Remove small sliver faces that can be collapsed to single edge +// if (nSmallCollapsed == 0 && nMerged == 0) +// { +// nSliversCollapsed = +// collapseHighAspectFaces +// ( +// mesh, +// boundaryPoint, +// processorPoints, +// 1E-9, +// 5, +// collapser +// ); +// } + + reduce(nSliversCollapsed, sumOp<label>()); + + Info<< indent << "Collapsing " << nSliversCollapsed << " small high aspect ratio faces" << endl; - } - // Simplify faces to quads wherever possible - //if (nCollapsed == 0) - //{ - // nCollapsed = - // simplifyFaces - // ( - // mesh, - // boundaryPoint, - // 4, // minimum size of face - // 0.2, // gap in edge lengths on face - // collapser - // ); - // Info<< "Collapsing " << nCollapsed << " polygonal faces" << endl; - //} + // Simplify faces to quads wherever possible + //if (nCollapsed == 0) + //{ + // nCollapsed = + // simplifyFaces + // ( + // mesh, + // boundaryPoint, + // 4, // minimum size of face + // 0.2, // gap in edge lengths on face + // collapser + // ); + // Info<< "Collapsing " << nCollapsed << " polygonal faces" + // << endl; + //} + + + label totalCollapsed = + nSmallCollapsed + + nMerged + + nSliversCollapsed; + + polyTopoChange meshMod(fvMeshRef); + + // Insert mesh refinement into polyTopoChange. + setRefinement(fvMeshRef, meshMod, allPointInfo); + + // Do all changes + Info<< indent << "Applying changes to the mesh" << nl + << decrIndent << endl; + + morphMap = meshMod.changeMesh(fvMeshRef, false); + +// // Contains new point label for old points +// const labelList& reversePointMap = morphMap().reversePointMap(); +// +// forAll(pointMap, pI) +// { +// const label originalPoint = pI; +// const label currentPoint = pointMap[pI]; +// +// if (currentPoint < reversePointMap.size()) +// { +// const label newPoint = reversePointMap[currentPoint]; +// +// if (newPoint != -1) +// { +// pointMap[originalPoint] = newPoint; +// } +// } +// } + + if (totalCollapsed == 0) + { + labelHashSet freezePoints = checkMeshQuality(fvMeshRef); - if (nCollapsed == 0) - { - break; - } + label nFreezePoints = freezePoints.size(); + reduce(nFreezePoints, sumOp<label>()); - polyTopoChange meshMod(mesh); + nFrozenEdges = nFreezePoints; - // Insert mesh refinement into polyTopoChange. - collapser.setRefinement(meshMod); + Info<< "Number of frozen points : " << nFreezePoints + << endl; - // Do all changes - Info<< "Morphing ..." << endl; + break; + } - autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false); + if (morphMap().hasMotionPoints()) + { + fvMeshRef.movePoints(morphMap().preMotionPoints()); + } - collapser.updateMesh(morphMap()); + meshChanged = true; - if (morphMap().hasMotionPoints()) + nIterations++; + } + + if (nFrozenEdges > 0) { - mesh.movePoints(morphMap().preMotionPoints()); + minLen *= minLenFactor; } - meshChanged = true; - } + reduce(nFrozenEdges, sumOp<label>()); + + Info<< "Number of frozen edges : " << nFrozenEdges << nl + << endl; + + if (nFrozenEdges == 0) + { + break; + } + + } while (true); + if (meshChanged) { @@ -587,14 +1416,21 @@ int main(int argc, char *argv[]) } else { - mesh.setInstance(oldInstance); + fvMeshPtr().setInstance(oldInstance); } - Info<< "Writing collapsed mesh to time " << runTime.timeName() << endl; + Info<< nl << "Writing collapsed mesh to time " + << runTime.timeName() << nl << endl; - mesh.write(); + fvMeshPtr().write(); } + Info<< "Final minimum length : " << minLen << " m" << nl << endl; + + Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + Info<< "End\n" << endl; return 0; diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C new file mode 100644 index 0000000000000000000000000000000000000000..6830e1c8f56ed18c006241ffb16388012da4b042 --- /dev/null +++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "pointEdgeCollapse.H" + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<< +( + Foam::Ostream& os, + const Foam::pointEdgeCollapse& wDist +) +{ + return os + << wDist.collapsePoint_ << wDist.collapseIndex_; +} + +Foam::Istream& Foam::operator>> +( + Foam::Istream& is, + Foam::pointEdgeCollapse& wDist +) +{ + return is + >> wDist.collapsePoint_ >> wDist.collapseIndex_; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H new file mode 100644 index 0000000000000000000000000000000000000000..cd8383b8b33b6ed6ddb0b03a172aa632205f3cee --- /dev/null +++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H @@ -0,0 +1,227 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::pointEdgeCollapse + +Description + Determines length of string of edges walked to point. + +SourceFiles + pointEdgeCollapseI.H + pointEdgeCollapse.C + +\*---------------------------------------------------------------------------*/ + +#ifndef pointEdgeCollapse_H +#define pointEdgeCollapse_H + +#include "point.H" +#include "tensor.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class polyPatch; +class polyMesh; + +/*---------------------------------------------------------------------------*\ + Class pointEdgeCollapse Declaration +\*---------------------------------------------------------------------------*/ + +class pointEdgeCollapse +{ + // Private data + + //- Collapse location + point collapsePoint_; + + //- Collapse string index + label collapseIndex_; + + + // Private Member Functions + + //- Evaluate distance to point. + template<class TrackingData> + inline bool update + ( + const pointEdgeCollapse& w2, + const scalar tol, + TrackingData& td + ); + + + //- Check for same coordinate + inline bool samePoint(const point& pt) const; + +public: + + // Constructors + + //- Construct null + inline pointEdgeCollapse(); + + //- Construct from components + inline pointEdgeCollapse + ( + const point& collapsePoint, + const label collapseIndex + ); + + + // Member Functions + + // Access + + inline const point& collapsePoint() const; + + inline label collapseIndex() const; + + + // Needed by meshWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value. + template<class TrackingData> + inline bool valid(TrackingData& td) const; + + //- Convert origin to relative vector to leaving point + // (= point coordinate) + template<class TrackingData> + inline void leaveDomain + ( + const polyPatch& patch, + const label patchPointI, + const point& pos, + TrackingData& td + ); + + //- Convert relative origin to absolute by adding entering point + template<class TrackingData> + inline void enterDomain + ( + const polyPatch& patch, + const label patchPointI, + const point& pos, + TrackingData& td + ); + + //- Apply rotation matrix to origin + template<class TrackingData> + inline void transform + ( + const tensor& rotTensor, + TrackingData& td + ); + + //- Influence of edge on point + template<class TrackingData> + inline bool updatePoint + ( + const polyMesh& mesh, + const label pointI, + const label edgeI, + const pointEdgeCollapse& edgeInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of different value on same point. + // Merge new and old info. + template<class TrackingData> + inline bool updatePoint + ( + const polyMesh& mesh, + const label pointI, + const pointEdgeCollapse& newPointInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of different value on same point. + // No information about current position whatsoever. + template<class TrackingData> + inline bool updatePoint + ( + const pointEdgeCollapse& newPointInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of point on edge. + template<class TrackingData> + inline bool updateEdge + ( + const polyMesh& mesh, + const label edgeI, + const label pointI, + const pointEdgeCollapse& pointInfo, + const scalar tol, + TrackingData& td + ); + + //- Same (like operator==) + template<class TrackingData> + inline bool equal(const pointEdgeCollapse&, TrackingData&) + const; + + + // Member Operators + + //Note: Used to determine whether to call update. + inline bool operator==(const pointEdgeCollapse&) const; + inline bool operator!=(const pointEdgeCollapse&) const; + + + // IOstream Operators + + friend Ostream& operator<<(Ostream&, const pointEdgeCollapse&); + friend Istream& operator>>(Istream&, pointEdgeCollapse&); +}; + + +//- Data associated with pointEdgeCollapse type are contiguous +template<> +inline bool contiguous<pointEdgeCollapse>() +{ + return true; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "pointEdgeCollapseI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H new file mode 100644 index 0000000000000000000000000000000000000000..a064bcf781dad9cbdbc24f7fc617269f6b355f7f --- /dev/null +++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H @@ -0,0 +1,313 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "polyMesh.H" +#include "transform.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Update this with w2. +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::update +( + const pointEdgeCollapse& w2, + const scalar tol, + TrackingData& td +) +{ + if (!w2.valid(td)) + { + FatalErrorIn("pointEdgeCollapse::update(..)") + << "problem." << abort(FatalError); + } + + if (w2.collapseIndex_ == -1) + { + // Not marked for collapse; only happens on edges. + return false; + } + + if (!valid(td)) + { + operator=(w2); + return true; + } + else + { + // Same coordinate. Same string? + if (w2.collapseIndex_ < collapseIndex_) + { + // Take over string index from w2 (and also coordinate but this + // was same) + operator=(w2); + return true; + } + else if (w2.collapseIndex_ == collapseIndex_) + { + bool identicalPoint = samePoint(w2.collapsePoint_); + bool nearer = magSqr(w2.collapsePoint_) < magSqr(collapsePoint_); + if (nearer) + { + operator=(w2); + } + if (identicalPoint) + { + return false; + } + else + { + return nearer; + } + } + else + { + return false; + } + +// if (samePoint(w2.collapsePoint_)) +// { +// // Same coordinate. Same string? +// if (w2.collapseIndex_ < collapseIndex_) +// { +// // Take over string index from w2 (and also coordinate but +// // this was same) +// operator=(w2); +// return true; +// } +// else +// { +// return false; +// } +// } +// else +// { +// // Find nearest coordinate +// if (magSqr(w2.collapsePoint_) < magSqr(collapsePoint_)) +// { +// operator=(w2); +// return true; +// } +// else +// { +// return false; +// } +// } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Null constructor +inline Foam::pointEdgeCollapse::pointEdgeCollapse() +: + collapsePoint_(GREAT, GREAT, GREAT), + collapseIndex_(-2) +{} + + +// Construct from origin, distance +inline Foam::pointEdgeCollapse::pointEdgeCollapse +( + const point& collapsePoint, + const label collapseIndex +) +: + collapsePoint_(collapsePoint), + collapseIndex_(collapseIndex) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::point& Foam::pointEdgeCollapse::collapsePoint() const +{ + return collapsePoint_; +} + + +inline Foam::label Foam::pointEdgeCollapse::collapseIndex() const +{ + return collapseIndex_; +} + + +inline bool Foam::pointEdgeCollapse::samePoint(const point& pt) const +{ + bool isLegal1 = (cmptMin(collapsePoint_) < 0.5*GREAT); + bool isLegal2 = (cmptMin(pt) < 0.5*GREAT); + + if (isLegal1 && isLegal2) + { + return mag(collapsePoint_ - pt) < 1e-9;//SMALL; + } + else + { + return isLegal1 == isLegal2; + } +} + + +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::valid(TrackingData& td) const +{ + return collapseIndex_ != -2; +} + + +template<class TrackingData> +inline void Foam::pointEdgeCollapse::leaveDomain +( + const polyPatch& patch, + const label patchPointI, + const point& coord, + TrackingData& td +) +{ + collapsePoint_ -= coord; +} + + +template<class TrackingData> +inline void Foam::pointEdgeCollapse::transform +( + const tensor& rotTensor, + TrackingData& td +) +{ + collapsePoint_ = Foam::transform(rotTensor, collapsePoint_); +} + + +// Update absolute geometric quantities. Note that distance (dist_) +// is not affected by leaving/entering domain. +template<class TrackingData> +inline void Foam::pointEdgeCollapse::enterDomain +( + const polyPatch& patch, + const label patchPointI, + const point& coord, + TrackingData& td +) +{ + // back to absolute form + collapsePoint_ += coord; +} + + +// Update this with information from connected edge +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::updatePoint +( + const polyMesh& mesh, + const label pointI, + const label edgeI, + const pointEdgeCollapse& edgeInfo, + const scalar tol, + TrackingData& td +) +{ + return update(edgeInfo, tol, td); +} + + +// Update this with new information on same point +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::updatePoint +( + const polyMesh& mesh, + const label pointI, + const pointEdgeCollapse& newPointInfo, + const scalar tol, + TrackingData& td +) +{ + return update(newPointInfo, tol, td); +} + + +// Update this with new information on same point. No extra information. +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::updatePoint +( + const pointEdgeCollapse& newPointInfo, + const scalar tol, + TrackingData& td +) +{ + return update(newPointInfo, tol, td); +} + + +// Update this with information from connected point +template<class TrackingData> +inline bool Foam::pointEdgeCollapse::updateEdge +( + const polyMesh& mesh, + const label edgeI, + const label pointI, + const pointEdgeCollapse& pointInfo, + const scalar tol, + TrackingData& td +) +{ + return update(pointInfo, tol, td); +} + + +template <class TrackingData> +inline bool Foam::pointEdgeCollapse::equal +( + const pointEdgeCollapse& rhs, + TrackingData& td +) const +{ + return operator==(rhs); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline bool Foam::pointEdgeCollapse::operator== +( + const Foam::pointEdgeCollapse& rhs +) const +{ + return + collapseIndex_ == rhs.collapseIndex_ + && samePoint(rhs.collapsePoint_); +} + + +inline bool Foam::pointEdgeCollapse::operator!= +( + const Foam::pointEdgeCollapse& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 0e760387ebad0b9b6410a15f21d8ef8c8bbb1592..074d5b8ee3405f62c8dd80b1d068a22a54d5a0c6 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -160,7 +160,7 @@ castellatedMeshControls // freedom in mesh motion // boundary : create free-standing boundary faces (baffles // but without the shared points) - //faceType internal; + //faceType baffle; } } diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C index d1a9b8f0288287e5f55d5e1fe6bbfed9d5d8c2ac..3393edb00d9ce90587892be0dc674abb716d8428 100644 --- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C +++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C @@ -71,6 +71,7 @@ Description #include "triSurface.H" #include "argList.H" #include "Time.H" +#include "featureEdgeMesh.H" #include "extendedFeatureEdgeMesh.H" #include "triSurfaceSearch.H" #include "OFstream.H" @@ -452,6 +453,30 @@ int main(int argc, char *argv[]) feMesh.writeObj(sFeatFileName); + { + // Write a featureEdgeMesh for backwards compatibility + featureEdgeMesh bfeMesh + ( + IOobject + ( + sFeatFileName + ".eMesh", // name + runTime.constant(), // instance + "triSurface", + runTime, // registry + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + feMesh.points(), + feMesh.edges() + ); + + Info<< nl << "Writing featureEdgeMesh to " + << bfeMesh.objectPath() << endl; + + bfeMesh.regIOobject::write(); + } + Info << "End\n" << endl; return 0; diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index 5aae0ede0ebf81adb979a655c6211425b9c95875..5dd7671b948c6286157ef84b8f138b1fde1c6c2a 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C @@ -913,8 +913,8 @@ int main(int argc, char *argv[]) // Info<< "span " << span << endl; - pointField start = searchSurf.faceCentres() - span*normals; - pointField end = searchSurf.faceCentres() + span*normals; + pointField start(searchSurf.faceCentres() - span*normals); + pointField end(searchSurf.faceCentres() + span*normals); const pointField& faceCentres = searchSurf.faceCentres(); List<List<pointIndexHit> > allHitInfo; diff --git a/applications/utilities/surface/surfaceSplitByTopology/surfaceSplitByTopology.C b/applications/utilities/surface/surfaceSplitByTopology/surfaceSplitByTopology.C index 118e608033962c1cfede4c38270ee2aa2d49891a..93a7caa3655a0df6d00ce01489690874e0732de5 100644 --- a/applications/utilities/surface/surfaceSplitByTopology/surfaceSplitByTopology.C +++ b/applications/utilities/surface/surfaceSplitByTopology/surfaceSplitByTopology.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -39,7 +39,7 @@ using namespace Foam; int main(int argc, char *argv[]) { argList::noParallel(); - argList::validArgs.clear(); + argList::validOptions.clear(); argList::validArgs.append("input surface file"); argList::validArgs.append("output surface file"); argList args(argc, argv); diff --git a/bin/foamClearPolyMesh b/bin/foamClearPolyMesh index 7e5ee1fe16b8ac66c7435214b4ed931e866c5334..a0c955a046bd172f6129e173f91e056ca781fd7e 100755 --- a/bin/foamClearPolyMesh +++ b/bin/foamClearPolyMesh @@ -3,7 +3,7 @@ # ========= | # \\ / F ield | OpenFOAM: The Open Source CFD Toolbox # \\ / O peration | -# \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation +# \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation # \\/ M anipulation | #------------------------------------------------------------------------------- # License @@ -133,6 +133,7 @@ for i in \ sets \ cellLevel \ pointLevel \ + level0Edge \ refinementHistory \ surfaceIndex \ ; diff --git a/etc/config/settings.csh b/etc/config/settings.csh index 94e5618bff1c69a1c2481be98138fa08ad2dd907..8c9ac17c7a822d3aa1f06d765fc61fa5be71650d 100644 --- a/etc/config/settings.csh +++ b/etc/config/settings.csh @@ -213,6 +213,13 @@ case ThirdParty: set mpfr_version=mpfr-3.1.0 set mpc_version=mpc-0.9 breaksw + case Gcc47: + case Gcc47++0x: + set gcc_version=gcc-4.7.0 + set gmp_version=gmp-5.0.4 + set mpfr_version=mpfr-3.1.0 + set mpc_version=mpc-0.9 + breaksw case Gcc45: case Gcc45++0x: set gcc_version=gcc-4.5.2 diff --git a/etc/config/settings.sh b/etc/config/settings.sh index dc39c480e1f09433e083ef98e0987b0825536915..0b12e4e453881a79484bfd374b58bc64709e51b5 100644 --- a/etc/config/settings.sh +++ b/etc/config/settings.sh @@ -233,7 +233,13 @@ case "${foamCompiler}" in OpenFOAM | ThirdParty) case "$WM_COMPILER" in Gcc | Gcc++0x | Gcc46 | Gcc46++0x) - gcc_version=gcc-4.6.2 + gcc_version=gcc-4.6.1 + gmp_version=gmp-5.0.4 + mpfr_version=mpfr-3.1.0 + mpc_version=mpc-0.9 + ;; + Gcc47 | Gcc47++0x) + gcc_version=gcc-4.7.0 gmp_version=gmp-5.0.4 mpfr_version=mpfr-3.1.0 mpc_version=mpc-0.9 diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C index 31a734740d0fb814dfc0971dfd6bc99f865637e4..030e7639f4b69d3e3b6ce60104efb46cec6fb4ce 100644 --- a/src/OSspecific/POSIX/fileMonitor.C +++ b/src/OSspecific/POSIX/fileMonitor.C @@ -32,6 +32,7 @@ License #include "regIOobject.H" // for fileModificationSkew symbol #ifdef FOAM_USE_INOTIFY +# include <unistd.h> # include <sys/inotify.h> # include <sys/ioctl.h> # include <errno.h> diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 3c054e2bfda92ea360122f0b79f3e29dca3188c1..614a7cdae59f39b501b7b296ec9d81e3f6cac746 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -240,7 +240,6 @@ lduMatrix = matrices/lduMatrix $(lduMatrix)/lduMatrix/lduMatrix.C $(lduMatrix)/lduMatrix/lduMatrixOperations.C $(lduMatrix)/lduMatrix/lduMatrixATmul.C -$(lduMatrix)/lduMatrix/lduMatrixTests.C $(lduMatrix)/lduMatrix/lduMatrixUpdateMatrixInterfaces.C $(lduMatrix)/lduMatrix/lduMatrixSolver.C $(lduMatrix)/lduMatrix/lduMatrixSmoother.C @@ -316,6 +315,7 @@ meshes/lduMesh/lduMesh.C LduMatrix = matrices/LduMatrix $(LduMatrix)/LduMatrix/lduMatrices.C +$(LduMatrix)/LduMatrix/solverPerformance.C $(LduMatrix)/LduMatrix/LduInterfaceField/LduInterfaceFields.C $(LduMatrix)/Smoothers/lduSmoothers.C $(LduMatrix)/Preconditioners/lduPreconditioners.C diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C index 30fc2fc0f69fe8fb9018dfc057788f6ffd9013ad..779f134de372034089edf1de3ed1ec0fc81768fc 100644 --- a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C @@ -57,8 +57,12 @@ autoPtr<interpolationWeights> interpolationWeights::New const scalarField& samples ) { - Info<< nl << "Selecting interpolationWeights " - << type << endl; + if (debug) + { + InfoIn("interpolationWeights::New") + << "Selecting interpolationWeights " + << type << endl; + } wordConstructorTable::iterator cstrIter = wordConstructorTablePtr_->find(type); @@ -87,33 +91,6 @@ interpolationWeights::~interpolationWeights() {} -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -//objectRegistry& interpolationWeights::registry -//( -// const objectRegistry& obr, -// const word& name -//) -//{ -// if (!obr.foundObject<objectRegistry>(name)) -// { -// objectRegistry* fieldsCachePtr = new objectRegistry -// ( -// IOobject -// ( -// name, -// obr.time().constant(), -// obr, -// IOobject::NO_READ, -// IOobject::NO_WRITE -// ) -// ); -// fieldsCachePtr->store(); -// } -// return const_cast<objectRegistry&>(obr.subRegistry(name)); -//} - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C index 97b2b2f40e9e34dca76899963571b3858f63eca6..fb131308b37e5773bae74f85525e529d5f11f574 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C @@ -382,7 +382,6 @@ Foam::Ostream& Foam::operator<< #include "LduMatrixOperations.C" #include "LduMatrixATmul.C" #include "LduMatrixUpdateMatrixInterfaces.C" -#include "LduMatrixTests.C" #include "LduMatrixPreconditioner.C" #include "LduMatrixSmoother.C" #include "LduMatrixSolver.C" diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H index e0ce6d28ba45f25bafe12b92bfaa1ae89ef11c63..c2414b0f8bc24a44a14a4421f8368a2a437d939f 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H @@ -54,6 +54,7 @@ SourceFiles #include "Field.H" #include "FieldField.H" #include "LduInterfaceFieldPtrsList.H" +#include "SolverPerformance.H" #include "typeInfo.H" #include "autoPtr.H" #include "runTimeSelectionTables.H" @@ -75,27 +76,6 @@ Ostream& operator<< const LduMatrix<Type, DType, LUType>& ); -template<class Type, class DType, class LUType> -typename LduMatrix<Type, DType, LUType>::solverPerformance max -( - const typename LduMatrix<Type, DType, LUType>::solverPerformance&, - const typename LduMatrix<Type, DType, LUType>::solverPerformance& -); - -template<class Type, class DType, class LUType> -Istream& operator>> -( - Istream&, - typename LduMatrix<Type, DType, LUType>::solverPerformance& -); - -template<class Type, class DType, class LUType> -Ostream& operator<< -( - Ostream&, - const typename LduMatrix<Type, DType, LUType>::solverPerformance& -); - /*---------------------------------------------------------------------------*\ Class LduMatrix Declaration @@ -127,168 +107,7 @@ class LduMatrix public: - //- Class returned by the solver - // containing performance statistics - class solverPerformance - { - word solverName_; - word fieldName_; - Type initialResidual_; - Type finalResidual_; - label noIterations_; - bool converged_; - FixedList<bool, pTraits<Type>::nComponents> singular_; - - - public: - - // Constructors - - solverPerformance() - : - initialResidual_(pTraits<Type>::zero), - finalResidual_(pTraits<Type>::zero), - noIterations_(0), - converged_(false), - singular_(false) - {} - - - solverPerformance - ( - const word& solverName, - const word& fieldName, - const Type& iRes = pTraits<Type>::zero, - const Type& fRes = pTraits<Type>::zero, - const label nIter = 0, - const bool converged = false, - const bool singular = false - ) - : - solverName_(solverName), - fieldName_(fieldName), - initialResidual_(iRes), - finalResidual_(fRes), - noIterations_(nIter), - converged_(converged), - singular_(singular) - {} - - - // Member functions - - //- Return solver name - const word& solverName() const - { - return solverName_; - } - - //- Return solver name - word& solverName() - { - return solverName_; - } - - - //- Return field name - const word& fieldName() const - { - return fieldName_; - } - - - //- Return initial residual - const Type& initialResidual() const - { - return initialResidual_; - } - - //- Return initial residual - Type& initialResidual() - { - return initialResidual_; - } - - - //- Return final residual - const Type& finalResidual() const - { - return finalResidual_; - } - - //- Return final residual - Type& finalResidual() - { - return finalResidual_; - } - - - //- Return number of iterations - label nIterations() const - { - return noIterations_; - } - - //- Return number of iterations - label& nIterations() - { - return noIterations_; - } - - - //- Has the solver converged? - bool converged() const - { - return converged_; - } - - //- Is the matrix singular? - bool singular() const; - - //- Check, store and return convergence - bool converged - ( - const Type& tolerance, - const Type& relTolerance - ); - - //- Singularity test - bool checkSingularity(const Type& residual); - - //- Print summary of solver performance to the given stream - void print(Ostream& os) const; - - - // Member Operators - - bool operator!=(const solverPerformance&) const; - - - // Friend functions - - //- Return the element-wise maximum of two solverPerformances - friend solverPerformance max <Type, DType, LUType> - ( - const solverPerformance&, - const solverPerformance& - ); - - - // Ostream Operator - - friend Istream& operator>> <Type, DType, LUType> - ( - Istream&, - solverPerformance& - ); - - friend Ostream& operator<< <Type, DType, LUType> - ( - Ostream&, - const solverPerformance& - ); - }; - + friend class SolverPerformance<Type>; //- Abstract base-class for LduMatrix solvers class solver @@ -417,7 +236,7 @@ public: //- Read and reset the solver parameters from the given dictionary virtual void read(const dictionary& solverDict); - virtual solverPerformance solve + virtual SolverPerformance<Type> solve ( Field<Type>& psi ) const = 0; @@ -641,15 +460,6 @@ public: // Declare name of the class and its debug switch ClassName("LduMatrix"); - //- Large Type for the use in solvers - static const scalar great_; - - //- Small Type for the use in solvers - static const scalar small_; - - //- Very small Type for the use in solvers - static const scalar vsmall_; - // Constructors @@ -863,16 +673,6 @@ typedef Foam::LduMatrix<Type, DType, LUType> \ defineNamedTemplateTypeNameAndDebug(ldu##Type##DType##LUType##Matrix, 0); \ \ \ -template<> \ -const scalar ldu##Type##DType##LUType##Matrix::great_(1e15); \ - \ -template<> \ -const scalar ldu##Type##DType##LUType##Matrix::small_(1e-15); \ - \ -template<> \ -const scalar ldu##Type##DType##LUType##Matrix::vsmall_(VSMALL); \ - \ - \ typedef LduMatrix<Type, DType, LUType>::smoother \ ldu##Type##DType##LUType##Smoother; \ \ diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C index c31af469f14b314751293dcecb81d1955d33c510..50831a0ec94afeac0f52aeca75926b7e31ad5a90 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C @@ -179,7 +179,7 @@ Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor return stabilise ( gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)), - matrix_.small_ + SolverPerformance<Type>::small_ ); // At convergence this simpler method is equivalent to the above diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C similarity index 70% rename from src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C rename to src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C index f5a675f0440eae2d0e70fd6b456a8ac6fe2c9d2f..5f307718eb200ef4f420a2bdb7222b337b123682 100644 --- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C @@ -23,27 +23,28 @@ License \*---------------------------------------------------------------------------*/ -#include "LduMatrix.H" +#include "SolverPerformance.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -template<class Type, class DType, class LUType> -bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::checkSingularity +template<class Type> +bool Foam::SolverPerformance<Type>::checkSingularity ( const Type& wApA ) { for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++) { - singular_[cmpt] = component(wApA, cmpt) < vsmall_; + singular_[cmpt] = + component(wApA, cmpt) < vsmall_; } return singular(); } -template<class Type, class DType, class LUType> -bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const +template<class Type> +bool Foam::SolverPerformance<Type>::singular() const { for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++) { @@ -54,8 +55,8 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const } -template<class Type, class DType, class LUType> -bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged +template<class Type> +bool Foam::SolverPerformance<Type>::checkConvergence ( const Type& Tolerance, const Type& RelTolerance @@ -73,7 +74,8 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged ( finalResidual_ < Tolerance || ( - RelTolerance > small_*pTraits<Type>::one + RelTolerance + > small_*pTraits<Type>::one && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_) ) ) @@ -89,16 +91,24 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged } -template<class Type, class DType, class LUType> -void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print +template<class Type> +void Foam::SolverPerformance<Type>::print ( Ostream& os ) const { for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++) { - os << solverName_ << ": Solving for " - << word(fieldName_ + pTraits<Type>::componentNames[cmpt]); + if (pTraits<Type>::nComponents == 1) + { + os << solverName_ << ": Solving for " << fieldName_; + + } + else + { + os << solverName_ << ": Solving for " + << word(fieldName_ + pTraits<Type>::componentNames[cmpt]); + } if (singular_[cmpt]) { @@ -115,10 +125,10 @@ void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print } -template<class Type, class DType, class LUType> -bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::operator!= +template<class Type> +bool Foam::SolverPerformance<Type>::operator!= ( - const LduMatrix<Type, DType, LUType>::solverPerformance& sp + const SolverPerformance<Type>& sp ) const { return @@ -134,14 +144,14 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::operator!= } -template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance Foam::max +template<class Type> +typename Foam::SolverPerformance<Type> Foam::max ( - const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp1, - const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp2 + const typename Foam::SolverPerformance<Type>& sp1, + const typename Foam::SolverPerformance<Type>& sp2 ) { - return lduMatrix::solverPerformance + return SolverPerformance<Type> ( sp1.solverName(), sp1.fieldName_, @@ -154,14 +164,14 @@ typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance Foam::max } -template<class Type, class DType, class LUType> +template<class Type> Foam::Istream& Foam::operator>> ( Istream& is, - typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp + typename Foam::SolverPerformance<Type>& sp ) { - is.readBeginList("LduMatrix<Type, DType, LUType>::solverPerformance"); + is.readBeginList("SolverPerformance<Type>"); is >> sp.solverName_ >> sp.fieldName_ >> sp.initialResidual_ @@ -169,17 +179,17 @@ Foam::Istream& Foam::operator>> >> sp.noIterations_ >> sp.converged_ >> sp.singular_; - is.readEndList("LduMatrix<Type, DType, LUType>::solverPerformance"); + is.readEndList("SolverPerformance<Type>"); return is; } -template<class Type, class DType, class LUType> +template<class Type> Foam::Ostream& Foam::operator<< ( Ostream& os, - const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp + const typename Foam::SolverPerformance<Type>& sp ) { os << token::BEGIN_LIST diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H new file mode 100644 index 0000000000000000000000000000000000000000..d39ecb10a80edf07b7afc15004eaad701cd440da --- /dev/null +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H @@ -0,0 +1,292 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::SolverPerformance + +Description + SolverPerformance is a general matrix class in which the coefficients are + stored as three arrays, one for the upper triangle, one for the + lower triangle and a third for the diagonal. + +SourceFiles + SolverPerformance.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SolverPerformance_H +#define SolverPerformance_H + +#include "word.H" +#include "FixedList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of friend functions and operators + +template<class Type> +class SolverPerformance; + +template<class Type> +SolverPerformance<Type> max +( + const SolverPerformance<Type>&, + const SolverPerformance<Type>& +); + +template<class Type> +Istream& operator>> +( + Istream&, + SolverPerformance<Type>& +); + +template<class Type> +Ostream& operator<< +( + Ostream&, + const SolverPerformance<Type>& +); + + + +/*---------------------------------------------------------------------------*\ + Class SolverPerformance Declaration +\*---------------------------------------------------------------------------*/ + +//- Class returned by the solver +// containing performance statistics +template<class Type> +class SolverPerformance +{ + word solverName_; + word fieldName_; + Type initialResidual_; + Type finalResidual_; + label noIterations_; + bool converged_; + FixedList<bool, pTraits<Type>::nComponents> singular_; + + +public: + + // Static data + + // Declare name of the class and its debug switch + ClassName("SolverPerformance"); + + //- Large Type for the use in solvers + static const scalar great_; + + //- Small Type for the use in solvers + static const scalar small_; + + //- Very small Type for the use in solvers + static const scalar vsmall_; + + + // Constructors + + SolverPerformance() + : + initialResidual_(pTraits<Type>::zero), + finalResidual_(pTraits<Type>::zero), + noIterations_(0), + converged_(false), + singular_(false) + {} + + + SolverPerformance + ( + const word& solverName, + const word& fieldName, + const Type& iRes = pTraits<Type>::zero, + const Type& fRes = pTraits<Type>::zero, + const label nIter = 0, + const bool converged = false, + const bool singular = false + ) + : + solverName_(solverName), + fieldName_(fieldName), + initialResidual_(iRes), + finalResidual_(fRes), + noIterations_(nIter), + converged_(converged), + singular_(singular) + {} + + + // Member functions + + //- Return solver name + const word& solverName() const + { + return solverName_; + } + + //- Return solver name + word& solverName() + { + return solverName_; + } + + + //- Return field name + const word& fieldName() const + { + return fieldName_; + } + + + //- Return initial residual + const Type& initialResidual() const + { + return initialResidual_; + } + + //- Return initial residual + Type& initialResidual() + { + return initialResidual_; + } + + + //- Return final residual + const Type& finalResidual() const + { + return finalResidual_; + } + + //- Return final residual + Type& finalResidual() + { + return finalResidual_; + } + + + //- Return number of iterations + label nIterations() const + { + return noIterations_; + } + + //- Return number of iterations + label& nIterations() + { + return noIterations_; + } + + + //- Has the solver converged? + bool converged() const + { + return converged_; + } + + //- Is the matrix singular? + bool singular() const; + + //- Check, store and return convergence + bool checkConvergence + ( + const Type& tolerance, + const Type& relTolerance + ); + + //- Singularity test + bool checkSingularity(const Type& residual); + + //- Print summary of solver performance to the given stream + void print(Ostream& os) const; + + + // Member Operators + + bool operator!=(const SolverPerformance<Type>&) const; + + + // Friend functions + + //- Return the element-wise maximum of two SolverPerformance<Type>s + friend SolverPerformance<Type> max <Type> + ( + const SolverPerformance<Type>&, + const SolverPerformance<Type>& + ); + + + // Ostream Operator + + friend Istream& operator>> <Type> + ( + Istream&, + SolverPerformance<Type>& + ); + + friend Ostream& operator<< <Type> + ( + Ostream&, + const SolverPerformance<Type>& + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#define makeSolverPerformance(Type) \ + \ +typedef Foam::SolverPerformance<Type> \ + solverPerformance##Type; \ + \ +defineNamedTemplateTypeNameAndDebug(solverPerformance##Type, 0); \ + \ +template<> \ +const scalar solverPerformance##Type::great_(1e20); \ + \ +template<> \ +const scalar solverPerformance##Type::small_(1e-20); \ + \ +template<> \ +const scalar solverPerformance##Type::vsmall_(VSMALL); \ + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "SolverPerformance.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.C new file mode 100644 index 0000000000000000000000000000000000000000..5bac12c784c9fd02909584a03377bdd30281c781 --- /dev/null +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.C @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "solverPerformance.H" +#include "fieldTypes.H" + +namespace Foam +{ + makeSolverPerformance(scalar); + makeSolverPerformance(vector); + makeSolverPerformance(sphericalTensor); + makeSolverPerformance(symmTensor); + makeSolverPerformance(tensor); +}; + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.H new file mode 100644 index 0000000000000000000000000000000000000000..139ab128b1c2a8db3ed883c57a3549e82a4b6e63 --- /dev/null +++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Typedef + Foam::solverPerformance + +Description + SolverPerformance instantiated for a scalar + +SourceFiles + solverPerformance.C + +\*---------------------------------------------------------------------------*/ + +#ifndef solverPerformance_H +#define solverPerformance_H + +#include "SolverPerformance.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef SolverPerformance<scalar> solverPerformance; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C index af128d249a8fab2e8617dfa8ee71c52ab60beee9..adbe9113a271c750af850a9443baf3190840d232 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C @@ -55,7 +55,7 @@ void Foam::DiagonalSolver<Type, DType, LUType>::read template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance +Foam::SolverPerformance<Type> Foam::DiagonalSolver<Type, DType, LUType>::solve ( Field<Type>& psi @@ -63,7 +63,7 @@ Foam::DiagonalSolver<Type, DType, LUType>::solve { psi = this->matrix_.source()/this->matrix_.diag(); - return typename LduMatrix<Type, DType, LUType>::solverPerformance + return SolverPerformance<Type> ( typeName, this->fieldName_, diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H index 400891b97776681ec9849d007605785213d9b79f..b9dc537c08afaf1ea2a771046842f84dc776aa56 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H @@ -83,10 +83,7 @@ public: void read(const dictionary& solverDict); //- Solve the matrix with this solver - typename LduMatrix<Type, DType, LUType>::solverPerformance solve - ( - Field<Type>& psi - ) const; + virtual SolverPerformance<Type> solve(Field<Type>& psi) const; }; diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C index 45453ff1d062095d2c3e08f23e467d6887517555..7ee6da05dcf6bd6c733164fe77e85773b9296469 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C @@ -47,7 +47,7 @@ Foam::PBiCCCG<Type, DType, LUType>::PBiCCCG // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance +Foam::SolverPerformance<Type> Foam::PBiCCCG<Type, DType, LUType>::solve ( Field<Type>& psi @@ -56,7 +56,7 @@ Foam::PBiCCCG<Type, DType, LUType>::solve word preconditionerName(this->controlDict_.lookup("preconditioner")); // --- Setup class containing solver performance data - typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf + SolverPerformance<Type> solverPerf ( preconditionerName + typeName, this->fieldName_ @@ -104,7 +104,7 @@ Foam::PBiCCCG<Type, DType, LUType>::solve solverPerf.finalResidual() = solverPerf.initialResidual(); // --- Check convergence, solve if not converged - if (!solverPerf.converged(this->tolerance_, this->relTol_)) + if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) { // --- Select and construct the preconditioner autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner> @@ -183,7 +183,7 @@ Foam::PBiCCCG<Type, DType, LUType>::solve } while ( solverPerf.nIterations()++ < this->maxIter_ - && !(solverPerf.converged(this->tolerance_, this->relTol_)) + && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) ); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H index e5363dc6510a33e8463d09173a7e56d286405143..39c0e56d44768d12e13879c872c6a97035db9fe9 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H @@ -87,10 +87,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve - ( - Field<Type>& psi - ) const; + virtual SolverPerformance<Type> solve(Field<Type>& psi) const; }; diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C index b3665f0617565e4b743c18b326ca1dfabee9fa74..260daedbf0e11c88a2004b7d53a4130b8c3587c6 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C @@ -47,13 +47,13 @@ Foam::PBiCICG<Type, DType, LUType>::PBiCICG // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance +Foam::SolverPerformance<Type> Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const { word preconditionerName(this->controlDict_.lookup("preconditioner")); // --- Setup class containing solver performance data - typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf + SolverPerformance<Type> solverPerf ( preconditionerName + typeName, this->fieldName_ @@ -75,7 +75,7 @@ Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Field<Type> wT(nCells); Type* __restrict__ wTPtr = wT.begin(); - Type wArT = this->matrix_.great_*pTraits<Type>::one; + Type wArT = solverPerf.great_*pTraits<Type>::one; Type wArTold = wArT; // --- Calculate A.psi and T.psi @@ -101,7 +101,7 @@ Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const solverPerf.finalResidual() = solverPerf.initialResidual(); // --- Check convergence, solve if not converged - if (!solverPerf.converged(this->tolerance_, this->relTol_)) + if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) { // --- Select and construct the preconditioner autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner> @@ -137,7 +137,7 @@ Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Type beta = cmptDivide ( wArT, - stabilise(wArTold, this->matrix_.vsmall_) + stabilise(wArTold, solverPerf.vsmall_) ); for (register label cell=0; cell<nCells; cell++) @@ -172,7 +172,7 @@ Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Type alpha = cmptDivide ( wArT, - stabilise(wApT, this->matrix_.vsmall_) + stabilise(wApT, solverPerf.vsmall_) ); for (register label cell=0; cell<nCells; cell++) @@ -188,7 +188,7 @@ Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const } while ( solverPerf.nIterations()++ < this->maxIter_ - && !(solverPerf.converged(this->tolerance_, this->relTol_)) + && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) ); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H index 340ae486678db5e6f1291734753c36b4945200e0..1265af2aac4b6197e9c7e2d7d0038ccb846a4a5a 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H @@ -87,10 +87,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve - ( - Field<Type>& psi - ) const; + virtual SolverPerformance<Type> solve(Field<Type>& psi) const; }; diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C index e6206e61163c47dc88cf1739bba98e4c03f267ba..cc8f1d1c72b1eb0398f6c35dd3775409150e7f1c 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C @@ -47,13 +47,13 @@ Foam::PCICG<Type, DType, LUType>::PCICG // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance +typename Foam::SolverPerformance<Type> Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const { word preconditionerName(this->controlDict_.lookup("preconditioner")); // --- Setup class containing solver performance data - typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf + SolverPerformance<Type> solverPerf ( preconditionerName + typeName, this->fieldName_ @@ -69,7 +69,7 @@ Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Field<Type> wA(nCells); Type* __restrict__ wAPtr = wA.begin(); - Type wArA = this->matrix_.great_*pTraits<Type>::one; + Type wArA = solverPerf.great_*pTraits<Type>::one; Type wArAold = wArA; // --- Calculate A.psi @@ -92,7 +92,7 @@ Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const solverPerf.finalResidual() = solverPerf.initialResidual(); // --- Check convergence, solve if not converged - if (!solverPerf.converged(this->tolerance_, this->relTol_)) + if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) { // --- Select and construct the preconditioner autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner> @@ -126,7 +126,7 @@ Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Type beta = cmptDivide ( wArA, - stabilise(wArAold, this->matrix_.vsmall_) + stabilise(wArAold, solverPerf.vsmall_) ); for (register label cell=0; cell<nCells; cell++) @@ -160,7 +160,7 @@ Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const Type alpha = cmptDivide ( wArA, - stabilise(wApA, this->matrix_.vsmall_) + stabilise(wApA, solverPerf.vsmall_) ); for (register label cell=0; cell<nCells; cell++) @@ -175,7 +175,7 @@ Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const } while ( solverPerf.nIterations()++ < this->maxIter_ - && !(solverPerf.converged(this->tolerance_, this->relTol_)) + && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) ); } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H index 250cc41b9a363dea56420fedfc723bd7ad86b019..cafeb7f61d3a79ae55861e4dea79b68942f83c38 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H @@ -87,10 +87,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve - ( - Field<Type>& psi - ) const; + virtual SolverPerformance<Type> solve(Field<Type>& psi) const; }; diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C index e6f0241308e5c2377c6e38b55f58693f370941bd..9cf5c81223b71386ed7f40a0737045a30766a57d 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C @@ -53,16 +53,16 @@ template<class Type, class DType, class LUType> void Foam::SmoothSolver<Type, DType, LUType>::readControls() { LduMatrix<Type, DType, LUType>::solver::readControls(); - readControl(this->controlDict_, nSweeps_, "nSweeps"); + this->readControl(this->controlDict_, nSweeps_, "nSweeps"); } template<class Type, class DType, class LUType> -typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance +Foam::SolverPerformance<Type> Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const { // --- Setup class containing solver performance data - typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf + SolverPerformance<Type> solverPerf ( typeName, this->fieldName_ @@ -113,7 +113,7 @@ Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const // Check convergence, solve if not converged - if (!solverPerf.converged(this->tolerance_, this->relTol_)) + if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_)) { autoPtr<typename LduMatrix<Type, DType, LUType>::smoother> smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New @@ -141,7 +141,7 @@ Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const } while ( (solverPerf.nIterations() += nSweeps_) < this->maxIter_ - && !(solverPerf.converged(this->tolerance_, this->relTol_)) + && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_)) ); } } diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H index 3d789420b390e3b67710f01fc9d0fba7ed6585fe..54db11be8c569bc80530e486f9b7cada29806ff8 100644 --- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H +++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H @@ -86,10 +86,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve - ( - Field<Type>& psi - ) const; + virtual SolverPerformance<Type> solve(Field<Type>& psi) const; }; diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C index 6fb1b68d0ac816674a0494e73bacb77866d952a6..16dd1a41881b1473634a842772148519571abbe2 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -30,9 +30,6 @@ License defineTypeNameAndDebug(Foam::lduMatrix, 1); -const Foam::scalar Foam::lduMatrix::great_ = 1.0e+20; -const Foam::scalar Foam::lduMatrix::small_ = 1.0e-20; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H index ae7438b937b4c14b924be39e6456466c97969146..ec5d7a89d5df3f7ea81406e1f2e2cf0ff78acf93 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H @@ -57,6 +57,7 @@ SourceFiles #include "typeInfo.H" #include "autoPtr.H" #include "runTimeSelectionTables.H" +#include "solverPerformance.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -86,166 +87,6 @@ class lduMatrix public: - //- Class returned by the solver, containing performance statistics - class solverPerformance - { - word solverName_; - word fieldName_; - scalar initialResidual_; - scalar finalResidual_; - label noIterations_; - bool converged_; - bool singular_; - - - public: - - // Constructors - - //- Construct null - solverPerformance() - : - initialResidual_(0), - finalResidual_(0), - noIterations_(0), - converged_(false), - singular_(false) - {} - - //- Construct from components - solverPerformance - ( - const word& solverName, - const word& fieldName, - const scalar iRes = 0, - const scalar fRes = 0, - const label nIter = 0, - const bool converged = false, - const bool singular = false - ) - : - solverName_(solverName), - fieldName_(fieldName), - initialResidual_(iRes), - finalResidual_(fRes), - noIterations_(nIter), - converged_(converged), - singular_(singular) - {} - - //- Construct from Istream - solverPerformance(Istream&); - - - // Member functions - - //- Return solver name - const word& solverName() const - { - return solverName_; - } - - //- Return solver name - word& solverName() - { - return solverName_; - } - - - //- Return field name - const word& fieldName() const - { - return fieldName_; - } - - - //- Return initial residual - scalar initialResidual() const - { - return initialResidual_; - } - - //- Return initial residual - scalar& initialResidual() - { - return initialResidual_; - } - - - //- Return final residual - scalar finalResidual() const - { - return finalResidual_; - } - - //- Return final residual - scalar& finalResidual() - { - return finalResidual_; - } - - - //- Return number of iterations - label nIterations() const - { - return noIterations_; - } - - //- Return number of iterations - label& nIterations() - { - return noIterations_; - } - - - //- Has the solver converged? - bool converged() const - { - return converged_; - } - - //- Is the matrix singular? - bool singular() const - { - return singular_; - } - - //- Check, store and return convergence - bool checkConvergence - ( - const scalar tolerance, - const scalar relTolerance - ); - - //- Singularity test - bool checkSingularity(const scalar residual); - - //- Print summary of solver performance - void print() const; - - - // Member Operators - - bool operator!=(const solverPerformance&) const; - - - // Friend functions - - //- Return the element-wise maximum of two solverPerformances - friend solverPerformance max - ( - const solverPerformance&, - const solverPerformance& - ); - - - // Ostream Operator - - friend Istream& operator>>(Istream&, solverPerformance&); - friend Ostream& operator<<(Ostream&, const solverPerformance&); - }; - - //- Abstract base-class for lduMatrix solvers class solver { @@ -670,12 +511,6 @@ public: // Declare name of the class and its debug switch ClassName("lduMatrix"); - //- Large scalar for the use in solvers - static const scalar great_; - - //- Small scalar for the use in solvers - static const scalar small_; - // Constructors diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C index db2b6d8bbb4f077a306466c040b6ce61405ecf19..723a99980646532ea78c7739fe89efa10204a8ec 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -295,4 +295,36 @@ Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual } +Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const +{ + tmp<scalarField > tH1 + ( + new scalarField(lduAddr().size(), 0.0) + ); + + if (lowerPtr_ || upperPtr_) + { + scalarField& H1_ = tH1(); + + scalar* __restrict__ H1Ptr = H1_.begin(); + + const label* __restrict__ uPtr = lduAddr().upperAddr().begin(); + const label* __restrict__ lPtr = lduAddr().lowerAddr().begin(); + + const scalar* __restrict__ lowerPtr = lower().begin(); + const scalar* __restrict__ upperPtr = upper().begin(); + + register const label nFaces = upper().size(); + + for (register label face=0; face<nFaces; face++) + { + H1Ptr[uPtr[face]] -= lowerPtr[face]; + H1Ptr[lPtr[face]] -= upperPtr[face]; + } + } + + return tH1; +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C index 07260b4a91009fdbdc7265ab6ad6765f75183c40..df1ba0d9455eaa0ef1b35dcf6bd8b9571b5cdbb4 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C @@ -351,36 +351,4 @@ void Foam::lduMatrix::operator*=(scalar s) } -Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const -{ - tmp<scalarField > tH1 - ( - new scalarField(lduAddr().size(), 0.0) - ); - - if (lowerPtr_ || upperPtr_) - { - scalarField& H1_ = tH1(); - - scalar* __restrict__ H1Ptr = H1_.begin(); - - const label* __restrict__ uPtr = lduAddr().upperAddr().begin(); - const label* __restrict__ lPtr = lduAddr().lowerAddr().begin(); - - const scalar* __restrict__ lowerPtr = lower().begin(); - const scalar* __restrict__ upperPtr = upper().begin(); - - register const label nFaces = upper().size(); - - for (register label face=0; face<nFaces; face++) - { - H1Ptr[uPtr[face]] -= lowerPtr[face]; - H1Ptr[lPtr[face]] -= upperPtr[face]; - } - } - - return tH1; -} - - // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C index 0381db1b89f836177634d26f044fd9e4621aee65..bedb55529eeef8ed05289d11343adece8b0d28f4 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -188,10 +188,12 @@ Foam::scalar Foam::lduMatrix::solver::normFactor matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_); tmpField *= gAverage(psi); - return gSum(mag(Apsi - tmpField) + mag(source - tmpField)) + matrix_.small_; + return + gSum(mag(Apsi - tmpField) + mag(source - tmpField)) + + solverPerformance::small_; // At convergence this simpler method is equivalent to the above - // return 2*gSumMag(source) + matrix_.small_; + // return 2*gSumMag(source) + solverPerformance::small_; } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTests.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTests.C deleted file mode 100644 index e9294d6a22729d321763413438c7280cff9b3769..0000000000000000000000000000000000000000 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTests.C +++ /dev/null @@ -1,190 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. - -Description - Convergence and singularity tests for solvers. - -\*---------------------------------------------------------------------------*/ - -#include "lduMatrix.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -Foam::lduMatrix::solverPerformance::solverPerformance(Istream& is) -{ - is >> *this; -} - - -bool Foam::lduMatrix::solverPerformance::checkConvergence -( - const scalar Tolerance, - const scalar RelTolerance -) -{ - if (debug >= 2) - { - Info<< solverName_ - << ": Iteration " << noIterations_ - << " residual = " << finalResidual_ - << endl; - } - - if - ( - finalResidual_ < Tolerance - || ( - RelTolerance > SMALL - && finalResidual_ <= RelTolerance*initialResidual_ - ) - // || (solverName == "symSolve" && iter == 0) - ) - { - converged_ = true; - } - else - { - converged_ = false; - } - - return converged_; -} - - -bool Foam::lduMatrix::solverPerformance::checkSingularity -( - const scalar residual -) -{ - if (residual > VSMALL) - { - singular_ = false; - } - else - { - singular_ = true; - } - - return singular_; -} - - -void Foam::lduMatrix::solverPerformance::print() const -{ - if (debug) - { - Info<< solverName_ << ": Solving for " << fieldName_; - - if (singular()) - { - Info<< ": solution singularity" << endl; - } - else - { - Info<< ", Initial residual = " << initialResidual_ - << ", Final residual = " << finalResidual_ - << ", No Iterations " << noIterations_ - << endl; - } - } -} - - -bool Foam::lduMatrix::solverPerformance::operator!= -( - const lduMatrix::solverPerformance& sp -) const -{ - return - ( - solverName() != sp.solverName() - || fieldName() != sp.fieldName() - || initialResidual() != sp.initialResidual() - || finalResidual() != sp.finalResidual() - || nIterations() != sp.nIterations() - || converged() != sp.converged() - || singular() != sp.singular() - ); -} - - -Foam::lduMatrix::solverPerformance Foam::max -( - const lduMatrix::solverPerformance& sp1, - const lduMatrix::solverPerformance& sp2 -) -{ - return lduMatrix::solverPerformance - ( - sp1.solverName(), - sp1.fieldName_, - max(sp1.initialResidual(), sp2.initialResidual()), - max(sp1.finalResidual(), sp2.finalResidual()), - max(sp1.nIterations(), sp2.nIterations()), - sp1.converged() && sp2.converged(), - sp1.singular() || sp2.singular() - ); -} - - -Foam::Istream& Foam::operator>> -( - Istream& is, - Foam::lduMatrix::solverPerformance& sp -) -{ - is.readBeginList("lduMatrix::solverPerformance"); - is >> sp.solverName_ - >> sp.fieldName_ - >> sp.initialResidual_ - >> sp.finalResidual_ - >> sp.noIterations_ - >> sp.converged_ - >> sp.singular_; - is.readEndList("lduMatrix::solverPerformance"); - - return is; -} - - -Foam::Ostream& Foam::operator<< -( - Ostream& os, - const Foam::lduMatrix::solverPerformance& sp -) -{ - os << token::BEGIN_LIST - << sp.solverName_ << token::SPACE - << sp.fieldName_ << token::SPACE - << sp.initialResidual_ << token::SPACE - << sp.finalResidual_ << token::SPACE - << sp.noIterations_ << token::SPACE - << sp.converged_ << token::SPACE - << sp.singular_ << token::SPACE - << token::END_LIST; - - return os; -} - - -// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H index 27094885e6fcdc9405eed78bce0b63eeaa639bc6..273d2c75b2171803ea82f515dfb84bb1d09d3534 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -229,7 +229,7 @@ public: // Member Functions //- Solve - virtual lduMatrix::solverPerformance solve + virtual solverPerformance solve ( scalarField& psi, const scalarField& source, diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C index e0574e5895390b27160cd59e26afb7fc5595fcfb..248c421f2e347dbf4a20cd3316a3c8395412b255 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -30,7 +30,7 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve +Foam::solverPerformance Foam::GAMGSolver::solve ( scalarField& psi, const scalarField& source, @@ -38,7 +38,7 @@ Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve ) const { // Setup class containing solver performance data - lduMatrix::solverPerformance solverPerf(typeName, fieldName_); + solverPerformance solverPerf(typeName, fieldName_); // Calculate A.psi used to calculate the initial residual scalarField Apsi(psi.size()); @@ -103,7 +103,7 @@ Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve if (debug >= 2) { - solverPerf.print(); + solverPerf.print(Info); } } while ( @@ -429,7 +429,7 @@ void Foam::GAMGSolver::solveCoarsestLevel { const label coarsestLevel = matrixLevels_.size() - 1; coarsestCorrField = 0; - lduMatrix::solverPerformance coarseSolverPerf; + solverPerformance coarseSolverPerf; if (matrixLevels_[coarsestLevel].asymmetric()) { @@ -468,7 +468,7 @@ void Foam::GAMGSolver::solveCoarsestLevel if (debug >= 2) { - coarseSolverPerf.print(); + coarseSolverPerf.print(Info); } } } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C index d47b606c242ba5e76a390e8fee290d19afb104fe..8f224f85e743fa19ef06184b0dfa9e812d85e830 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C @@ -62,7 +62,7 @@ Foam::PBiCG::PBiCG // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::lduMatrix::solverPerformance Foam::PBiCG::solve +Foam::solverPerformance Foam::PBiCG::solve ( scalarField& psi, const scalarField& source, @@ -70,7 +70,7 @@ Foam::lduMatrix::solverPerformance Foam::PBiCG::solve ) const { // --- Setup class containing solver performance data - lduMatrix::solverPerformance solverPerf + solverPerformance solverPerf ( lduMatrix::preconditioner::getName(controlDict_) + typeName, fieldName_ @@ -92,7 +92,7 @@ Foam::lduMatrix::solverPerformance Foam::PBiCG::solve scalarField wT(nCells); scalar* __restrict__ wTPtr = wT.begin(); - scalar wArT = matrix_.great_; + scalar wArT = solverPerf.great_; scalar wArTold = wArT; // --- Calculate A.psi and T.psi diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.H b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.H index 7479f75cc5533ea0e12f8509bfdc3d4edeac6678..fda7d7d314d895c938ef360111d0bf505ff6eaff 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -88,7 +88,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual lduMatrix::solverPerformance solve + virtual solverPerformance solve ( scalarField& psi, const scalarField& source, diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C index 37233afb719aee9470dad7bc9fcc767862212141..b4a06ec8a3df10db2a9c074c27e319b8f10be584 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -62,7 +62,7 @@ Foam::PCG::PCG // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::lduMatrix::solverPerformance Foam::PCG::solve +Foam::solverPerformance Foam::PCG::solve ( scalarField& psi, const scalarField& source, @@ -70,7 +70,7 @@ Foam::lduMatrix::solverPerformance Foam::PCG::solve ) const { // --- Setup class containing solver performance data - lduMatrix::solverPerformance solverPerf + solverPerformance solverPerf ( lduMatrix::preconditioner::getName(controlDict_) + typeName, fieldName_ @@ -86,7 +86,7 @@ Foam::lduMatrix::solverPerformance Foam::PCG::solve scalarField wA(nCells); scalar* __restrict__ wAPtr = wA.begin(); - scalar wArA = matrix_.great_; + scalar wArA = solverPerf.great_; scalar wArAold = wArA; // --- Calculate A.psi diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H index 7aca6131848fc743c9890a5f8ca676a861ec0146..74b3abd17c8fecbc05fd9e77cb8ba4e9be5465b4 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -88,7 +88,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual lduMatrix::solverPerformance solve + virtual solverPerformance solve ( scalarField& psi, const scalarField& source, diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C index fd03af39e7a6b3c1d1d710d3b4807780a6e85713..73e386f5eddbf34cdafc8cd4b35be1ef7e539c87 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -56,7 +56,7 @@ Foam::diagonalSolver::diagonalSolver // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -Foam::lduMatrix::solverPerformance Foam::diagonalSolver::solve +Foam::solverPerformance Foam::diagonalSolver::solve ( scalarField& psi, const scalarField& source, @@ -65,7 +65,7 @@ Foam::lduMatrix::solverPerformance Foam::diagonalSolver::solve { psi = source/matrix_.diag(); - return lduMatrix::solverPerformance + return solverPerformance ( typeName, fieldName_, diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H index 1d67155e8cbe3e51d90fe9944285ad037c743a42..85908b33496b8a0bf45612c47248a6aa23b85f4b 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -86,7 +86,7 @@ public: {} //- Solve the matrix with this solver - lduMatrix::solverPerformance solve + solverPerformance solve ( scalarField& psi, const scalarField& source, diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C index 05597a3b66b6780f90ef1f1975503665573155bb..4287eb56999feeb249a9d3dfe544b4323669124d 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,7 +74,7 @@ void Foam::smoothSolver::readControls() } -Foam::lduMatrix::solverPerformance Foam::smoothSolver::solve +Foam::solverPerformance Foam::smoothSolver::solve ( scalarField& psi, const scalarField& source, @@ -82,7 +82,7 @@ Foam::lduMatrix::solverPerformance Foam::smoothSolver::solve ) const { // Setup class containing solver performance data - lduMatrix::solverPerformance solverPerf(typeName, fieldName_); + solverPerformance solverPerf(typeName, fieldName_); // If the nSweeps_ is negative do a fixed number of sweeps if (nSweeps_ < 0) diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H index ec62305803c2aba27c9cdd186d68ea3d3e60bbd2..2009768401218f28214dc3b42925f0ba1fd9bef8 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -92,7 +92,7 @@ public: // Member Functions //- Solve the matrix with this solver - virtual lduMatrix::solverPerformance solve + virtual solverPerformance solve ( scalarField& psi, const scalarField& source, diff --git a/src/OpenFOAM/meshes/data/data.C b/src/OpenFOAM/meshes/data/data.C index 7540482cfe17a741a7b4bcb487a72b6d2246e036..02366479aa7bbc4c4a3c5bf1dbba6abab703e55d 100644 --- a/src/OpenFOAM/meshes/data/data.C +++ b/src/OpenFOAM/meshes/data/data.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,7 @@ License #include "data.H" #include "Time.H" -#include "lduMatrix.H" +#include "solverPerformance.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -63,12 +63,12 @@ const Foam::dictionary& Foam::data::solverPerformanceDict() const void Foam::data::setSolverPerformance ( const word& name, - const lduMatrix::solverPerformance& sp + const solverPerformance& sp ) const { dictionary& dict = const_cast<dictionary&>(solverPerformanceDict()); - List<lduMatrix::solverPerformance> perfs; + List<solverPerformance> perfs; if (prevTimeIndex_ != this->time().timeIndex()) { @@ -90,7 +90,7 @@ void Foam::data::setSolverPerformance void Foam::data::setSolverPerformance ( - const lduMatrix::solverPerformance& sp + const solverPerformance& sp ) const { setSolverPerformance(sp.fieldName(), sp); diff --git a/src/OpenFOAM/meshes/data/data.H b/src/OpenFOAM/meshes/data/data.H index 31ce391acbac39c9c0d8be84926f36684029787a..e4fd40fb6a1041efb036ce4ab3112770828bb996 100644 --- a/src/OpenFOAM/meshes/data/data.H +++ b/src/OpenFOAM/meshes/data/data.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -39,7 +39,7 @@ SourceFiles #define data_H #include "IOdictionary.H" -#include "lduMatrix.H" +#include "solverPerformance.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -94,13 +94,13 @@ public: void setSolverPerformance ( const word& name, - const lduMatrix::solverPerformance& + const solverPerformance& ) const; //- Add/set the solverPerformance entry, using its fieldName void setSolverPerformance ( - const lduMatrix::solverPerformance& + const solverPerformance& ) const; }; diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C index e59c5ec898d3427c074141e7404c7d71dcb9937d..4b52d26699939fb872ab7e68244303b7f1fc32b9 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C +++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C @@ -45,6 +45,7 @@ const Foam::interpolationWeights& Foam::TableBase<Type>::interpolator() const tableSamples_ ); } + return interpolatorPtr_(); } @@ -172,6 +173,7 @@ Foam::TableBase<Type>::outOfBounds { boundsHandling prev = boundsHandling_; boundsHandling_ = bound; + return prev; } @@ -366,22 +368,8 @@ Type Foam::TableBase<Type>::value(const scalar x) const { t += currentWeights_[i]*table_[currentIndices_[i]].second(); } - return t; - //// Find i such that x(i) < xDash < x(i+1) - //label i = 0; - //while ((table_[i+1].first() < xDash) && (i+1 < table_.size())) - //{ - // i++; - //} - // - //// Linear interpolation to find value - //return Type - //( - // (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first()) - // * (table_[i+1].second() - table_[i].second()) - // + table_[i].second() - //); + return t; } @@ -396,66 +384,8 @@ Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const { sum += currentWeights_[i]*table_[currentIndices_[i]].second(); } - return sum; - - //// Initialise return value - //Type sum = pTraits<Type>::zero; - // - //// Return zero if out of bounds - //if ((x1 > table_.last().first()) || (x2 < table_[0].first())) - //{ - // return sum; - //} - // - //// Find next index greater than x1 - //label id1 = 0; - //while ((table_[id1].first() < x1) && (id1 < table_.size())) - //{ - // id1++; - //} - // - //// Find next index less than x2 - //label id2 = table_.size() - 1; - //while ((table_[id2].first() > x2) && (id2 >= 1)) - //{ - // id2--; - //} - // - //if ((id1 - id2) == 1) - //{ - // // x1 and x2 lie within 1 interval - // sum = 0.5*(value(x1) + value(x2))*(x2 - x1); - //} - //else - //{ - // // x1 and x2 cross multiple intervals - // - // // Integrate table body - // for (label i=id1; i<id2; i++) - // { - // sum += - // (table_[i].second() + table_[i+1].second()) - // * (table_[i+1].first() - table_[i].first()); - // } - // sum *= 0.5; - // - // // Add table ends (partial segments) - // if (id1 > 0) - // { - // sum += 0.5 - // * (value(x1) + table_[id1].second()) - // * (table_[id1].first() - x1); - // } - // if (id2 < table_.size() - 1) - // { - // sum += 0.5 - // * (table_[id2].second() + value(x2)) - // * (x2 - table_[id2].first()); - // } - //} - // - //return sum; + return sum; } @@ -481,9 +411,9 @@ Foam::dimensioned<Type> Foam::TableBase<Type>::dimIntegrate ); } + // * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * // #include "TableBaseIO.C" - // ************************************************************************* // diff --git a/src/conversion/ensight/part/ensightPartCells.C b/src/conversion/ensight/part/ensightPartCells.C index 4c44abedcb1ba97f2dca3d861df1d6b67dfd2385..80f0cdc4ce61126c71f9f35b7a422fcac139fcfe 100644 --- a/src/conversion/ensight/part/ensightPartCells.C +++ b/src/conversion/ensight/part/ensightPartCells.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -313,6 +313,7 @@ void Foam::ensightPartCells::writeConnectivity if (key == "nfaced") { const faceList& meshFaces = mesh_.faces(); + const labelUList& owner = mesh_.faceOwner(); // write the number of faces per element forAll(idList, i) @@ -345,16 +346,33 @@ void Foam::ensightPartCells::writeConnectivity const label id = idList[i] + offset_; const labelUList& cFace = mesh_.cells()[id]; - forAll(cFace, faceI) + forAll(cFace, cFaceI) { - const face& cf = meshFaces[cFace[faceI]]; + const label faceId = cFace[cFaceI]; + const face& cf = meshFaces[faceId]; + + // convert global -> local index + // (note: Ensight indices start with 1) - forAll(cf, ptI) + // ensight >= 9 needs consistently oriented nfaced cells + if (id == owner[faceId]) { - // convert global -> local index - // (note: Ensight indices start with 1) - os.write(pointMap[cf[ptI]] + 1); + forAll(cf, ptI) + { + os.write(pointMap[cf[ptI]] + 1); + } } + else + { + // as per face::reverseFace(), but without copying + + os.write(pointMap[cf[0]] + 1); + for (label ptI = cf.size()-1; ptI > 0; --ptI) + { + os.write(pointMap[cf[ptI]] + 1); + } + } + os.newline(); } } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C index 415c6facdca40204dfe3e097e10452ec756a9933..22872abf47d6b1c5f87150831ef4ce1d22a82b3a 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C @@ -27,6 +27,10 @@ License #include "polyMesh.H" #include "polyTopoChange.H" #include "ListOps.H" +#include "globalMeshData.H" +#include "OFstream.H" +#include "meshTools.H" +#include "syncTools.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -121,6 +125,7 @@ void Foam::edgeCollapser::filterFace(const label faceI, face& f) const } } + // Check for pinched face. Tries to correct // - consecutive duplicate vertex. Removes duplicate vertex. // - duplicate vertex with one other vertex in between (spike). @@ -190,15 +195,15 @@ void Foam::edgeCollapser::printRegions() const if (master != -1) { - Info<< "Region:" << regionI << nl + Pout<< "Region:" << regionI << nl << " master:" << master - << ' ' << mesh_.points()[master] << nl; + << ' ' << pointRegionMasterLocation_[regionI] << nl; forAll(pointRegion_, pointI) { if (pointRegion_[pointI] == regionI && pointI != master) { - Info<< " slave:" << pointI + Pout<< " slave:" << pointI << ' ' << mesh_.points()[pointI] << nl; } } @@ -272,6 +277,7 @@ Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh) : mesh_(mesh), pointRegion_(mesh.nPoints(), -1), + pointRegionMasterLocation_(mesh.nPoints() / 100), pointRegionMaster_(mesh.nPoints() / 100), freeRegions_() {} @@ -289,6 +295,8 @@ bool Foam::edgeCollapser::unaffectedEdge(const label edgeI) const bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master) { + const pointField& points = mesh_.points(); + const edge& e = mesh_.edges()[edgeI]; label pointRegion0 = pointRegion_[e[0]]; @@ -310,7 +318,7 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master) { FatalErrorIn ("edgeCollapser::collapseEdge(const label, const label)") - << "Problem : freeed region :" << freeRegion + << "Problem : freed region :" << freeRegion << " has already master " << pointRegionMaster_[freeRegion] << abort(FatalError); @@ -327,13 +335,22 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master) pointRegion_[e[1]] = freeRegion; pointRegionMaster_(freeRegion) = master; + pointRegionMasterLocation_(freeRegion) = points[master]; } else { // e[1] is part of collapse network, e[0] not. Add e0 to e1 region. pointRegion_[e[0]] = pointRegion1; - pointRegionMaster_[pointRegion1] = master; + if + ( + pointRegionMaster_[pointRegion1] == e[0] + || pointRegionMaster_[pointRegion1] == e[1] + ) + { + pointRegionMaster_[pointRegion1] = master; + pointRegionMasterLocation_[pointRegion1] = points[master]; + } } } else @@ -343,7 +360,15 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master) // e[0] is part of collapse network. Add e1 to e0 region pointRegion_[e[1]] = pointRegion0; - pointRegionMaster_[pointRegion0] = master; + if + ( + pointRegionMaster_[pointRegion0] == e[0] + || pointRegionMaster_[pointRegion0] == e[1] + ) + { + pointRegionMaster_[pointRegion0] = master; + pointRegionMasterLocation_[pointRegion0] = points[master]; + } } else if (pointRegion0 != pointRegion1) { @@ -356,6 +381,9 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master) // Use minRegion as region for combined net, free maxRegion. pointRegionMaster_[minRegion] = master; pointRegionMaster_[maxRegion] = -1; + pointRegionMasterLocation_[minRegion] = points[master]; + pointRegionMasterLocation_[maxRegion] = point(0, 0, 0); + freeRegions_.insert(maxRegion); if (minRegion != pointRegion0) @@ -380,12 +408,61 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod) const labelList& faceNeighbour = mesh_.faceNeighbour(); const labelListList& pointFaces = mesh_.pointFaces(); const labelListList& cellEdges = mesh_.cellEdges(); + const pointZoneMesh& pointZones = mesh_.pointZones(); + - // Print regions: - //printRegions() bool meshChanged = false; + // Synchronise pointRegionMasterLocation_ + const globalMeshData& globalData = mesh_.globalData(); + const mapDistribute& map = globalData.globalPointSlavesMap(); + const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch(); + const labelList& meshPoints = coupledPatch.meshPoints(); + const Map<label>& meshPointMap = coupledPatch.meshPointMap(); + + + List<point> newPoints = coupledPatch.localPoints(); + + for (label pI = 0; pI < coupledPatch.nPoints(); ++pI) + { + const label pointRegionMaster = pointRegion_[meshPoints[pI]]; + + if (pointRegionMaster != -1) + { + newPoints[pI] + = pointRegionMasterLocation_[pointRegionMaster]; + } + } + + globalData.syncData + ( + newPoints, + globalData.globalPointSlaves(), + globalData.globalPointTransformedSlaves(), + map, + minMagSqrEqOp<point>() + ); + + OFstream str1("newPoints_" + name(Pstream::myProcNo()) + ".obj"); + forAll(pointRegion_, pI) + { + if (meshPointMap.found(pI)) + { + meshTools::writeOBJ(str1, newPoints[meshPointMap[pI]]); + } + } + + for (label pI = 0; pI < coupledPatch.nPoints(); ++pI) + { + const label pointRegionMaster = pointRegion_[meshPoints[pI]]; + + if (pointRegionMaster != -1) + { + pointRegionMasterLocation_[pointRegionMaster] + = newPoints[pI]; + } + } // Current faces (is also collapseStatus: f.size() < 3) faceList newFaces(mesh_.faces()); @@ -393,7 +470,6 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod) // Current cellCollapse status boolList cellRemoved(mesh_.nCells(), false); - do { // Update face collapse from edge collapses @@ -521,12 +597,49 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod) } } + // Modify the point location of the remaining points + forAll(pointRegion_, pointI) + { + const label pointRegion = pointRegion_[pointI]; + + if + ( + !pointRemoved(pointI) + && meshPointMap.found(pointI) + ) + { + meshMod.modifyPoint + ( + pointI, + newPoints[meshPointMap[pointI]], + pointZones.whichZone(pointI), + false + ); + } + else if + ( + pointRegion != -1 + && !pointRemoved(pointI) + && !meshPointMap.found(pointI) + ) + { + const point& collapsePoint + = pointRegionMasterLocation_[pointRegion]; + + meshMod.modifyPoint + ( + pointI, + collapsePoint, + pointZones.whichZone(pointI), + false + ); + } + } const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh(); const faceZoneMesh& faceZones = mesh_.faceZones(); - // Renumber faces that use points forAll(pointRegion_, pointI) { @@ -585,6 +698,9 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod) } } + // Print regions: +// printRegions(); + return meshChanged; } @@ -593,8 +709,10 @@ void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map) { pointRegion_.setSize(mesh_.nPoints()); pointRegion_ = -1; + // Reset count, do not remove underlying storage pointRegionMaster_.clear(); + pointRegionMasterLocation_.clear(); freeRegions_.clear(); } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H index 98d8ef30d5bce148647499c65cbcbde299df9fc1..3a10a9319749081af5960967cb6eaf09ec3e17d4 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H @@ -39,6 +39,7 @@ SourceFiles #include "labelList.H" #include "DynamicList.H" +#include "point.H" #include "typeInfo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -66,6 +67,10 @@ class edgeCollapser //- For every point -1 or region number labelList pointRegion_; + //- Actual location of the point to collapse to for every region master + // point. This will be forced to be consistent across processors + DynamicList<point> pointRegionMasterLocation_; + //- -1 or master vertex for region number DynamicList<label> pointRegionMaster_; diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C index 827913a1c129e97a45871b7086e86e695737c633..261945476db5ca2e746cd78afa44b49b14591d61 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C @@ -42,6 +42,7 @@ License #include "mapDistributePolyMesh.H" #include "refinementData.H" #include "refinementDistanceData.H" +#include "degenerateMatcher.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -1750,6 +1751,184 @@ void Foam::hexRef8::setInstance(const fileName& inst) } +void Foam::hexRef8::collectLevelPoints +( + const labelList& f, + const label level, + DynamicList<label>& points +) const +{ + forAll(f, fp) + { + if (pointLevel_[f[fp]] <= level) + { + points.append(f[fp]); + } + } +} + + +void Foam::hexRef8::collectLevelPoints +( + const labelList& meshPoints, + const labelList& f, + const label level, + DynamicList<label>& points +) const +{ + forAll(f, fp) + { + label pointI = meshPoints[f[fp]]; + if (pointLevel_[pointI] <= level) + { + points.append(pointI); + } + } +} + + +// Return true if we've found 6 quads. faces guaranteed to be outwards pointing. +bool Foam::hexRef8::matchHexShape +( + const label cellI, + const label cellLevel, + DynamicList<face>& quads +) const +{ + const cell& cFaces = mesh_.cells()[cellI]; + + // Work arrays + DynamicList<label> verts(4); + quads.clear(); + + + // 1. pick up any faces with four cellLevel points + + forAll(cFaces, i) + { + label faceI = cFaces[i]; + const face& f = mesh_.faces()[faceI]; + + verts.clear(); + collectLevelPoints(f, cellLevel, verts); + if (verts.size() == 4) + { + if (mesh_.faceOwner()[faceI] != cellI) + { + reverse(verts); + } + quads.append(face(0)); + labelList& quadVerts = quads.last(); + quadVerts.transfer(verts); + } + } + + + if (quads.size() < 6) + { + Map<labelList> pointFaces(2*cFaces.size()); + + forAll(cFaces, i) + { + label faceI = cFaces[i]; + const face& f = mesh_.faces()[faceI]; + + // Pick up any faces with only one level point. + // See if there are four of these where the commont point + // is a level+1 point. This common point is then the mid of + // a split face. + + verts.clear(); + collectLevelPoints(f, cellLevel, verts); + if (verts.size() == 1) + { + // Add to pointFaces for any level+1 point (this might be + // a midpoint of a split face) + forAll(f, fp) + { + label pointI = f[fp]; + if (pointLevel_[pointI] == cellLevel+1) + { + Map<labelList>::iterator iter = + pointFaces.find(pointI); + if (iter != pointFaces.end()) + { + labelList& pFaces = iter(); + if (findIndex(pFaces, faceI) == -1) + { + pFaces.append(faceI); + } + } + else + { + pointFaces.insert + ( + pointI, + labelList(1, faceI) + ); + } + } + } + } + } + + // 2. Check if we've collected any midPoints. + forAllConstIter(Map<labelList>, pointFaces, iter) + { + const labelList& pFaces = iter(); + + if (pFaces.size() == 4) + { + // Collect and orient. + faceList fourFaces(pFaces.size()); + forAll(pFaces, pFaceI) + { + label faceI = pFaces[pFaceI]; + const face& f = mesh_.faces()[faceI]; + if (mesh_.faceOwner()[faceI] == cellI) + { + fourFaces[pFaceI] = f; + } + else + { + fourFaces[pFaceI] = f.reverseFace(); + } + } + + primitivePatch bigFace + ( + SubList<face>(fourFaces, fourFaces.size()), + mesh_.points() + ); + const labelListList& edgeLoops = bigFace.edgeLoops(); + + if (edgeLoops.size() == 1) + { + // Collect the 4 cellLevel points + verts.clear(); + collectLevelPoints + ( + bigFace.meshPoints(), + bigFace.edgeLoops()[0], + cellLevel, + verts + ); + + if (verts.size() == 4) + { + quads.append(face(0)); + labelList& quadVerts = quads.last(); + quadVerts.transfer(verts); + } + } + } + } + } + + return (quads.size() == 6); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from mesh, read refinement data @@ -4337,6 +4516,9 @@ void Foam::hexRef8::updateMesh // Update face removal engine faceRemover_.updateMesh(map); + + // Clear cell shapes + cellShapesPtr_.clear(); } @@ -4421,6 +4603,9 @@ void Foam::hexRef8::subset // Nothing needs doing to faceRemover. //faceRemover_.subset(pointMap, faceMap, cellMap); + + // Clear cell shapes + cellShapesPtr_.clear(); } @@ -4447,6 +4632,9 @@ void Foam::hexRef8::distribute(const mapDistributePolyMesh& map) // Update face removal engine faceRemover_.distribute(map); + + // Clear cell shapes + cellShapesPtr_.clear(); } @@ -4922,6 +5110,66 @@ void Foam::hexRef8::checkRefinementLevels } +const Foam::cellShapeList& Foam::hexRef8::cellShapes() const +{ + if (cellShapesPtr_.empty()) + { + if (debug) + { + Pout<< "hexRef8::cellShapes() : calculating splitHex cellShapes." + << " cellLevel:" << cellLevel_.size() + << " pointLevel:" << pointLevel_.size() + << endl; + } + + const cellShapeList& meshShapes = mesh_.cellShapes(); + cellShapesPtr_.reset(new cellShapeList(meshShapes)); + + label nSplitHex = 0; + label nUnrecognised = 0; + + forAll(cellLevel_, cellI) + { + if (meshShapes[cellI].model().index() == 0) + { + label level = cellLevel_[cellI]; + + // Return true if we've found 6 quads + DynamicList<face> quads; + bool haveQuads = matchHexShape + ( + cellI, + level, + quads + ); + + if (haveQuads) + { + faceList faces(quads.xfer()); + cellShapesPtr_()[cellI] = degenerateMatcher::match(faces); + nSplitHex++; + } + else + { + nUnrecognised++; + } + } + } + if (debug) + { + Pout<< "hexRef8::cellShapes() :" + << " nCells:" << mesh_.nCells() << " of which" << nl + << " primitive:" << (mesh_.nCells()-nSplitHex-nUnrecognised) + << nl + << " split-hex:" << nSplitHex << nl + << " poly :" << nUnrecognised << nl + << endl; + } + } + return cellShapesPtr_(); +} + + // // Unrefinement diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H index 503412dd38cd7812e8bfda8faccb570ff0f32fe1..31e2a5587dda72626ce01d96ffaa7e1de70d6a78 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H @@ -44,6 +44,7 @@ SourceFiles #include "refinementHistory.H" #include "PackedBoolList.H" #include "uniformDimensionedFields.H" +#include "cellShapeList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -89,6 +90,9 @@ class hexRef8 //- Level of saved cells Map<label> savedCellLevel_; + //- cell shapes when seen as split hexes + mutable autoPtr<cellShapeList> cellShapesPtr_; + // Private Member Functions @@ -303,6 +307,34 @@ class hexRef8 void checkWantedRefinementLevels(const labelList&) const; + // Cellshape recognition + + //- Collect all points on face of certain level + void collectLevelPoints + ( + const labelList& f, + const label level, + DynamicList<label>& points + ) const; + + //- Collect all points on face (in local numbering) of certain level + void collectLevelPoints + ( + const labelList& meshPoints, + const labelList& f, + const label level, + DynamicList<label>& points + ) const; + + //- Collect all faces with four corner points and return true if + // hex was matched (6 faces of each four corner points) + bool matchHexShape + ( + const label cellI, + const label cellLevel, + DynamicList<face>& quads + ) const; + //- Disallow default bitwise copy construct hexRef8(const hexRef8&); @@ -497,6 +529,10 @@ public: const labelList& pointsToCheck ) const; + //- Utility: get hexes as cell shapes + const cellShapeList& cellShapes() const; + + // Unrefinement (undoing refinement, not arbitrary coarsening) //- Return the points at the centre of top-level split cells diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C index e8948c4faded46d41c48759e29e6935f0cf7f562..bdbd851527f9ef90334029cb0b4e489c0c89482e 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C @@ -696,16 +696,16 @@ void Foam::tetDecomposer::setRefinement void Foam::tetDecomposer::updateMesh(const mapPolyMesh& map) { - inplaceRenumber(map.pointMap(), cellToPoint_); - inplaceRenumber(map.pointMap(), faceToPoint_); + inplaceRenumber(map.reversePointMap(), cellToPoint_); + inplaceRenumber(map.reversePointMap(), faceToPoint_); forAll(faceOwnerCells_, faceI) { - inplaceRenumber(map.cellMap(), faceOwnerCells_[faceI]); + inplaceRenumber(map.reverseCellMap(), faceOwnerCells_[faceI]); } forAll(faceNeighbourCells_, faceI) { - inplaceRenumber(map.cellMap(), faceNeighbourCells_[faceI]); + inplaceRenumber(map.reverseCellMap(), faceNeighbourCells_[faceI]); } } diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C b/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C index 865ed3ab62ec7459c76ad59b6e7bb1d417c233f1..fb695078255d7f5015a3c85e9ce96a80473914ab 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C +++ b/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -48,7 +48,7 @@ Foam::scalar Foam::seriesProfile::evaluate forAll(values, i) { - result += values[i]*cos((i+1)*xIn); + result += values[i]*sin((i + 1)*xIn); } return result; diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H b/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H index 73b474aae1f3bb22bbdbacc1b29ea35a914051fb..7c4793e28e9fbe954dd6bf72af0a22c314038665 100644 --- a/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H +++ b/src/fieldSources/basicSource/rotorDiskSource/profileModel/series/seriesProfile.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,8 +28,8 @@ Description Series-up based profile data - drag and lift coefficients computed as sum of cosine series - Cd = sum_i(CdCoeff)*cos(i*AOA) - Cl = sum_i(ClCoeff)*cos(i*AOA) + Cd = sum_i(CdCoeff)*sin(i*AOA) + Cl = sum_i(ClCoeff)*sin(i*AOA) where: AOA = angle of attack [deg] converted to [rad] internally diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C index 4cd984f0ccc59999e923e9214516db3f18cc238b..2c89cc4cae0a0eb08e80800de12ca2a12b3d1a43 100644 --- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C +++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C @@ -70,7 +70,7 @@ bool Foam::pimpleControl::criteriaSatisfied() const label fieldI = applyToField(variableName); if (fieldI != -1) { - const List<lduMatrix::solverPerformance> sp(iter().stream()); + const List<solverPerformance> sp(iter().stream()); const scalar residual = sp.last().initialResidual(); checked = true; diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C index b944f38b8f4a17eff3705a96389462a8f71e0e4d..27557eac3c6b3a736f252992eb6525b86fafc26e 100644 --- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C +++ b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,7 +59,7 @@ bool Foam::simpleControl::criteriaSatisfied() const label fieldI = applyToField(variableName); if (fieldI != -1) { - const List<lduMatrix::solverPerformance> sp(iter().stream()); + const List<solverPerformance> sp(iter().stream()); const scalar residual = sp.first().initialResidual(); checked = true; diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index a69f7f5add6d6ee9c1d3a6c2ffd52f3da59873ad..87a6647ca9cda019b499c4757498da6f2c452c62 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -841,25 +841,23 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const ); volScalarField& H1_ = tH1(); - // Loop over field components - /* - for (direction cmpt=0; cmpt<Type::nComponents; cmpt++) - { - scalarField psiCmpt(psi_.internalField().component(cmpt)); + H1_.internalField() = lduMatrix::H1(); - scalarField boundaryDiagCmpt(psi_.size(), 0.0); - addBoundaryDiag(boundaryDiagCmpt, cmpt); - boundaryDiagCmpt.negate(); - addCmptAvBoundaryDiag(boundaryDiagCmpt); + forAll(psi_.boundaryField(), patchI) + { + const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI]; - H1_.internalField().replace(cmpt, boundaryDiagCmpt); + if (ptf.coupled() && ptf.size()) + { + addToInternalField + ( + lduAddr().patchAddr(patchI), + boundaryCoeffs_[patchI].component(0), + H1_ + ); + } } - H1_.internalField() += lduMatrix::H1(); - */ - - H1_.internalField() = lduMatrix::H1(); - H1_.internalField() /= psi_.mesh().V(); H1_.correctBoundaryConditions(); @@ -1377,7 +1375,7 @@ void Foam::checkMethod template<class Type> -Foam::lduMatrix::solverPerformance Foam::solve +Foam::solverPerformance Foam::solve ( fvMatrix<Type>& fvm, const dictionary& solverControls @@ -1387,13 +1385,13 @@ Foam::lduMatrix::solverPerformance Foam::solve } template<class Type> -Foam::lduMatrix::solverPerformance Foam::solve +Foam::solverPerformance Foam::solve ( const tmp<fvMatrix<Type> >& tfvm, const dictionary& solverControls ) { - lduMatrix::solverPerformance solverPerf = + solverPerformance solverPerf = const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls); tfvm.clear(); @@ -1403,15 +1401,15 @@ Foam::lduMatrix::solverPerformance Foam::solve template<class Type> -Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm) +Foam::solverPerformance Foam::solve(fvMatrix<Type>& fvm) { return fvm.solve(); } template<class Type> -Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm) +Foam::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm) { - lduMatrix::solverPerformance solverPerf = + solverPerformance solverPerf = const_cast<fvMatrix<Type>&>(tfvm()).solve(); tfvm.clear(); diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H index f567254f8655c1b5762fc6450297a93555d5a851..ce214e60c801e9af8346e4c5b7671476bc9045c4 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H @@ -239,11 +239,11 @@ public: //- Solve returning the solution statistics. // Use the given solver controls - lduMatrix::solverPerformance solve(const dictionary&); + solverPerformance solve(const dictionary&); //- Solve returning the solution statistics. // Solver controls read from fvSolution - lduMatrix::solverPerformance solve(); + solverPerformance solve(); }; @@ -387,13 +387,21 @@ public: // Solver controls read from fvSolution autoPtr<fvSolver> solver(); - //- Solve returning the solution statistics. + //- Solve segregated or coupled returning the solution statistics. + // Use the given solver controls + solverPerformance solve(const dictionary&); + + //- Solve segregated returning the solution statistics. + // Use the given solver controls + solverPerformance solveSegregated(const dictionary&); + + //- Solve coupled returning the solution statistics. // Use the given solver controls - lduMatrix::solverPerformance solve(const dictionary&); + solverPerformance solveCoupled(const dictionary&); //- Solve returning the solution statistics. // Solver controls read from fvSolution - lduMatrix::solverPerformance solve(); + solverPerformance solve(); //- Return the matrix residual tmp<Field<Type> > residual() const; @@ -542,14 +550,14 @@ void checkMethod //- Solve returning the solution statistics given convergence tolerance // Use the given solver controls template<class Type> -lduMatrix::solverPerformance solve(fvMatrix<Type>&, const dictionary&); +solverPerformance solve(fvMatrix<Type>&, const dictionary&); //- Solve returning the solution statistics given convergence tolerance, // deleting temporary matrix after solution. // Use the given solver controls template<class Type> -lduMatrix::solverPerformance solve +solverPerformance solve ( const tmp<fvMatrix<Type> >&, const dictionary& @@ -559,14 +567,14 @@ lduMatrix::solverPerformance solve //- Solve returning the solution statistics given convergence tolerance // Solver controls read fvSolution template<class Type> -lduMatrix::solverPerformance solve(fvMatrix<Type>&); +solverPerformance solve(fvMatrix<Type>&); //- Solve returning the solution statistics given convergence tolerance, // deleting temporary matrix after solution. // Solver controls read fvSolution template<class Type> -lduMatrix::solverPerformance solve(const tmp<fvMatrix<Type> >&); +solverPerformance solve(const tmp<fvMatrix<Type> >&); //- Return the correction form of the given matrix diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C index bc3fc66d018ee99dfd90359723478396b76a9514..0cd882fe5463391e6c22fa1715279e7563232e0c 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C @@ -23,6 +23,9 @@ License \*---------------------------------------------------------------------------*/ +#include "LduMatrix.H" +#include "diagTensorField.H" + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class Type> @@ -50,7 +53,7 @@ void Foam::fvMatrix<Type>::setComponentReference template<class Type> -Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve +Foam::solverPerformance Foam::fvMatrix<Type>::solve ( const dictionary& solverControls ) @@ -62,12 +65,51 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve << endl; } + word type(solverControls.lookupOrDefault<word>("type", "segregated")); + + if (type == "segregated") + { + return solveSegregated(solverControls); + } + else if (type == "coupled") + { + return solveCoupled(solverControls); + } + else + { + FatalIOErrorIn + ( + "fvMatrix<Type>::solve(const dictionary& solverControls)", + solverControls + ) << "Unknown type " << type + << "; currently supported solver types are segregated and coupled" + << exit(FatalIOError); + + return solverPerformance(); + } +} + + +template<class Type> +Foam::solverPerformance Foam::fvMatrix<Type>::solveSegregated +( + const dictionary& solverControls +) +{ + if (debug) + { + Info<< "fvMatrix<Type>::solveSegregated" + "(const dictionary& solverControls) : " + "solving fvMatrix<Type>" + << endl; + } + GeometricField<Type, fvPatchField, volMesh>& psi = const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_); - lduMatrix::solverPerformance solverPerfVec + solverPerformance solverPerfVec ( - "fvMatrix<Type>::solve", + "fvMatrix<Type>::solveSegregated", psi.name() ); @@ -134,7 +176,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve cmpt ); - lduMatrix::solverPerformance solverPerf; + solverPerformance solverPerf; // Solver call solverPerf = lduMatrix::solver::New @@ -147,7 +189,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve solverControls )->solve(psiCmpt, sourceCmpt, cmpt); - solverPerf.print(); + solverPerf.print(Info); solverPerfVec = max(solverPerfVec, solverPerf); solverPerfVec.solverName() = solverPerf.solverName(); @@ -164,6 +206,62 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve } +template<class Type> +Foam::solverPerformance Foam::fvMatrix<Type>::solveCoupled +( + const dictionary& solverControls +) +{ + if (debug) + { + Info<< "fvMatrix<Type>::solveCoupled" + "(const dictionary& solverControls) : " + "solving fvMatrix<Type>" + << endl; + } + + GeometricField<Type, fvPatchField, volMesh>& psi = + const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_); + + LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh()); + coupledMatrix.diag() = diag(); + coupledMatrix.upper() = upper(); + coupledMatrix.lower() = lower(); + coupledMatrix.source() = source(); + + addBoundaryDiag(coupledMatrix.diag(), 0); + addBoundarySource(coupledMatrix.source(), false); + + coupledMatrix.interfaces() = psi.boundaryField().interfaces(); + coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0); + coupledMatrix.interfacesLower() = internalCoeffs().component(0); + + autoPtr<typename LduMatrix<Type, scalar, scalar>::solver> + coupledMatrixSolver + ( + LduMatrix<Type, scalar, scalar>::solver::New + ( + psi.name(), + coupledMatrix, + solverControls + ) + ); + + SolverPerformance<Type> solverPerf + ( + coupledMatrixSolver->solve(psi) + ); + + solverPerf.print(Info); + + psi.correctBoundaryConditions(); + + // psi.mesh().setSolverPerformance(psi.name(), solverPerf); + + return solverPerformance(); +} + + template<class Type> Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver> Foam::fvMatrix<Type>::solver() @@ -183,7 +281,7 @@ Foam::fvMatrix<Type>::solver() template<class Type> -Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve() +Foam::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve() { return solve ( @@ -200,7 +298,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve() template<class Type> -Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve() +Foam::solverPerformance Foam::fvMatrix<Type>::solve() { return solve ( diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C index 23a259d574552e48cd73d2214ce8628122ad1182..bd75991ca9cebca2d7031b0f04a5c022055e34ad 100644 --- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C +++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C @@ -93,7 +93,7 @@ Foam::fvMatrix<Foam::scalar>::solver template<> -Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve +Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve ( const dictionary& solverControls ) @@ -111,13 +111,13 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve // assign new solver controls solver_->read(solverControls); - lduMatrix::solverPerformance solverPerf = solver_->solve + solverPerformance solverPerf = solver_->solve ( psi.internalField(), totalSource ); - solverPerf.print(); + solverPerf.print(Info); fvMat_.diag() = saveDiag; @@ -130,14 +130,15 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve template<> -Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve +Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated ( const dictionary& solverControls ) { if (debug) { - Info<< "fvMatrix<scalar>::solve(const dictionary& solverControls) : " + Info<< "fvMatrix<scalar>::solveSegregated" + "(const dictionary& solverControls) : " "solving fvMatrix<scalar>" << endl; } @@ -152,7 +153,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve addBoundarySource(totalSource, false); // Solver call - lduMatrix::solverPerformance solverPerf = lduMatrix::solver::New + solverPerformance solverPerf = lduMatrix::solver::New ( psi.name(), *this, @@ -162,7 +163,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve solverControls )->solve(psi.internalField(), totalSource); - solverPerf.print(); + solverPerf.print(Info); diag() = saveDiag; diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H index 072c9f2cd09e0de014046fb5dfa78b41546bd239..dc0d89db4d5a2bb5da4c87c588cfb6506c53ee5e 100644 --- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H +++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -62,13 +62,13 @@ autoPtr<fvMatrix<scalar>::fvSolver> fvMatrix<scalar>::solver ); template<> -lduMatrix::solverPerformance fvMatrix<scalar>::fvSolver::solve +solverPerformance fvMatrix<scalar>::fvSolver::solve ( const dictionary& ); template<> -lduMatrix::solverPerformance fvMatrix<scalar>::solve +solverPerformance fvMatrix<scalar>::solveSegregated ( const dictionary& ); diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C index d7134f5cd370c780767ae55fd113b24995fef91a..e65329296b73a05ad79ea6fa170164943be6a222 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C @@ -70,8 +70,7 @@ void Foam::KinematicCloud<CloudType>::setModels() SurfaceFilmModel<KinematicCloud<CloudType> >::New ( subModelProperties_, - *this, - g_ + *this ).ptr() ); diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H index ab0e7ed940587014b345ede6e14d34103450bf01..c7db0661878739b505930d6e03b335e57caa4399 100644 --- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H +++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H @@ -241,7 +241,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::massInSystem() const forAllConstIter(typename KinematicCloud<CloudType>, *this, iter) { const parcelType& p = iter(); - sysMass += p.mass()*p.nParticle(); + sysMass += p.nParticle()*p.mass(); } return sysMass; @@ -258,7 +258,7 @@ Foam::KinematicCloud<CloudType>::linearMomentumOfSystem() const { const parcelType& p = iter(); - linearMomentum += p.mass()*p.U(); + linearMomentum += p.nParticle()*p.mass()*p.U(); } return linearMomentum; @@ -275,7 +275,7 @@ Foam::KinematicCloud<CloudType>::linearKineticEnergyOfSystem() const { const parcelType& p = iter(); - linearKineticEnergy += 0.5*p.mass()*(p.U() & p.U()); + linearKineticEnergy += p.nParticle()*0.5*p.mass()*(p.U() & p.U()); } return linearKineticEnergy; @@ -293,7 +293,7 @@ Foam::KinematicCloud<CloudType>::rotationalKineticEnergyOfSystem() const const parcelType& p = iter(); rotationalKineticEnergy += - 0.5*p.momentOfInertia()*(p.omega() & p.omega()); + p.nParticle()*0.5*p.momentOfInertia()*(p.omega() & p.omega()); } return rotationalKineticEnergy; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.H index 100e33e9e240c2deedf109ff713ae31d9ce3546e..05efc427a49db5ae7d652ec7ca3d25d8a7e1b3c6 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -62,7 +62,7 @@ class CellZoneInjection // Private data //- Name of cell zone - const word& cellZoneName_; + const word cellZoneName_; //- Number density const scalar numberDensity_; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.C index 86c721efc1423e3fca536f88ac02b769773bdcba..399730c95088ccfb0ce125cf0d0046291ec101e5 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,8 +32,7 @@ template<class CloudType> Foam::NoSurfaceFilm<CloudType>::NoSurfaceFilm ( const dictionary&, - CloudType& owner, - const dimensionedVector& + CloudType& owner ) : SurfaceFilmModel<CloudType>(owner) diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.H index 672feb967e6455e90964a530f888794cb6f8fe45..b67136aa6af5d77f49ee387f0a225b5c49d9d728 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/NoSurfaceFilm/NoSurfaceFilm.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,7 +68,7 @@ public: // Constructors //- Construct from dictionary - NoSurfaceFilm(const dictionary&, CloudType&, const dimensionedVector&); + NoSurfaceFilm(const dictionary&, CloudType&); //- Construct copy NoSurfaceFilm(const NoSurfaceFilm<CloudType>& dm); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C index 46fb9296dc1624edc38fbe0475181bc85f4d0d66..7eb6e8cfeeff6aa1c369c357a7a998fe2230d997 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,7 +36,7 @@ template<class CloudType> Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel(CloudType& owner) : SubModelBase<CloudType>(owner), - g_(dimensionedVector("zero", dimAcceleration, vector::zero)), + g_(owner.g()), ejectedParcelType_(0), massParcelPatch_(0), diameterParcelPatch_(0), @@ -53,12 +53,11 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel ( const dictionary& dict, CloudType& owner, - const dimensionedVector& g, const word& type ) : SubModelBase<CloudType>(owner, dict, typeName, type), - g_(g), + g_(owner.g()), ejectedParcelType_ ( this->coeffDict().lookupOrDefault("ejectedParcelType", -1) diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H index 967e748f56c8385483aca4dd3c2f4f280fa42da6..b822ee644a6d2c34c4bcdd8458ce67a57170ae2d 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -140,10 +140,9 @@ public: dictionary, ( const dictionary& dict, - CloudType& owner, - const dimensionedVector& g + CloudType& owner ), - (dict, owner, g) + (dict, owner) ); @@ -157,7 +156,6 @@ public: ( const dictionary& dict, CloudType& owner, - const dimensionedVector& g, const word& type ); @@ -182,8 +180,7 @@ public: static autoPtr<SurfaceFilmModel<CloudType> > New ( const dictionary& dict, - CloudType& owner, - const dimensionedVector& g + CloudType& owner ); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelNew.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelNew.C index 49702268a1604016eda953c49992aead8988e792..b2e7d479667a5e53d3aa6d910f5a204bcacc9007 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelNew.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelNew.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,8 +32,7 @@ Foam::autoPtr<Foam::SurfaceFilmModel<CloudType> > Foam::SurfaceFilmModel<CloudType>::New ( const dictionary& dict, - CloudType& owner, - const dimensionedVector& g + CloudType& owner ) { const word modelType(dict.lookup("surfaceFilmModel")); @@ -59,7 +58,7 @@ Foam::SurfaceFilmModel<CloudType>::New << exit(FatalError); } - return autoPtr<SurfaceFilmModel<CloudType> >(cstrIter()(dict, owner, g)); + return autoPtr<SurfaceFilmModel<CloudType> >(cstrIter()(dict, owner)); } diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C index 6a6c60db3984390169d526e356537dd40f29575a..51ab413fee44dc71f32e745e156fd8c8c881bcd4 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -487,11 +487,10 @@ template<class CloudType> Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm ( const dictionary& dict, - CloudType& owner, - const dimensionedVector& g + CloudType& owner ) : - SurfaceFilmModel<CloudType>(dict, owner, g, typeName), + SurfaceFilmModel<CloudType>(dict, owner, typeName), rndGen_(owner.rndGen()), thermo_ ( diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H index bddbfcc66e9f28034472ff045a00e984850fc896..5868916c2adea66de5a7e72dd1b5c00e66095a8e 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -248,12 +248,7 @@ public: // Constructors //- Construct from components - ThermoSurfaceFilm - ( - const dictionary& dict, - CloudType& owner, - const dimensionedVector& g - ); + ThermoSurfaceFilm(const dictionary& dict, CloudType& owner); //- Construct copy ThermoSurfaceFilm(const ThermoSurfaceFilm<CloudType>& sfm); diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C index 58c74f4c4b55ce3f81a6ed3e8407d4bbe4d219e0..5ff0f7bf2e00989d8d81195c9d35017c23c5fac2 100644 --- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C +++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C @@ -379,12 +379,12 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField { const labelList& faceMap = faceProcAddressing_[procI]; - // Addressing into original field - labelList curAddr(faceMap.size()); // Correctly oriented copy of internal field Field<Type> procInternalField(procField.internalField()); + // Addressing into original field + labelList curAddr(procInternalField.size()); - forAll(faceMap, addrI) + forAll(procInternalField, addrI) { curAddr[addrI] = mag(faceMap[addrI])-1; if (faceMap[addrI] < 0) diff --git a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C index c7c10d0ae933be5ba70426030959091eaeb8b421..0d27e370505b0022588b729711b0f1c88062e76f 100644 --- a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C +++ b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C @@ -320,7 +320,6 @@ void LaunderSharmaKE::correct() C1_*G*epsilon_/k_ - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_) - fvm::Sp(C2_*f2()*rho_*epsilon_/k_, epsilon_) - //+ 0.75*1.5*flameKproduction*epsilon_/k_ + E ); @@ -340,7 +339,6 @@ void LaunderSharmaKE::correct() == G - fvm::SuSp(2.0/3.0*rho_*divU, k_) - fvm::Sp(rho_*(epsilon_ + D)/k_, k_) - //+ flameKproduction ); kEqn().relax(); diff --git a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes index 56aebd7a94fc31d1ae3a7d20453629f08aac63c2..de66185079e101550bf082894052c6f679fa3047 100644 --- a/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes +++ b/tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes @@ -42,7 +42,7 @@ laplacianSchemes { default none; laplacian(muEff,U) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(alphaEff,h) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DBEff,B) Gauss linear corrected; diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/0/U b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/0/U index 0d8701320a12fd4bb641f6b61c3f21c254b28bba..84b8f75f3cfe486af95eb409be5724fed1af8f31 100644 --- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/0/U +++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/0/U @@ -48,7 +48,7 @@ boundaryField } outlet { - type fluxCorrectedVelocity; //inletOutlet; + type pressureInletOutletVelocity; value uniform (0 0 0); inletValue uniform (0 0 0); } diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes index 1772c2a3fd2de453b29a4d86681ed374d8c0286d..391c95e6dc3d311d1bd639b18527b57799b74ceb 100644 --- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes +++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes @@ -52,7 +52,7 @@ laplacianSchemes laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DomegaEff,omega) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(alphaEff,h) Gauss linear corrected; } diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution index 859a3d046413714e9d9bfe3cc04be1a0e22d11b4..084c382e84f96dcde779b055cbbcbd6d39117bc7 100644 --- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution @@ -52,6 +52,7 @@ solvers PIMPLE { momentumPredictor yes; + transonic no; nOuterCorrectors 50; nCorrectors 1; nNonOrthogonalCorrectors 0; @@ -66,6 +67,8 @@ PIMPLE tolerance 0.0001; } } + + turbOnFinalIterOnly off; } relaxationFactors diff --git a/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes index da4ffa1eaab98542875e2f7dc266e492b8620f12..0fa6ebaf3327be761be8e19e8622e7f1ffb22a20 100644 --- a/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes +++ b/tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes @@ -52,7 +52,7 @@ laplacianSchemes laplacian(DepsilonEff,epsilon) Gauss linear orthogonal; laplacian(DREff,R) Gauss linear orthogonal; laplacian(DomegaEff,omega) Gauss linear orthogonal; - laplacian((rho*(1|A(U))),p) Gauss linear orthogonal; + laplacian(Dp,p) Gauss linear orthogonal; laplacian(alphaEff,h) Gauss linear orthogonal; } diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/T b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/T new file mode 100644 index 0000000000000000000000000000000000000000..561ca236be3bbed252ef63b29c24063a3bea4c6b --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/T @@ -0,0 +1,54 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 293; + +boundaryField +{ + + front + { + type zeroGradient; + } + back + { + type zeroGradient; + } + wall + { + type zeroGradient; + } + porosityWall + { + type zeroGradient; + } + + inlet + { + type fixedValue; + value $internalField; + } + outlet + { + type inletOutlet; + value $internalField; + inletValue $internalField; + } +} + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/U b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/U new file mode 100644 index 0000000000000000000000000000000000000000..84b8f75f3cfe486af95eb409be5724fed1af8f31 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/U @@ -0,0 +1,57 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + front + { + type fixedValue; + value uniform (0 0 0); + } + back + { + type fixedValue; + value uniform (0 0 0); + } + wall + { + type fixedValue; + value uniform (0 0 0); + } + porosityWall + { + type slip; + value uniform (0 0 0); + } + inlet + { + type flowRateInletVelocity; + flowRate constant 0.1; + value uniform (0 0 0); + } + outlet + { + type pressureInletOutletVelocity; + value uniform (0 0 0); + inletValue uniform (0 0 0); + } +} + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/alphat b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/alphat new file mode 100644 index 0000000000000000000000000000000000000000..f609fa6c14c41bba0ec2726af2cf56a9ca8c949a --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/alphat @@ -0,0 +1,57 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object alphat; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type alphatWallFunction; + value uniform 0; + } + back + { + type alphatWallFunction; + value uniform 0; + } + wall + { + type alphatWallFunction; + value uniform 0; + } + porosityWall + { + type alphatWallFunction; + value uniform 0; + } + inlet + { + type calculated; + value uniform 0; + } + outlet + { + type calculated; + value uniform 0; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/epsilon b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/epsilon new file mode 100644 index 0000000000000000000000000000000000000000..e4dccfe57777d09934d55a61d67f94e2917585f0 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/epsilon @@ -0,0 +1,59 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object epsilon; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -3 0 0 0 0]; + +internalField uniform 200; + +boundaryField +{ + front + { + type compressible::epsilonWallFunction; + value uniform 200; + } + back + { + type compressible::epsilonWallFunction; + value uniform 200; + } + wall + { + type compressible::epsilonWallFunction; + value uniform 200; + } + porosityWall + { + type compressible::epsilonWallFunction; + value uniform 200; + } + inlet + { + type compressible::turbulentMixingLengthDissipationRateInlet; + mixingLength 0.005; + value uniform 200; + } + outlet + { + type inletOutlet; + inletValue uniform 200; + value uniform 200; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/k b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/k new file mode 100644 index 0000000000000000000000000000000000000000..655a91bb45e3aa597a5280c13d1a8262524a0fbd --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/k @@ -0,0 +1,59 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object k; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + front + { + type compressible::kqRWallFunction; + value uniform 1; + } + back + { + type compressible::kqRWallFunction; + value uniform 1; + } + wall + { + type compressible::kqRWallFunction; + value uniform 1; + } + porosityWall + { + type compressible::kqRWallFunction; + value uniform 1; + } + inlet + { + type turbulentIntensityKineticEnergyInlet; + intensity 0.05; + value uniform 1; + } + outlet + { + type inletOutlet; + inletValue uniform 1; + value uniform 1; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/mut b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/mut new file mode 100644 index 0000000000000000000000000000000000000000..0cea2db2d2e1aa076b70254b0e314b9959cc0d78 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/mut @@ -0,0 +1,57 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object mut; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + front + { + type mutkWallFunction; + value uniform 0; + } + back + { + type mutkWallFunction; + value uniform 0; + } + wall + { + type mutkWallFunction; + value uniform 0; + } + porosityWall + { + type mutkWallFunction; + value uniform 0; + } + inlet + { + type calculated; + value uniform 0; + } + outlet + { + type calculated; + value uniform 0; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/0/p b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/p new file mode 100644 index 0000000000000000000000000000000000000000..21db04d6109a0348229d019cb3a13d728044c86d --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/0/p @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1.0e5; + +boundaryField +{ + front + { + type zeroGradient; + } + back + { + type zeroGradient; + } + wall + { + type zeroGradient; + } + porosityWall + { + type zeroGradient; + } + + inlet + { + type zeroGradient; + } + outlet + { + type fixedValue; + value $internalField; + } +} + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/Allrun b/tutorials/compressible/rhoPimplecFoam/angledDuct/Allrun new file mode 100755 index 0000000000000000000000000000000000000000..c9fa98255753a468f44e925754d110ecf3d56332 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/Allrun @@ -0,0 +1,10 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +runApplication blockMesh +runApplication `getApplication` diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/RASProperties b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/RASProperties new file mode 100644 index 0000000000000000000000000000000000000000..a4937b503a46850b2626f0d301e4a07b9f691507 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/RASProperties @@ -0,0 +1,25 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object RASProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +RASModel kEpsilon; + +turbulence on; + +printCoeffs on; + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/blockMeshDict.m4 b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/blockMeshDict.m4 new file mode 100644 index 0000000000000000000000000000000000000000..79da11e10ac8f655aa16fc6287c9496875d8672f --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/blockMeshDict.m4 @@ -0,0 +1,165 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + `format' ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// block definition for a porosity with an angled inlet/outlet +// the porosity is not aligned with the main axes +// +dnl> ----------------------------------------------------------------- +dnl> <STANDARD DEFINTIONS> +dnl> +changecom(//)changequote([,]) dnl> +define(calc, [esyscmd(perl -e 'print ($1)')]) dnl> +define(VCOUNT, 0) dnl> +define(vlabel, [[// ]pt VCOUNT ($1) define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))]) dnl> +dnl> +define(hex2D, hex ($1b $2b $3b $4b $1f $2f $3f $4f)) dnl> +define(quad2D, ($1f $1b $2b $2f)) dnl> +define(frontQuad, ($1f $2f $3f $4f)) dnl> +define(backQuad, ($4b $3b $2b $1b)) dnl> +dnl> +dnl> </STANDARD DEFINTIONS> +dnl> ----------------------------------------------------------------- +dnl> +define(ncells, 20) dnl> +define(ninlet, 15) dnl> +define(nporo, 20) dnl> +define(noutlet, 20) dnl> +dnl> +define(x0,0) dnl> +define(y0,0) dnl> +define(y0,0) dnl> +define(Cos,0.7071067812) dnl> == cos(45) +define(Sin,0.7071067812) dnl> == sin(45) +dnl> +define(width,50) dnl> +define(zBack,calc(-width/2)) dnl> +define(zFront,calc(width/2)) dnl> +define(leninlet,150)dnl> +define(lenporo,100)dnl> +define(lenoutlet,100)dnl> +dnl> +define(xhyp,calc(Sin*width)) dnl> +define(yhyp,calc(Cos*width)) dnl> +define(xinlet,leninlet)dnl> +define(xporo,calc(Cos*lenporo)) dnl> +define(yporo,calc(Sin*lenporo)) dnl> +define(xoutlet,calc(xporo + Cos*lenoutlet)) dnl> +define(youtlet,calc(yporo + Sin*lenoutlet)) dnl> +dnl> + +convertToMeters 0.001; + +vertices +( + // inlet region + ( -xinlet y0 zBack ) vlabel(in1b) + ( -xinlet yhyp zBack ) vlabel(in2b) + ( -xinlet y0 zFront ) vlabel(in1f) + ( -xinlet yhyp zFront ) vlabel(in2f) + + // join inlet->outlet + ( x0 y0 zBack ) vlabel(join1b) + ( -xhyp yhyp zBack ) vlabel(join2b) + ( x0 y0 zFront ) vlabel(join1f) + ( -xhyp yhyp zFront ) vlabel(join2f) + + // porosity ends ->outlet + ( xporo yporo zBack ) vlabel(poro1b) + ( calc(xporo - xhyp) calc(yporo + yhyp) zBack ) vlabel(poro2b) + ( xporo yporo zFront ) vlabel(poro1f) + ( calc(xporo - xhyp) calc(yporo + yhyp) zFront ) vlabel(poro2f) + + // outlet + ( xoutlet youtlet zBack ) vlabel(out1b) + ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zBack ) vlabel(out2b) + ( xoutlet youtlet zFront ) vlabel(out1f) + ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zFront ) vlabel(out2f) +); + +blocks +( + // inlet block + hex2D(in1, join1, join2, in2) + inlet ( ninlet ncells ncells ) simpleGrading (1 1 1) + + // porosity block + hex2D(join1, poro1, poro2, join2) + porosity ( nporo ncells ncells ) simpleGrading (1 1 1) + + // outlet block + hex2D(poro1, out1, out2, poro2) + outlet ( noutlet ncells ncells ) simpleGrading (1 1 1) +); + +edges +( +); + +patches +( + // is there no way of defining all my 'defaultFaces' to be 'wall'? + wall front + ( + // inlet block + frontQuad(in1, join1, join2, in2) + // outlet block + frontQuad(poro1, out1, out2, poro2) + ) + + wall back + ( + // inlet block + backQuad(in1, join1, join2, in2) + // outlet block + backQuad(poro1, out1, out2, poro2) + ) + + wall wall + ( + // inlet block + quad2D(in1, join1) + quad2D(join2, in2) + // outlet block + quad2D(poro1, out1) + quad2D(out2, poro2) + ) + + wall porosityWall + ( + // porosity block + frontQuad(join1, poro1, poro2, join2) + // porosity block + backQuad(join1, poro1, poro2, join2) + // porosity block + quad2D(join1, poro1) + quad2D(poro2, join2) + ) + + patch inlet + ( + quad2D(in2, in1) + ) + + patch outlet + ( + quad2D(out2, out1) + ) +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/boundary b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/boundary new file mode 100644 index 0000000000000000000000000000000000000000..0abd1608aba0dcb6aa66c9488133a3c4b51c7588 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/polyMesh/boundary @@ -0,0 +1,58 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +6 +( + front + { + type wall; + nFaces 700; + startFace 63400; + } + back + { + type wall; + nFaces 700; + startFace 64100; + } + wall + { + type wall; + nFaces 1400; + startFace 64800; + } + porosityWall + { + type wall; + nFaces 1600; + startFace 66200; + } + inlet + { + type patch; + nFaces 400; + startFace 67800; + } + outlet + { + type patch; + nFaces 400; + startFace 68200; + } +) + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones new file mode 100644 index 0000000000000000000000000000000000000000..afeb6461e13249ef04a0d449173de79ecdd2111f --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/porousZones @@ -0,0 +1,36 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object porousZones; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +1 +( + porosity + { + coordinateSystem + { + e1 (0.70710678 0.70710678 0); + e2 (0 0 1); + } + + Darcy + { + d d [0 -2 0 0 0 0 0] (5e7 -1000 -1000); + f f [0 -1 0 0 0 0 0] (0 0 0); + } + } +) + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/thermophysicalProperties b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/thermophysicalProperties new file mode 100644 index 0000000000000000000000000000000000000000..d6d597d433f7d69e5e2329e6c716ed66b8a99675 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/thermophysicalProperties @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object thermophysicalProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>; + +mixture +{ + specie + { + nMoles 1; + molWeight 28.9; + } + thermodynamics + { + Cp 1007; + Hf 0; + } + transport + { + As 1.4792e-06; + Ts 116; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/turbulenceProperties b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/turbulenceProperties new file mode 100644 index 0000000000000000000000000000000000000000..3721a46a2ead37eb2bf10434bcde59afa9fe9bf6 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/constant/turbulenceProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType RASModel; + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/system/controlDict b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/controlDict new file mode 100644 index 0000000000000000000000000000000000000000..bbcbeb2543fde4ea39ad5df1cb6f14005cc74202 --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/controlDict @@ -0,0 +1,55 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application rhoPimplecFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 10; + +deltaT 1; + +writeControl adjustableRunTime; + +writeInterval 10; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + +adjustTimeStep no; + +maxCo 10; + +maxDeltaT 1; + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSchemes new file mode 100644 index 0000000000000000000000000000000000000000..391c95e6dc3d311d1bd639b18527b57799b74ceb --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSchemes @@ -0,0 +1,76 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; + grad(p) Gauss linear; +} + +divSchemes +{ + default none; + div(phi,U) Gauss upwind; + div(phid,p) Gauss upwind; + div(phi,K) Gauss linear; + div(phi,h) Gauss upwind; + div(phi,k) Gauss upwind; + div(phi,epsilon) Gauss upwind; + div(phi,R) Gauss upwind; + div(phi,omega) Gauss upwind; + div((rho*R)) Gauss linear; + div(R) Gauss linear; + div(U) Gauss linear; + div((muEff*dev2(T(grad(U))))) Gauss linear; +} + +laplacianSchemes +{ + default none; + laplacian(muEff,U) Gauss linear corrected; + laplacian(mut,U) Gauss linear corrected; + laplacian(DkEff,k) Gauss linear corrected; + laplacian(DepsilonEff,epsilon) Gauss linear corrected; + laplacian(DREff,R) Gauss linear corrected; + laplacian(DomegaEff,omega) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; + laplacian(alphaEff,h) Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSolution new file mode 100644 index 0000000000000000000000000000000000000000..4a126ed8e44ac2adea221069d86969dc2c51456a --- /dev/null +++ b/tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSolution @@ -0,0 +1,89 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0.01; + } + + pFinal + { + $p; + tolerance 1e-06; + relTol 0; + } + + "(rho|U|h|k|epsilon|omega)" + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0.1; + } + + "(rho|U|h|k|epsilon|omega)Final" + { + $U; + tolerance 1e-05; + relTol 0; + } + +} + +PIMPLE +{ + momentumPredictor yes; + transonic no; + nOuterCorrectors 50; + nCorrectors 1; + nNonOrthogonalCorrectors 0; + rhoMin rhoMin [ 1 -3 0 0 0 ] 0.1; + rhoMax rhoMax [ 1 -3 0 0 0 ] 3.0; + + residualControl + { + "(U|k|epsilon)" + { + relTol 0; + tolerance 0.0001; + } + } + + turbOnFinalIterOnly off; +} + +relaxationFactors +{ + fields + { + "p.*" 1; + "rho.*" 1; + } + equations + { + "(U|h|k|epsilon|omega).*" 0.85; + "p.*" 1; + } +} + + +// ************************************************************************* // diff --git a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes index 00981f61fcb6a70666b768c6177a089b00cfc7c5..62fadfb8ca3ee3c2eca193b571f5637905547972 100644 --- a/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes +++ b/tutorials/compressible/rhoPorousMRFLTSPimpleFoam/angledDuct/system/fvSchemes @@ -52,7 +52,7 @@ laplacianSchemes laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DomegaEff,omega) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(alphaEff,h) Gauss linear corrected; } diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes index dbca5c1607ad223ed5c440bd179c7708125c49b0..d49880685b11175bebc06136d7a7de0372e73f8c 100644 --- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes +++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes @@ -40,7 +40,7 @@ laplacianSchemes { default none; laplacian(muEff,U) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(alphaEff,h) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes index d237bd73617ff269aa6435fc68b5269a8bca0fec..711e4e3233a2e700ada115f8ffeb80db86b6950d 100644 --- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes +++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicit/system/fvSchemes @@ -39,7 +39,6 @@ laplacianSchemes { laplacian(muEff,U) Gauss linear corrected; laplacian(alphaEff,h) Gauss linear corrected; - laplacian((rho|A(U)),p) Gauss linear corrected; laplacian((rho*rAU),p) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary index 0abd1608aba0dcb6aa66c9488133a3c4b51c7588..a664893343e7216491ca9913ee0171ee31980aba 100644 --- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary +++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary @@ -8,7 +8,7 @@ FoamFile { version 2.0; - format ascii; + format binary; class polyBoundaryMesh; location "constant/polyMesh"; object boundary; diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/U b/tutorials/compressible/rhoSimplecFoam/squareBend/0/U index 70988607c87df6f0cbf431c7736807fcbb9568cd..d5bd0a74cf75eac8797c6d9337e2ea656734ba4d 100644 --- a/tutorials/compressible/rhoSimplecFoam/squareBend/0/U +++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/U @@ -28,7 +28,7 @@ boundaryField inlet { type flowRateInletVelocity; - flowRate constant 0.5; //0.75; + flowRate constant 0.5; value uniform (0 0 0); } outlet diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/decomposeParDict b/tutorials/compressible/rhoSimplecFoam/squareBend/system/decomposeParDict new file mode 100644 index 0000000000000000000000000000000000000000..9e844a62ba7901b16f2d0c8c470ae7b5acb2504a --- /dev/null +++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/decomposeParDict @@ -0,0 +1,45 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 8; + +method hierarchical; + +simpleCoeffs +{ + n (8 1 1); + delta 0.001; +} + +hierarchicalCoeffs +{ + n (4 2 1); + delta 0.001; + order xyz; +} + +manualCoeffs +{ + dataFile ""; +} + +distributed no; + +roots ( ); + + +// ************************************************************************* // diff --git a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes index b4cf8d2e38d798122b9909b0a8e91b78f1c87e23..51799627c881c59a26a6af79531717e41daf7497 100644 --- a/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes +++ b/tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes @@ -40,7 +40,7 @@ divSchemes laplacianSchemes { default none; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(muEff,U) Gauss linear corrected; laplacian(alphaEff,e) Gauss linear corrected; } diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes index 616ee4aa1bb4895620df483460be5b4f0a89e45b..52cd9397231b9b09a6a718aa674299f30e7e9202 100644 --- a/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes +++ b/tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes @@ -40,7 +40,7 @@ divSchemes laplacianSchemes { default none; - laplacian((rho*(1|A(U))),p) Gauss linear orthogonal; + laplacian(Dp,p) Gauss linear orthogonal; laplacian(muEff,U) Gauss linear orthogonal; laplacian(alphaEff,e) Gauss linear orthogonal; } diff --git a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes index b12dc45dcf6d0ba6aaad7e83fbbbdcac2bffd7c9..4f1434c4797986ad7e8fcf905340a7c01a2423c7 100644 --- a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes +++ b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes @@ -46,7 +46,7 @@ laplacianSchemes laplacian(DkEff,k) Gauss linear limited 0.5; laplacian(DREff,R) Gauss linear limited 0.5; laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5; - laplacian((rho*(1|A(U))),p) Gauss linear limited 0.5; + laplacian(Dp,p) Gauss linear limited 0.5; laplacian(alphaEff,e) Gauss linear limited 0.5; } diff --git a/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes b/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes index e7b5a41a604f719f59b61610b77d2b2de68a5e4c..6aa61d831af8b7e2f1f8c47350bfba0ff2b592f3 100644 --- a/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes +++ b/tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes @@ -46,7 +46,7 @@ laplacianSchemes laplacian(DkEff,k) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; laplacian(alphaEff,e) Gauss linear corrected; } diff --git a/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSchemes b/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSchemes index d7aad508a2c0d987406d7ab0b34a888da769ac57..d1a9bb7aa1310b3bee441fbc1800724e6f19a0a3 100644 --- a/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSchemes +++ b/tutorials/compressible/sonicLiquidFoam/decompressionTank/system/fvSchemes @@ -37,7 +37,7 @@ laplacianSchemes { default none; laplacian(mu,U) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian(Dp,p) Gauss linear corrected; } interpolationSchemes diff --git a/wmake/rules/linux64Gcc47/c b/wmake/rules/linux64Gcc47/c new file mode 100644 index 0000000000000000000000000000000000000000..f4114be3143d1210ffea500a2b361008910abed0 --- /dev/null +++ b/wmake/rules/linux64Gcc47/c @@ -0,0 +1,16 @@ +.SUFFIXES: .c .h + +cWARN = -Wall + +cc = gcc -m64 + +include $(RULES)/c$(WM_COMPILE_OPTION) + +cFLAGS = $(GFLAGS) $(cWARN) $(cOPT) $(cDBUG) $(LIB_HEADER_DIRS) -fPIC + +ctoo = $(WM_SCHEDULER) $(cc) $(cFLAGS) -c $$SOURCE -o $@ + +LINK_LIBS = $(cDBUG) + +LINKLIBSO = $(cc) -shared +LINKEXE = $(cc) -Xlinker --add-needed -Xlinker -z -Xlinker nodefs diff --git a/wmake/rules/linux64Gcc47/c++ b/wmake/rules/linux64Gcc47/c++ new file mode 100644 index 0000000000000000000000000000000000000000..98b25ed1fea2a1baa0ad749c09a76bd877ea4a6d --- /dev/null +++ b/wmake/rules/linux64Gcc47/c++ @@ -0,0 +1,21 @@ +.SUFFIXES: .C .cxx .cc .cpp + +c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast + +CC = g++ -m64 + +include $(RULES)/c++$(WM_COMPILE_OPTION) + +ptFLAGS = -DNoRepository -ftemplate-depth-100 + +c++FLAGS = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC + +Ctoo = $(WM_SCHEDULER) $(CC) $(c++FLAGS) -c $$SOURCE -o $@ +cxxtoo = $(Ctoo) +cctoo = $(Ctoo) +cpptoo = $(Ctoo) + +LINK_LIBS = $(c++DBUG) + +LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed +LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed diff --git a/wmake/rules/linux64Gcc47/c++Debug b/wmake/rules/linux64Gcc47/c++Debug new file mode 100644 index 0000000000000000000000000000000000000000..19bdb9c3346fc7a69380dfedd6e7911fe220a965 --- /dev/null +++ b/wmake/rules/linux64Gcc47/c++Debug @@ -0,0 +1,2 @@ +c++DBUG = -ggdb3 -DFULLDEBUG +c++OPT = -O0 -fdefault-inline diff --git a/wmake/rules/linux64Gcc47/c++Opt b/wmake/rules/linux64Gcc47/c++Opt new file mode 100644 index 0000000000000000000000000000000000000000..2aedabd6280a3476bc58db13139a0a3aa579502b --- /dev/null +++ b/wmake/rules/linux64Gcc47/c++Opt @@ -0,0 +1,2 @@ +c++DBUG = +c++OPT = -O3 diff --git a/wmake/rules/linux64Gcc47/c++Prof b/wmake/rules/linux64Gcc47/c++Prof new file mode 100644 index 0000000000000000000000000000000000000000..3bda4dad55e898a8198f6e8bfe21e8d829d7230a --- /dev/null +++ b/wmake/rules/linux64Gcc47/c++Prof @@ -0,0 +1,2 @@ +c++DBUG = -pg +c++OPT = -O2 diff --git a/wmake/rules/linux64Gcc47/cDebug b/wmake/rules/linux64Gcc47/cDebug new file mode 100644 index 0000000000000000000000000000000000000000..72b638f458220e329d52b59e3566a3c807101f9d --- /dev/null +++ b/wmake/rules/linux64Gcc47/cDebug @@ -0,0 +1,2 @@ +cDBUG = -ggdb -DFULLDEBUG +cOPT = -O1 -fdefault-inline -finline-functions diff --git a/wmake/rules/linux64Gcc47/cOpt b/wmake/rules/linux64Gcc47/cOpt new file mode 100644 index 0000000000000000000000000000000000000000..17318709f1fa39e6bf89cbe87778bc6fa459de17 --- /dev/null +++ b/wmake/rules/linux64Gcc47/cOpt @@ -0,0 +1,2 @@ +cDBUG = +cOPT = -O3 diff --git a/wmake/rules/linux64Gcc47/cProf b/wmake/rules/linux64Gcc47/cProf new file mode 100644 index 0000000000000000000000000000000000000000..ca3ac9bf5f0cd61fe99e0f05fa1bd4bdf9fa6cf7 --- /dev/null +++ b/wmake/rules/linux64Gcc47/cProf @@ -0,0 +1,2 @@ +cDBUG = -pg +cOPT = -O2 diff --git a/wmake/rules/linux64Gcc47/general b/wmake/rules/linux64Gcc47/general new file mode 100644 index 0000000000000000000000000000000000000000..4a42b11b1ee6aebe596e3cd2e2633d9c70cb5e9a --- /dev/null +++ b/wmake/rules/linux64Gcc47/general @@ -0,0 +1,8 @@ +CPP = cpp -traditional-cpp + +PROJECT_LIBS = -l$(WM_PROJECT) -ldl + +include $(GENERAL_RULES)/standard + +include $(RULES)/c +include $(RULES)/c++ diff --git a/wmake/rules/linux64Gcc47/mplibHPMPI b/wmake/rules/linux64Gcc47/mplibHPMPI new file mode 100644 index 0000000000000000000000000000000000000000..574492a236a32f7d87d00bf0e3507a5ac8e54f55 --- /dev/null +++ b/wmake/rules/linux64Gcc47/mplibHPMPI @@ -0,0 +1,3 @@ +PFLAGS = +PINC = -I$(MPI_ARCH_PATH)/include -D_MPICC_H +PLIBS = -L$(MPI_ARCH_PATH)/lib/linux_amd64 -lmpi diff --git a/wmake/rules/linux64Gcc47/mplibINTELMPI b/wmake/rules/linux64Gcc47/mplibINTELMPI new file mode 100644 index 0000000000000000000000000000000000000000..cf80ec2eaf68d1c2f6adf208964b6490c4c8fd36 --- /dev/null +++ b/wmake/rules/linux64Gcc47/mplibINTELMPI @@ -0,0 +1,3 @@ +PFLAGS = -DMPICH_SKIP_MPICXX +PINC = -I$(MPI_ARCH_PATH)/include64 +PLIBS = -L$(MPI_ARCH_PATH)/lib64 -lmpi diff --git a/wmake/rules/linuxGcc47/c b/wmake/rules/linuxGcc47/c new file mode 100644 index 0000000000000000000000000000000000000000..d914fcd37d084b82a1833722d6ad7a0db3dd1c93 --- /dev/null +++ b/wmake/rules/linuxGcc47/c @@ -0,0 +1,16 @@ +.SUFFIXES: .c .h + +cWARN = -Wall + +cc = gcc -m32 + +include $(RULES)/c$(WM_COMPILE_OPTION) + +cFLAGS = $(GFLAGS) $(cWARN) $(cOPT) $(cDBUG) $(LIB_HEADER_DIRS) -fPIC + +ctoo = $(WM_SCHEDULER) $(cc) $(cFLAGS) -c $$SOURCE -o $@ + +LINK_LIBS = $(cDBUG) + +LINKLIBSO = $(cc) -shared +LINKEXE = $(cc) -Xlinker --add-needed -Xlinker -z -Xlinker nodefs diff --git a/wmake/rules/linuxGcc47/c++ b/wmake/rules/linuxGcc47/c++ new file mode 100644 index 0000000000000000000000000000000000000000..357f4106e10d7d1108a6713802e6f0d01cd8be9a --- /dev/null +++ b/wmake/rules/linuxGcc47/c++ @@ -0,0 +1,21 @@ +.SUFFIXES: .C .cxx .cc .cpp + +c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast + +CC = g++ -m32 + +include $(RULES)/c++$(WM_COMPILE_OPTION) + +ptFLAGS = -DNoRepository -ftemplate-depth-100 + +c++FLAGS = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC + +Ctoo = $(WM_SCHEDULER) $(CC) $(c++FLAGS) -c $$SOURCE -o $@ +cxxtoo = $(Ctoo) +cctoo = $(Ctoo) +cpptoo = $(Ctoo) + +LINK_LIBS = $(c++DBUG) + +LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed +LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed diff --git a/wmake/rules/linuxGcc47/c++Debug b/wmake/rules/linuxGcc47/c++Debug new file mode 100644 index 0000000000000000000000000000000000000000..19bdb9c3346fc7a69380dfedd6e7911fe220a965 --- /dev/null +++ b/wmake/rules/linuxGcc47/c++Debug @@ -0,0 +1,2 @@ +c++DBUG = -ggdb3 -DFULLDEBUG +c++OPT = -O0 -fdefault-inline diff --git a/wmake/rules/linuxGcc47/c++Opt b/wmake/rules/linuxGcc47/c++Opt new file mode 100644 index 0000000000000000000000000000000000000000..2aedabd6280a3476bc58db13139a0a3aa579502b --- /dev/null +++ b/wmake/rules/linuxGcc47/c++Opt @@ -0,0 +1,2 @@ +c++DBUG = +c++OPT = -O3 diff --git a/wmake/rules/linuxGcc47/c++Prof b/wmake/rules/linuxGcc47/c++Prof new file mode 100644 index 0000000000000000000000000000000000000000..3bda4dad55e898a8198f6e8bfe21e8d829d7230a --- /dev/null +++ b/wmake/rules/linuxGcc47/c++Prof @@ -0,0 +1,2 @@ +c++DBUG = -pg +c++OPT = -O2 diff --git a/wmake/rules/linuxGcc47/cDebug b/wmake/rules/linuxGcc47/cDebug new file mode 100644 index 0000000000000000000000000000000000000000..72b638f458220e329d52b59e3566a3c807101f9d --- /dev/null +++ b/wmake/rules/linuxGcc47/cDebug @@ -0,0 +1,2 @@ +cDBUG = -ggdb -DFULLDEBUG +cOPT = -O1 -fdefault-inline -finline-functions diff --git a/wmake/rules/linuxGcc47/cOpt b/wmake/rules/linuxGcc47/cOpt new file mode 100644 index 0000000000000000000000000000000000000000..17318709f1fa39e6bf89cbe87778bc6fa459de17 --- /dev/null +++ b/wmake/rules/linuxGcc47/cOpt @@ -0,0 +1,2 @@ +cDBUG = +cOPT = -O3 diff --git a/wmake/rules/linuxGcc47/cProf b/wmake/rules/linuxGcc47/cProf new file mode 100644 index 0000000000000000000000000000000000000000..ca3ac9bf5f0cd61fe99e0f05fa1bd4bdf9fa6cf7 --- /dev/null +++ b/wmake/rules/linuxGcc47/cProf @@ -0,0 +1,2 @@ +cDBUG = -pg +cOPT = -O2 diff --git a/wmake/rules/linuxGcc47/general b/wmake/rules/linuxGcc47/general new file mode 100644 index 0000000000000000000000000000000000000000..4b051e6b9840df29cac2e8ced14c7d18b0b337d5 --- /dev/null +++ b/wmake/rules/linuxGcc47/general @@ -0,0 +1,9 @@ +CPP = cpp -traditional-cpp +LD = ld -melf_i386 + +PROJECT_LIBS = -l$(WM_PROJECT) -ldl + +include $(GENERAL_RULES)/standard + +include $(RULES)/c +include $(RULES)/c++ diff --git a/wmake/rules/linuxGcc47/mplibHPMPI b/wmake/rules/linuxGcc47/mplibHPMPI new file mode 100644 index 0000000000000000000000000000000000000000..8aff40632bd23af9607d63c4eb675a8de0cd287c --- /dev/null +++ b/wmake/rules/linuxGcc47/mplibHPMPI @@ -0,0 +1,3 @@ +PFLAGS = +PINC = -I$(MPI_ARCH_PATH)/include -D_MPICC_H +PLIBS = -L$(MPI_ARCH_PATH)/lib/linux_ia32 -lmpi