diff --git a/README.txt b/README.txt new file mode 100644 index 0000000000..5e0378f4a9 --- /dev/null +++ b/README.txt @@ -0,0 +1,10 @@ +2013-01 + +- to UPstream added allocation of communicators +- added communicators to lduMesh,lduInterface and (indirectly) + to lduInterfaceField +- gamg agglomeration allocates new communicator for every level +- in all linear solvers/smoothers/preconditioners make they use + communicator +- added lots of warnings if using unexpected communicator (UPstream::warnComm) +- did LUScalarMatrix for 'directSolveCoarsest' diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C index d636481b31..71662dec34 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C +++ b/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C index 24752a166a..c0b4b73afa 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoPimpleFoam/Allwmake b/applications/solvers/compressible/rhoPimpleFoam/Allwmake index ac06b7350a..1b272c17f0 100755 --- a/applications/solvers/compressible/rhoPimpleFoam/Allwmake +++ b/applications/solvers/compressible/rhoPimpleFoam/Allwmake @@ -5,5 +5,6 @@ set -x wmake wmake rhoPimplecFoam wmake rhoLTSPimpleFoam +wmake rhoPimpleDyMFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/files b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/files new file mode 100644 index 0000000000..034b5c2b1b --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/files @@ -0,0 +1,3 @@ +rhoPimpleDyMFoam.C + +EXE = $(FOAM_APPBIN)/rhoPimpleDyMFoam diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/options new file mode 100644 index 0000000000..2093b620b5 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/options @@ -0,0 +1,27 @@ +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 \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +EXE_LIBS = \ + -lfluidThermophysicalModels \ + -lspecie \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lfvOptions \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H new file mode 100644 index 0000000000..24316a1685 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H @@ -0,0 +1,59 @@ +{ + if (mesh.changing()) + { + forAll(U.boundaryField(), patchi) + { + if (U.boundaryField()[patchi].fixesValue()) + { + U.boundaryField()[patchi].initEvaluate(); + } + } + + forAll(U.boundaryField(), patchi) + { + if (U.boundaryField()[patchi].fixesValue()) + { + U.boundaryField()[patchi].evaluate(); + + phi.boundaryField()[patchi] = + rho.boundaryField()[patchi] + *( + U.boundaryField()[patchi] + & mesh.Sf().boundaryField()[patchi] + ); + } + } + } + + volScalarField pcorr + ( + IOobject + ( + "pcorr", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("pcorr", p.dimensions(), 0.0), + pcorrTypes + ); + + dimensionedScalar Dp("Dp", dimTime, 1.0); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pcorrEqn + ( + fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divrhoU + ); + + pcorrEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi -= pcorrEqn.flux(); + } + } +} diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H new file mode 100644 index 0000000000..a602fd4843 --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H @@ -0,0 +1,13 @@ + wordList pcorrTypes + ( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + for (label i=0; i("correctPhi", true); + + bool checkMeshCourantNo = + pimple.dict().lookupOrDefault("checkMeshCourantNo", false); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C new file mode 100644 index 0000000000..fb74f2b39a --- /dev/null +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + rhoPimpleFoam + +Description + Transient solver for laminar or turbulent flow of compressible fluids + for HVAC and similar applications. + + Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and + pseudo-transient simulations. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "psiThermo.H" +#include "turbulenceModel.H" +#include "bound.H" +#include "pimpleControl.H" +#include "fvIOoptionList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "initContinuityErrs.H" + + pimpleControl pimple(mesh); + + #include "readControls.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "createPcorrTypes.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // Create old-time absolute flux for ddtPhiCorr + surfaceScalarField phiAbs("phiAbs", phi); + + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readControls.H" + #include "compressibleCourantNo.H" + + // Make the fluxes absolute before mesh-motion + fvc::makeAbsolute(phi, rho, U); + + // Update absolute flux for ddtPhiCorr + phiAbs = phi; + + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + { + // Store divrhoU from the previous time-step/mesh for the correctPhi + volScalarField divrhoU(fvc::div(phi)); + + // Do any mesh changes + mesh.update(); + + if (mesh.changing() && correctPhi) + { + #include "correctPhi.H" + } + } + + // Make the fluxes relative to the mesh-motion + fvc::makeRelative(phi, rho, U); + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + + if (pimple.nCorrPIMPLE() <= 1) + { + #include "rhoEqn.H" + Info<< "rhoEqn max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + } + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "EEqn.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/incompressible/pimpleFoam/pimpleDyMFoam/Make/options b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/Make/options index 5cc44fea11..263fd7375f 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/Make/options +++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/Make/options @@ -1,25 +1,24 @@ EXE_INC = \ -I.. \ - -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/fvOptions/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude - + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ EXE_LIBS = \ - -ldynamicFvMesh \ - -ltopoChangerFvMesh \ - -ldynamicMesh \ - -lmeshTools \ -lincompressibleTransportModels \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleLESModels \ -lfiniteVolume \ -lfvOptions \ - -lsampling + -lsampling \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -ldynamicMesh \ + -lmeshTools diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H deleted file mode 100644 index eaaa3e60c9..0000000000 --- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H +++ /dev/null @@ -1,23 +0,0 @@ -// Solve the Momentum equation - -tmp UEqn -( - fvm::ddt(U) - + fvm::div(phi, U) - + turbulence->divDevReff(U) - == - fvOptions(U) -); - -UEqn().relax(); - -fvOptions.constrain(UEqn()); - -rAU = 1.0/UEqn().A(); - -if (pimple.momentumPredictor()) -{ - solve(UEqn() == -fvc::grad(p)); - - fvOptions.correct(U); -} diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H index 9bed803d1e..e1e897c18b 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H +++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H @@ -22,20 +22,6 @@ } } - wordList pcorrTypes - ( - p.boundaryField().size(), - zeroGradientFvPatchScalarField::typeName - ); - - forAll(p.boundaryField(), patchI) - { - if (p.boundaryField()[patchI].fixesValue()) - { - pcorrTypes[patchI] = fixedValueFvPatchScalarField::typeName; - } - } - volScalarField pcorr ( IOobject @@ -51,11 +37,13 @@ pcorrTypes ); + dimensionedScalar Dp("Dp", dimTime, 1.0); + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pcorrEqn ( - fvm::laplacian(rAU, pcorr) == fvc::div(phi) + fvm::laplacian(Dp, pcorr) == fvc::div(phi) ); pcorrEqn.setReference(pRefCell, pRefValue); @@ -68,6 +56,4 @@ } } -phi.oldTime() = phi; - #include "continuityErrs.H" diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H index 16b3bd977d..082ec86ea8 100644 --- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H +++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H @@ -40,19 +40,3 @@ ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); - - Info<< "Reading field rAU if present\n" << endl; - volScalarField rAU - ( - IOobject - ( - "rAU", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - mesh, - runTime.deltaT(), - zeroGradientFvPatchScalarField::typeName - ); diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H new file mode 100644 index 0000000000..a602fd4843 --- /dev/null +++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H @@ -0,0 +1,13 @@ + wordList pcorrTypes + ( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + for (label i=0; i 1 && pimple.nCorrPIMPLE() != 1) - { - FatalErrorIn(args.executable()) - << "Sub-cycling alpha is only allowed for PISO, " - "i.e. when the number of outer-correctors = 1" - << exit(FatalError); - } - bool correctPhi = pimple.dict().lookupOrDefault("correctPhi", true); diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index a42d911980..ed5faa7c86 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C @@ -56,7 +56,7 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); - #include "readControls.H" + #include "readTimeControls.H" #include "initContinuityErrs.H" #include "createFields.H" #include "CourantNo.H" @@ -68,7 +68,7 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readControls.H" + #include "readTimeControls.H" #include "CourantNo.H" #include "setDeltaT.H" diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 73babb08f0..1319378821 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -85,8 +85,8 @@ if (pimple.finalNonOrthogonalIter()) { - //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); - //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; + p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); + p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh; dgdt = ( @@ -102,7 +102,7 @@ } } - p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); + // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin); // Update densities from change in p_rgh rho1 += psi1*(p_rgh - p_rgh_0); diff --git a/applications/solvers/multiphase/compressibleInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H deleted file mode 100644 index 87a055d641..0000000000 --- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H +++ /dev/null @@ -1,13 +0,0 @@ - #include "readTimeControls.H" - - label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); - - label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); - - if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1) - { - FatalErrorIn(args.executable()) - << "Sub-cycling alpha is only allowed for PISO operation, " - "i.e. when the number of outer-correctors = 1" - << exit(FatalError); - } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H index 1f12fc9a12..1810bdd4d9 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H @@ -28,8 +28,8 @@ + ( he1.name() == thermo1.phasePropertyName("e") - ? fvc::div(alphaPhi1, p) - : -dalpha1pdt + ? fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p) + : -alpha1*dpdt )/rho1 - fvm::laplacian(k1, he1) @@ -50,8 +50,8 @@ + ( he2.name() == thermo2.phasePropertyName("e") - ? fvc::div(alphaPhi2, p) - : -dalpha2pdt + ? fvc::ddt(alpha2)*p + fvc::div(alphaPhi2, p) + : -alpha2*dpdt )/rho2 - fvm::laplacian(k2, he2) diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H index cc5b3ebe21..5b030b95ea 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H @@ -275,8 +275,8 @@ ); - Info<< "Creating field dalpha1pdt\n" << endl; - volScalarField dalpha1pdt + Info<< "Creating field dpdt\n" << endl; + volScalarField dpdt ( IOobject ( @@ -285,21 +285,9 @@ mesh ), mesh, - dimensionedScalar("dalpha1pdt", p.dimensions()/dimTime, 0) + dimensionedScalar("dpdt", p.dimensions()/dimTime, 0) ); - Info<< "Creating field dalpha2pdt\n" << endl; - volScalarField dalpha2pdt - ( - IOobject - ( - "dpdt", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar("dalpha2pdt", p.dimensions()/dimTime, 0) - ); Info<< "Creating field kinetic energy K\n" << endl; volScalarField K1("K" + phase1Name, 0.5*magSqr(U1)); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index 02e49b5661..756a231fc2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -208,13 +208,8 @@ K1 = 0.5*magSqr(U1); K2 = 0.5*magSqr(U2); - if (thermo1.dpdt()) + if (thermo1.dpdt() || thermo2.dpdt()) { - dalpha1pdt = fvc::ddt(alpha1, p); - } - - if (thermo2.dpdt()) - { - dalpha2pdt = fvc::ddt(alpha2, p); + dpdt = fvc::ddt(p); } } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H index cb8b4efc3b..9913595cf2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H @@ -1,4 +1,2 @@ #include "readTimeControls.H" - - int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr"))); - int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles"))); + #include "alphaControls.H" diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H index 57c78027a4..428876fd22 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H @@ -1,5 +1,4 @@ -label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); -label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); +#include "alphaControls.H" if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H index 8525a8bf67..f419153506 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -130,10 +130,7 @@ << ", " << gMax(1/rDeltaT.internalField()) << endl; } - label nAlphaSubCycles - ( - readLabel(pimpleDict.lookup("nAlphaSubCycles")) - ); + #include "alphaControls.H" rSubDeltaT = rDeltaT*nAlphaSubCycles; } diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H index 57c78027a4..428876fd22 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H @@ -1,5 +1,4 @@ -label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); -label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); +#include "alphaControls.H" if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H index c4cdbc044b..8d3e60da40 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H @@ -50,11 +50,11 @@ { phi -= pcorrEqn.flux(); phiAbs = phi; + phiAbs.oldTime() = phi; fvc::makeRelative(phi, U); + phi.oldTime() = phi; } } } -phi.oldTime() = phi; - #include "continuityErrs.H" diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index fa16f17141..3ccfa73fee 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -51,10 +51,11 @@ int main(int argc, char *argv[]) #include "createDynamicFvMesh.H" #include "initContinuityErrs.H" #include "createFields.H" - #include "readTimeControls.H" pimpleControl pimple(mesh); + #include "readTimeControls.H" + surfaceScalarField phiAbs("phiAbs", phi); fvc::makeAbsolute(phiAbs, U); diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C index 51f082b616..dccf3f8b28 100644 --- a/applications/solvers/multiphase/interFoam/interFoam.C +++ b/applications/solvers/multiphase/interFoam/interFoam.C @@ -63,7 +63,6 @@ int main(int argc, char *argv[]) #include "CourantNo.H" #include "setInitialDeltaT.H" - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H index 97a09ce017..9e3f4f6bd4 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H @@ -1,6 +1,4 @@ -const dictionary& pimpleDict = pimple.dict(); -label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr"))); -label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); +#include "alphaControls.H" if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C index e16f83bd2e..f812b472e0 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -174,8 +174,10 @@ Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties ( readScalar ( - mixture.U().mesh().solutionDict().subDict("PIMPLE"). - lookup("cAlpha") + mixture.U().mesh().solverDict + ( + mixture_.alpha1().name() + ).lookup("cAlpha") ) ), sigma12_(mixture.lookup("sigma12")), diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options index 58b340781f..d4ead14721 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options @@ -5,7 +5,10 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \ - -I$(LIB_SRC)/finiteVolume/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude\ + -I$(LIB_SRC)/sampling/lnInclude EXE_LIBS = \ -ltwoPhaseMixture \ @@ -15,4 +18,7 @@ EXE_LIBS = \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleLESModels \ - -lfiniteVolume + -lfiniteVolume \ + -lmeshTools \ + -lfvOptions \ + -lsampling diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index 9ce20db849..9e5461614b 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -4,9 +4,41 @@ surfaceScalarField phir("phir", phic*interface.nHatf()); - for (int gCorr=0; gCorr > vDotAlphal = + twoPhaseProperties->vDotAlphal(); + const volScalarField& vDotcAlphal = vDotAlphal[0](); + const volScalarField& vDotvAlphal = vDotAlphal[1](); + const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal); + + tmp tphiAlpha; + + if (MULESCorr) { - surfaceScalarField phiAlpha + fvScalarMatrix alpha1Eqn + ( + fvm::ddt(alpha1) + + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1) + == + fvm::Sp(vDotvmcAlphal, alpha1) + + vDotcAlphal + ); + + alpha1Eqn.solve(); + + Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + + tphiAlpha = alpha1Eqn.flux(); + } + + volScalarField alpha10("alpha10", alpha1); + + for (int aCorr=0; aCorr tphiAlphaCorr ( fvc::flux ( @@ -16,71 +48,58 @@ ) + fvc::flux ( - -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); - Pair > vDotAlphal = - twoPhaseProperties->vDotAlphal(); - const volScalarField& vDotcAlphal = vDotAlphal[0](); - const volScalarField& vDotvAlphal = vDotAlphal[1](); + if (MULESCorr) + { + tphiAlphaCorr() -= tphiAlpha(); - volScalarField Sp - ( - IOobject + + volScalarField alpha100("alpha100", alpha10); + alpha10 = alpha1; + + MULES::correct ( - "Sp", - runTime.timeName(), - mesh - ), - vDotvAlphal - vDotcAlphal - ); + geometricOneField(), + alpha1, + tphiAlphaCorr(), + vDotvmcAlphal, + ( + divU*(alpha10 - alpha100) + - vDotvmcAlphal*alpha10 + )(), + 1, + 0 + ); - volScalarField Su - ( - IOobject + tphiAlpha() += tphiAlphaCorr(); + } + else + { + MULES::explicitSolve ( - "Su", - runTime.timeName(), - mesh - ), - // Divergence term is handled explicitly to be - // consistent with the explicit transport solution - divU*alpha1 - + vDotcAlphal - ); + geometricOneField(), + alpha1, + phi, + tphiAlphaCorr(), + vDotvmcAlphal, + (divU*alpha1 + vDotcAlphal)(), + 1, + 0 + ); - //MULES::explicitSolve - //( - // geometricOneField(), - // alpha1, - // phi, - // phiAlpha, - // Sp, - // Su, - // 1, - // 0 - //); + tphiAlpha = tphiAlphaCorr; + } - MULES::implicitSolve - ( - geometricOneField(), - alpha1, - phi, - phiAlpha, - Sp, - Su, - 1, - 0 - ); - - rhoPhi += - (runTime.deltaT()/totalDeltaT) - *(phiAlpha*(rho1 - rho2) + phi*rho2); + alpha2 = 1.0 - alpha1; } + rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2; + Info<< "Liquid phase volume fraction = " << alpha1.weightedAverage(mesh.V()).value() << " Min(alpha1) = " << min(alpha1).value() diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H index ad2cf243d2..a4338b907f 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H @@ -11,21 +11,18 @@ surfaceScalarField rhoPhi ); { - const dictionary& pimpleDict = pimple.dict(); - - label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr"))); - - label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); + #include "alphaControls.H" surfaceScalarField phic(mag(phi/mesh.magSf())); phic = min(interface.cAlpha()*phic, max(phic)); volScalarField divU(fvc::div(phi)); - dimensionedScalar totalDeltaT = runTime.deltaT(); - if (nAlphaSubCycles > 1) { + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi); + for ( subCycle alphaSubCycle(alpha1, nAlphaSubCycles); @@ -33,12 +30,15 @@ surfaceScalarField rhoPhi ) { #include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; } + + rhoPhi = rhoPhiSum; } else { #include "alphaEqn.H" } - rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2; + rho == alpha1*rho1 + alpha2*rho2; } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H index 60f15ab553..0e45fa495e 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H @@ -12,20 +12,6 @@ mesh ); - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -42,14 +28,19 @@ #include "createPhi.H" + Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl; autoPtr twoPhaseProperties = phaseChangeTwoPhaseMixture::New(U, phi); + volScalarField& alpha1(twoPhaseProperties->alpha1()); + volScalarField& alpha2(twoPhaseProperties->alpha2()); + const dimensionedScalar& rho1 = twoPhaseProperties->rho1(); const dimensionedScalar& rho2 = twoPhaseProperties->rho2(); const dimensionedScalar& pSat = twoPhaseProperties->pSat(); + // Need to store rho for ddt(rho, U) volScalarField rho ( @@ -60,11 +51,12 @@ mesh, IOobject::READ_IF_PRESENT ), - alpha1*rho1 + (scalar(1) - alpha1)*rho2, + alpha1*rho1 + alpha2*rho2, alpha1.boundaryField().types() ); rho.oldTime(); + // Construct interface from alpha1 distribution interfaceProperties interface(alpha1, U, twoPhaseProperties()); @@ -113,3 +105,5 @@ ); p_rgh = p - rho*gh; } + + fv::IOoptionList fvOptions(mesh); diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C index 00a7976a2e..bc6af055f0 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,6 +47,7 @@ Description #include "phaseChangeTwoPhaseMixture.H" #include "turbulenceModel.H" #include "pimpleControl.H" +#include "fvIOoptionList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -80,14 +81,10 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; + twoPhaseProperties->correct(); + #include "alphaEqnSubCycle.H" - - if (pimple.nCorrPIMPLE() == 1) - { - interface.correct(); - } - - turbulence->correct(); + interface.correct(); // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) @@ -99,9 +96,12 @@ int main(int argc, char *argv[]) { #include "pEqn.H" } - } - twoPhaseProperties->correct(); + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } runTime.write(); diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H index 817db5012c..7d9b71669a 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H @@ -46,6 +46,7 @@ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); + fvOptions.correct(U); } } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index acad49cd96..a19b6a3d9b 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -29,7 +29,10 @@ License #include "Time.H" #include "subCycle.H" #include "MULES.H" +#include "surfaceInterpolate.H" +#include "fvcGrad.H" #include "fvcSnGrad.H" +#include "fvcDiv.H" #include "fvcFlux.H" #include "fvcAverage.H" @@ -809,9 +812,8 @@ void Foam::multiphaseSystem::solve() const Time& runTime = mesh_.time(); - const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE"); - - label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); + const dictionary& alphaControls = mesh_.solverDict(phases_.first().name()); + label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C index fb6283d1ef..76b4dab8dd 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.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-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,7 +28,10 @@ License #include "Time.H" #include "subCycle.H" #include "MULES.H" +#include "surfaceInterpolate.H" +#include "fvcGrad.H" #include "fvcSnGrad.H" +#include "fvcDiv.H" #include "fvcFlux.H" // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * // @@ -246,15 +249,12 @@ void Foam::multiphaseMixture::solve() const Time& runTime = mesh_.time(); - const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE"); - - label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); - - scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha"))); - - volScalarField& alpha = phases_.first(); + const dictionary& alphaControls = mesh_.solverDict(alpha.name()); + label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); + scalar cAlpha(readScalar(alphaControls.lookup("cAlpha"))); + if (nAlphaSubCycles > 1) { surfaceScalarField rhoPhiSum(0.0*rhoPhi_); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H index 57c78027a4..428876fd22 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H @@ -1,5 +1,4 @@ -label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr"))); -label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles"))); +#include "alphaControls.H" if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H index 1abe3ef10d..813fdd8f10 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H @@ -1,5 +1,4 @@ #include "readTimeControls.H" + #include "alphaControls.H" - int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr"))); - int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles"))); Switch correctAlpha(pimple.dict().lookup("correctAlpha")); diff --git a/applications/test/laplacianFoam-communicators/Make/files b/applications/test/laplacianFoam-communicators/Make/files new file mode 100644 index 0000000000..67b30cabf3 --- /dev/null +++ b/applications/test/laplacianFoam-communicators/Make/files @@ -0,0 +1,3 @@ +laplacianFoam.C + +EXE = $(FOAM_USER_APPBIN)/laplacianFoam-communicators diff --git a/applications/test/laplacianFoam-communicators/Make/options b/applications/test/laplacianFoam-communicators/Make/options new file mode 100644 index 0000000000..d9745f69a0 --- /dev/null +++ b/applications/test/laplacianFoam-communicators/Make/options @@ -0,0 +1,3 @@ +EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = -lfiniteVolume diff --git a/applications/test/laplacianFoam-communicators/createFields.H b/applications/test/laplacianFoam-communicators/createFields.H new file mode 100644 index 0000000000..616afe1a88 --- /dev/null +++ b/applications/test/laplacianFoam-communicators/createFields.H @@ -0,0 +1,37 @@ + Info<< "Reading field T\n" << endl; + + volScalarField T + ( + IOobject + ( + "T", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Reading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + + Info<< "Reading diffusivity DT\n" << endl; + + dimensionedScalar DT + ( + transportProperties.lookup("DT") + ); diff --git a/applications/test/laplacianFoam-communicators/laplacianFoam.C b/applications/test/laplacianFoam-communicators/laplacianFoam.C new file mode 100644 index 0000000000..8a54cb256c --- /dev/null +++ b/applications/test/laplacianFoam-communicators/laplacianFoam.C @@ -0,0 +1,817 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + laplacianFoam + +Description + Solves a simple Laplace equation, e.g. for thermal diffusion in a solid. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "simpleControl.H" +#include "globalIndex.H" +#include "lduPrimitiveMesh.H" +#include "processorGAMGInterface.H" +#include "GAMGInterfaceField.H" +#include "processorLduInterfaceField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void checkUpperTriangular +( + const label size, + const labelUList& l, + const labelUList& u +) +{ + forAll(l, faceI) + { + if (u[faceI] < l[faceI]) + { + FatalErrorIn + ( + "checkUpperTriangular" + "(const label, const labelUList&, const labelUList&)" + ) << "Reversed face. Problem at face " << faceI + << " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError); + } + if (l[faceI] < 0 || u[faceI] < 0 || u[faceI] >= size) + { + FatalErrorIn + ( + "checkUpperTriangular" + "(const label, const labelUList&, const labelUList&)" + ) << "Illegal cell label. Problem at face " << faceI + << " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError); + } + } + + for (label faceI=1; faceI < l.size(); faceI++) + { + if (l[faceI-1] > l[faceI]) + { + FatalErrorIn + ( + "checkUpperTriangular" + "(const label, const labelUList&, const labelUList&)" + ) << "Lower not in incremental cell order." + << " Problem at face " << faceI + << " l:" << l[faceI] << " u:" << u[faceI] + << " previous l:" << l[faceI-1] << abort(FatalError); + } + else if (l[faceI-1] == l[faceI]) + { + // Same cell. + if (u[faceI-1] > u[faceI]) + { + FatalErrorIn + ( + "checkUpperTriangular" + "(const label, const labelUList&, const labelUList&)" + ) << "Upper not in incremental cell order." + << " Problem at face " << faceI + << " l:" << l[faceI] << " u:" << u[faceI] + << " previous u:" << u[faceI-1] << abort(FatalError); + } + } + } +} + + +void print(const string& msg, const lduMesh& mesh) +{ + const lduAddressing& addressing = mesh.lduAddr(); + const lduInterfacePtrsList interfaces = mesh.interfaces(); + + Pout<< "Mesh:" << msg.c_str() << nl + << " cells:" << addressing.size() << nl + << " faces:" << addressing.lowerAddr().size() << nl + << " patches:" << interfaces.size() << nl; + + + const labelUList& l = addressing.lowerAddr(); + const labelUList& startAddr = addressing.losortStartAddr(); + const labelUList& addr = addressing.losortAddr(); + + forAll(addressing, cellI) + { + Pout<< " cell:" << cellI << nl; + + label start = startAddr[cellI]; + label end = startAddr[cellI+1]; + + for (label index = start; index < end; index++) + { + Pout<< " nbr:" << l[addr[index]] << nl; + } + } + + Pout<< " Patches:" << nl; + forAll(interfaces, i) + { + if (interfaces.set(i)) + { + if (isA(interfaces[i])) + { + const processorLduInterface& pldui = + refCast(interfaces[i]); + Pout<< " " << i + << " me:" << pldui.myProcNo() + << " nbr:" << pldui.neighbProcNo() + << " comm:" << pldui.comm() + << " tag:" << pldui.tag() + << nl; + } + + { + Pout<< " " << i << " addressing:" << nl; + const labelUList& faceCells = interfaces[i].faceCells(); + forAll(faceCells, i) + { + Pout<< "\t\t" << i << '\t' << faceCells[i] << nl; + } + } + } + } +} + + +template +lduSchedule nonBlockingSchedule +( + const lduInterfacePtrsList& interfaces +) +{ + lduSchedule schedule(2*interfaces.size()); + label slotI = 0; + + forAll(interfaces, i) + { + if (interfaces.set(i) && !isA(interfaces[i])) + { + schedule[slotI].patch = i; + schedule[slotI].init = true; + slotI++; + schedule[slotI].patch = i; + schedule[slotI].init = false; + slotI++; + } + } + + forAll(interfaces, i) + { + if (interfaces.set(i) && isA(interfaces[i])) + { + schedule[slotI].patch = i; + schedule[slotI].init = true; + slotI++; + } + } + + forAll(interfaces, i) + { + if (interfaces.set(i) && isA(interfaces[i])) + { + schedule[slotI].patch = i; + schedule[slotI].init = false; + slotI++; + } + } + + return schedule; +} + + +void sendReceive +( + const label comm, + const label tag, + const globalIndex& offsets, + const scalarField& field, + + scalarField& allField +) +{ + label nProcs = Pstream::nProcs(comm); + + if (Pstream::master(comm)) + { + allField.setSize(offsets.size()); + + // Assign master slot + SubList + ( + allField, + offsets.localSize(0), + offsets.offset(0) + ).assign(field); + + // Assign slave slots + for (label procI = 1; procI < nProcs; procI++) + { + SubList procSlot + ( + allField, + offsets.localSize(procI), + offsets.offset(procI) + ); + + Pout<< "Receiving allField from " << procI + << " at offset:" << offsets.offset(procI) + << " size:" << offsets.size() + << endl; + + IPstream::read + ( + Pstream::nonBlocking, + procI, + reinterpret_cast(procSlot.begin()), + procSlot.byteSize(), + tag, + comm + ); + } + } + else + { + OPstream::write + ( + Pstream::nonBlocking, + 0, // master + reinterpret_cast(field.begin()), + field.byteSize(), + tag, + comm + ); + } +} + + +void sendReceive +( + const label comm, + const label tag, + const globalIndex& offsets, + const FieldField& field, + + FieldField& allField +) +{ + PstreamBuffers pBufs(Pstream::nonBlocking, Pstream::msgType(), comm); + + if (!Pstream::master(comm)) + { + UOPstream toMaster(Pstream::masterNo(), pBufs); + + Pout<< "To 0 sending " << field.size() + << " fields." << endl; + + forAll(field, intI) + { + toMaster << field[intI]; + } + } + pBufs.finishedSends(); + if (Pstream::master(comm)) + { + allField.setSize(offsets.size()); + forAll(allField, i) + { + allField.set(i, new scalarField(0)); + } + + // Insert master values + forAll(field, intI) + { + allField[intI] = field[intI]; + } + + + // Receive and insert slave values + label nProcs = Pstream::nProcs(comm); + + for (label procI = 1; procI < nProcs; procI++) + { + UIPstream fromSlave(procI, pBufs); + + label nSlaveInts = offsets.localSize(procI); + + Pout<< "From " << procI << " receiving " + << nSlaveInts << " fields." << endl; + + for (label intI = 0; intI < nSlaveInts; intI++) + { + label slotI = offsets.toGlobal(procI, intI); + + Pout<< " int:" << intI << " goes into slot " << slotI + << endl; + + fromSlave >> allField[slotI]; + } + } + } +} + + + +void collect +( + const label comm, + const globalIndex& cellOffsets, + const globalIndex& faceOffsets, + + const scalarField& diagonal, + const scalarField& upper, + const scalarField& lower, + + scalarField& allDiagonal, + scalarField& allUpper, + scalarField& allLower +) +{ + label nOutstanding = Pstream::nRequests(); + int allDiagonalTag = Pstream::allocateTag("allDiagonal:" __FILE__); + int allUpperTag = Pstream::allocateTag("allUpper:" __FILE__); + int allLowerTag = Pstream::allocateTag("allLower:" __FILE__); + + + sendReceive + ( + comm, + allDiagonalTag, + cellOffsets, + diagonal, + allDiagonal + ); + + sendReceive + ( + comm, + allUpperTag, + faceOffsets, + upper, + allUpper + ); + + sendReceive + ( + comm, + allLowerTag, + faceOffsets, + lower, + allLower + ); + + Pstream::waitRequests(nOutstanding); + + Pstream::freeTag("allDiagonal:" __FILE__, allDiagonalTag); + Pstream::freeTag("allUpper:" __FILE__, allUpperTag); + Pstream::freeTag("allLower:" __FILE__, allLowerTag); +} + + +void setCommunicator(fvMesh& mesh, const label newComm) +{ + const polyBoundaryMesh& pbm = mesh.boundaryMesh(); + + // The current state is consistent with the mesh so check where the new + // communicator is and adjust accordingly. + + forAll(pbm, patchI) + { + if (isA(pbm[patchI])) + { + processorPolyPatch& ppp = const_cast + ( + refCast + < + const processorPolyPatch + >(pbm[patchI]) + ); + + label thisRank = UPstream::procNo + ( + newComm, + ppp.comm(), + ppp.myProcNo() + ); + label nbrRank = UPstream::procNo + ( + newComm, + ppp.comm(), + ppp.neighbProcNo() + ); + + //ppp.comm() = newComm; + ppp.myProcNo() = thisRank; + ppp.neighbProcNo() = nbrRank; + } + } + mesh.polyMesh::comm() = newComm; +} + + +namespace Foam +{ + typedef UPtrList GAMGInterfaceFieldPtrsList; +} + + +// Gather matrices from processors procIDs[1..] on procIDs[0] +void gatherMatrices +( + const labelList& procIDs, + const PtrList& procMeshes, + + const lduMatrix& mat, + const FieldField& interfaceBouCoeffs, + const FieldField& interfaceIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + + PtrList& otherMats, + PtrList >& otherBouCoeffs, + PtrList >& otherIntCoeffs, + PtrList& otherInterfaces +) +{ + const label meshComm = mat.mesh().comm(); + + //lduInterfacePtrsList interfaces(mesh.interfaces()); + + if (Pstream::myProcNo(meshComm) == procIDs[0]) + { + // Master. + otherMats.setSize(procIDs.size()-1); + otherBouCoeffs.setSize(procIDs.size()-1); + otherIntCoeffs.setSize(procIDs.size()-1); + otherInterfaces.setSize(procIDs.size()-1); + + for (label i = 1; i < procIDs.size(); i++) + { + const lduMesh& procMesh = procMeshes[i]; + const lduInterfacePtrsList procInterfaces = procMesh.interfaces(); + + IPstream fromSlave + ( + Pstream::scheduled, + procIDs[i], + 0, // bufSize + Pstream::msgType(), + meshComm + ); + + otherMats.set(i-1, new lduMatrix(procMesh, fromSlave)); + + // Receive number of/valid interfaces + boolList validTransforms(fromSlave); + List validRanks(fromSlave); + + // Size coefficients + otherBouCoeffs.set + ( + i-1, + new FieldField(validTransforms.size()) + ); + otherIntCoeffs.set + ( + i-1, + new FieldField(validTransforms.size()) + ); + otherInterfaces.set + ( + i-1, + new GAMGInterfaceFieldPtrsList(validTransforms.size()) + ); + + forAll(validTransforms, intI) + { + if (validTransforms[intI]) + { + const processorGAMGInterface& interface = + refCast + ( + procInterfaces[intI] + ); + + + otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave)); + otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave)); + otherInterfaces[i-1].set + ( + intI, + GAMGInterfaceField::New + ( + interface, //procInterfaces[intI], + validTransforms[intI], + validRanks[intI] + ).ptr() + ); + } + } + } + } + else + { + // Send to master + OPstream toMaster + ( + Pstream::scheduled, + procIDs[0], + 0, + Pstream::msgType(), + meshComm + ); + + + // Count valid interfaces + boolList validTransforms(interfaceBouCoeffs.size(), false); + List validRanks(interfaceBouCoeffs.size(), -1); + forAll(interfaces, intI) + { + if (interfaces.set(intI)) + { + const processorLduInterfaceField& interface = + refCast + ( + interfaces[intI] + ); + + validTransforms[intI] = interface.doTransform(); + validRanks[intI] = interface.rank(); + } + } + + toMaster << mat << validTransforms << validRanks; + forAll(validTransforms, intI) + { + if (validTransforms[intI]) + { + toMaster + << interfaceBouCoeffs[intI] + << interfaceIntCoeffs[intI]; + } + } + } +} + + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "createFields.H" + + simpleControl simple(mesh); + + //const polyBoundaryMesh& pbm = mesh.boundaryMesh(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nCalculating temperature distribution\n" << endl; + + + // Get a subset of processors + labelList subProcs(3); + subProcs[0] = 0; + subProcs[1] = 1; + subProcs[2] = 2; + + + const UPstream::communicator newComm + ( + UPstream::worldComm, + subProcs, + true + ); + + + Pout<< "procIDs world :" << UPstream::procID(UPstream::worldComm) << endl; + Pout<< "procIDs newComm:" << UPstream::procID(newComm) << endl; + + +//// On ALL processors: get the interfaces: +//lduInterfacePtrsList interfaces(mesh.interfaces()); +//PtrList procMeshes; +// +//if (Pstream::myProcNo(newComm) != -1) +//{ +// print("InitialMesh", mesh); +// +// labelList procIDs(3); +// procIDs[0] = 0; +// procIDs[1] = 1; +// procIDs[2] = 2; +// +////XXXXXX +// // Collect meshes from procs 0,1 (in newComm) on 1. +// lduPrimitiveMesh::gather(mesh, procIDs, procMeshes); +// +// if (Pstream::myProcNo(newComm) == procIDs[0]) +// { +// // Print a bit +// forAll(procMeshes, i) +// { +// const lduMesh& pMesh = procMeshes[i]; +// print("procMesh" + Foam::name(i), pMesh); +// +// const lduAddressing& addr = pMesh.lduAddr(); +// checkUpperTriangular +// ( +// addr.size(), +// addr.lowerAddr(), +// addr.upperAddr() +// ); +// } +// +// +// // Combine +// labelList cellOffsets; +// labelListList faceMap; +// labelListList boundaryMap; +// labelListListList boundaryFaceMap; +// //autoPtr allMeshPtr = combineMeshes +// //( +// // newComm, +// // procIDs, +// // procMeshes, +// // +// // cellOffsets, // per mesh the starting cell +// // faceMap, // per mesh the starting face +// // boundaryMap, // per mesh,per interface the starting face +// // boundaryFaceMap +// //); +// //const lduPrimitiveMesh& allMesh = allMeshPtr(); +// const lduPrimitiveMesh allMesh +// ( +// newComm, +// procIDs, +// procMeshes, +// +// cellOffsets, +// faceMap, +// boundaryMap, +// boundaryFaceMap +// ); +// +// +// print("ALLMESH", allMesh); +// +// forAll(procMeshes, procMeshI) +// { +// const lduMesh& pMesh = procMeshes[procMeshI]; +// //const lduAddressing& pAddressing = pMesh.lduAddr(); +// +// Pout<< "procMesh" << procMeshI << endl +// << " cells start at:" << cellOffsets[procMeshI] << endl +// << " faces to to:" << faceMap[procMeshI] << endl; +// +// lduInterfacePtrsList interfaces = pMesh.interfaces(); +// forAll(interfaces, intI) +// { +// Pout<< " patch:" << intI +// << " becomes patch:" << boundaryMap[procMeshI][intI] +// << endl; +// +// Pout<< " patch:" << intI +// << " faces become faces:" +// << boundaryFaceMap[procMeshI][intI] +// << endl; +// } +// } +// } +// +// +// // Construct test data +// // ~~~~~~~~~~~~~~~~~~~ +// +// GAMGInterfaceFieldPtrsList interfaces(interfaces.size()); +// FieldField interfaceBouCoeffs(interfaces.size()); +// FieldField interfaceIntCoeffs(interfaces.size()); +// +// forAll(interfaces, intI) +// { +// if (interfaces.set(intI)) +// { +// label size = interfaces[intI].size(); +// +// interfaces.set +// ( +// intI, +// GAMGInterfaceField::New +// ( +// mesh.interfaces()[intI], +// interfaces[intI] +// ) +// ); +// interfaceBouCoeffs.set(intI, new scalarField(size, 111)); +// interfaceIntCoeffs.set(intI, new scalarField(size, 222)); +// } +// } +// +// +// PtrList otherMats; +// PtrList > otherBouCoeffs; +// PtrList > otherIntCoeffs; +// PtrList otherInterfaces; +// gatherMatrices +// ( +// procIDs, +// procMeshes, +// +// mat, +// interfaceBouCoeffs, +// interfaceIntCoeffs, +// interfaces, +// +// otherMats, +// otherBouCoeffs, +// otherIntCoeffs, +// otherInterfaces +// ); +////XXXXXX +//} + + + + { + Pout<< "World:" << UPstream::worldComm + << " procID:" << 2 + << " subComm:" << newComm + << " rank1:" << UPstream::procNo(newComm, UPstream::worldComm, 1) + << " rank2:" << UPstream::procNo(newComm, UPstream::worldComm, 2) + << endl; + } + + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix Teqn + ( + //fvm::ddt(T) - fvm::laplacian(DT, T) + fvm::laplacian(DT, T) + ); + + + { + label oldWarn = UPstream::warnComm; + UPstream::warnComm = newComm; + + label oldComm = mesh.comm(); + setCommunicator(mesh, newComm); + Pout<< "** oldcomm:" << oldComm + << " newComm:" << mesh.comm() << endl; + + if (Pstream::myProcNo(mesh.comm()) != -1) + { + solve(Teqn); + } + + setCommunicator(mesh, oldComm); + Pout<< "** reset mesh to:" << mesh.comm() << endl; + + UPstream::warnComm = oldWarn; + } + } + + #include "write.H" + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Pout<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/laplacianFoam-communicators/write.H b/applications/test/laplacianFoam-communicators/write.H new file mode 100644 index 0000000000..47aa182c0a --- /dev/null +++ b/applications/test/laplacianFoam-communicators/write.H @@ -0,0 +1,46 @@ + if (runTime.outputTime()) + { + volVectorField gradT(fvc::grad(T)); + + volScalarField gradTx + ( + IOobject + ( + "gradTx", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::X) + ); + + volScalarField gradTy + ( + IOobject + ( + "gradTy", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::Y) + ); + + volScalarField gradTz + ( + IOobject + ( + "gradTz", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + gradT.component(vector::Z) + ); + + + runTime.write(); + } diff --git a/applications/test/parallel-communicators/Make/files b/applications/test/parallel-communicators/Make/files new file mode 100644 index 0000000000..ab261d4da9 --- /dev/null +++ b/applications/test/parallel-communicators/Make/files @@ -0,0 +1,3 @@ +Test-parallel-communicators.C + +EXE = $(FOAM_USER_APPBIN)/Test-parallel-communicators diff --git a/applications/test/parallel-communicators/Make/options b/applications/test/parallel-communicators/Make/options new file mode 100644 index 0000000000..e69de29bb2 diff --git a/applications/test/parallel-communicators/Test-parallel-communicators.C b/applications/test/parallel-communicators/Test-parallel-communicators.C new file mode 100644 index 0000000000..2b9277fb6e --- /dev/null +++ b/applications/test/parallel-communicators/Test-parallel-communicators.C @@ -0,0 +1,203 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + Test-parallel-communicators + +Description + Checks communication using user-defined communicators + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "IPstream.H" +#include "OPstream.H" +#include "vector.H" +#include "IOstreams.H" +#include "PstreamReduceOps.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +scalar sumReduce +( + const label comm, + const scalar localValue +) +{ + scalar sum; + if (Pstream::parRun()) + { + if (UPstream::master(comm)) + { + // Add master value and all slaves + sum = localValue; + + for + ( + int slave=Pstream::firstSlave(); + slave<=Pstream::lastSlave(comm); + slave++ + ) + { + scalar slaveValue; + UIPstream::read + ( + Pstream::blocking, + slave, + reinterpret_cast(&slaveValue), + sizeof(scalar), + UPstream::msgType(), // tag + comm // communicator + ); + + sum += slaveValue; + } + + // Send back to slaves + + for + ( + int slave=UPstream::firstSlave(); + slave<=UPstream::lastSlave(comm); + slave++ + ) + { + UOPstream::write + ( + UPstream::blocking, + slave, + reinterpret_cast(&sum), + sizeof(scalar), + UPstream::msgType(), // tag + comm // communicator + ); + } + } + else + { + { + UOPstream::write + ( + UPstream::blocking, + UPstream::masterNo(), + reinterpret_cast(&localValue), + sizeof(scalar), + UPstream::msgType(), // tag + comm // communicator + ); + } + + { + UIPstream::read + ( + UPstream::blocking, + UPstream::masterNo(), + reinterpret_cast(&sum), + sizeof(scalar), + UPstream::msgType(), // tag + comm // communicator + ); + } + } + } + return sum; +} + + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + // Allocate a communicator + label n = Pstream::nProcs(UPstream::worldComm); + + DynamicList