diff --git a/applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H b/applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H index 049b6c0b3610160c88fe1f6942456f34addc8fec..1ef67b292e159c0b75283bda44d2af47043d85f0 100644 --- a/applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H +++ b/applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H @@ -55,7 +55,7 @@ ( dimensionedScalar::lookupOrAddToDict ( - "alphaEps", + "alphak", kEpsilonDict, 1.0 ) diff --git a/applications/solvers/multiphase/interFoam/Allwclean b/applications/solvers/multiphase/interFoam/Allwclean index 350d4b268b070f686f65a6e959f2f5abdfde608d..2c54bfa3218478eef96fce9e9bbc8f23612add83 100755 --- a/applications/solvers/multiphase/interFoam/Allwclean +++ b/applications/solvers/multiphase/interFoam/Allwclean @@ -6,5 +6,6 @@ wclean wclean interDyMFoam wclean MRFInterFoam wclean porousInterFoam +wclean LTSInterFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interFoam/Allwmake b/applications/solvers/multiphase/interFoam/Allwmake index b25e4440e60bc084516b2829565268642a7ea2cc..8044426582eacdd44b4277674eb8d715c0365ed0 100755 --- a/applications/solvers/multiphase/interFoam/Allwmake +++ b/applications/solvers/multiphase/interFoam/Allwmake @@ -6,5 +6,6 @@ wmake wmake interDyMFoam wmake MRFInterFoam wmake porousInterFoam +wmake LTSInterFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C new file mode 100644 index 0000000000000000000000000000000000000000..9fcb19a7cc242f590087ea7b62ad268606516e43 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 + interFoam + +Description + Solver for 2 incompressible, isothermal immiscible fluids using a VOF + (volume of fluid) phase-fraction based interface capturing approach. + + The momentum and other fluid properties are of the "mixture" and a single + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + For a two-fluid approach see twoPhaseEulerFoam. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "MULES.H" +#include "subCycle.H" +#include "interfaceProperties.H" +#include "twoPhaseMixture.H" +#include "turbulenceModel.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readPISOControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "correctPhi.H" + #include "CourantNo.H" + #include "setInitialrDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readPISOControls.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + #include "setrDeltaT.H" + + twoPhaseProperties.correct(); + + #include "alphaEqnSubCycle.H" + turbulence->correct(); + #include "UEqn.H" + + // --- PISO loop + for (int corr=0; corr<nCorr; corr++) + { + #include "pEqn.H" + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C new file mode 100644 index 0000000000000000000000000000000000000000..e281ae575cbfd41e0000cb1f2c1d178417047c63 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 "MULES.H" +#include "upwind.H" +#include "uncorrectedSnGrad.H" +#include "gaussConvectionScheme.H" +#include "gaussLaplacianScheme.H" +#include "uncorrectedSnGrad.H" +#include "surfaceInterpolate.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "syncTools.H" + +#include "fvCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void Foam::MULES::explicitLTSSolve +( + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +) +{ + explicitLTSSolve + ( + geometricOneField(), + psi, + phi, + phiPsi, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +void Foam::MULES::implicitSolve +( + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +) +{ + implicitSolve + ( + geometricOneField(), + psi, + phi, + phiPsi, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H new file mode 100644 index 0000000000000000000000000000000000000000..729b041aa62d0b7870d28a613ab3f1a332d96bae --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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/>. + +Global + MULES + +Description + MULES: Multidimensional universal limiter with explicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter. + + Parameters are the variable to solve, the normal convective flux and the + actual explicit flux of the variable which is also used to return limited + flux used in the bounded-solution. + +SourceFiles + MULES.C + +\*---------------------------------------------------------------------------*/ + +#ifndef MULES_H +#define MULES_H + +#include "volFields.H" +#include "surfaceFieldsFwd.H" +#include "primitiveFieldsFwd.H" +#include "zeroField.H" +#include "geometricOneField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace MULES +{ + +template<class RhoType, class SpType, class SuType> +void explicitLTSSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void explicitLTSSolve +( + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +); + +template<class RhoType, class SpType, class SuType> +void implicitSolve +( + const RhoType& rho, + volScalarField& gamma, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void implicitSolve +( + volScalarField& gamma, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + +template<class RhoType, class SpType, class SuType> +void limiter +( + scalarField& allLambda, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiBD, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + +} // End namespace MULES + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "MULESTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..a6d3a3f6b246e2463979604d84ad33dd108c5c01 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C @@ -0,0 +1,602 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 "MULES.H" +#include "upwind.H" +#include "uncorrectedSnGrad.H" +#include "gaussConvectionScheme.H" +#include "gaussLaplacianScheme.H" +#include "uncorrectedSnGrad.H" +#include "surfaceInterpolate.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "syncTools.H" + +#include "fvCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::explicitLTSSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + Info<< "MULES: Solving for " << psi.name() << endl; + + const fvMesh& mesh = psi.mesh(); + psi.correctBoundaryConditions(); + + surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi); + + surfaceScalarField& phiCorr = phiPsi; + phiCorr -= phiBD; + + scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + limiter + ( + allLambda, + rho, + psi, + phiBD, + phiCorr, + Sp.field(), + Su.field(), + psiMax, + psiMin, + 3 + ); + + phiPsi = phiBD + lambda*phiCorr; + + scalarField& psiIf = psi; + const scalarField& psi0 = psi.oldTime(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT"); + + psiIf = 0.0; + fvc::surfaceIntegrate(psiIf, phiPsi); + + if (mesh.moving()) + { + psiIf = + ( + mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc() + + Su.field() + - psiIf + )/(rho*rDeltaT - Sp.field()); + } + else + { + psiIf = + ( + rho.oldTime()*psi0*rDeltaT + + Su.field() + - psiIf + )/(rho*rDeltaT - Sp.field()); + } + + psi.correctBoundaryConditions(); +} + + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::implicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label maxIter + ( + readLabel(MULEScontrols.lookup("maxIter")) + ); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + scalar maxUnboundedness + ( + readScalar(MULEScontrols.lookup("maxUnboundedness")) + ); + + scalar CoCoeff + ( + readScalar(MULEScontrols.lookup("CoCoeff")) + ); + + scalarField allCoLambda(mesh.nFaces()); + + { + surfaceScalarField Cof = + mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs() + *mag(phi)/mesh.magSf(); + + slicedSurfaceScalarField CoLambda + ( + IOobject + ( + "CoLambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allCoLambda, + false // Use slices for the couples + ); + + CoLambda == 1.0/max(CoCoeff*Cof, scalar(1)); + } + + scalarField allLambda(allCoLambda); + //scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + linear<scalar> CDs(mesh); + upwind<scalar> UDs(mesh, phi); + //fv::uncorrectedSnGrad<scalar> snGrads(mesh); + + fvScalarMatrix psiConvectionDiffusion + ( + fvm::ddt(rho, psi) + + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi) + //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads) + //.fvmLaplacian(Dpsif, psi) + - fvm::Sp(Sp, psi) + - Su + ); + + surfaceScalarField phiBD = psiConvectionDiffusion.flux(); + + surfaceScalarField& phiCorr = phiPsi; + phiCorr -= phiBD; + + for (label i=0; i<maxIter; i++) + { + if (i != 0 && i < 4) + { + allLambda = allCoLambda; + } + + limiter + ( + allLambda, + rho, + psi, + phiBD, + phiCorr, + Sp.field(), + Su.field(), + psiMax, + psiMin, + nLimiterIter + ); + + solve + ( + psiConvectionDiffusion + fvc::div(lambda*phiCorr), + MULEScontrols + ); + + scalar maxPsiM1 = gMax(psi.internalField()) - 1.0; + scalar minPsi = gMin(psi.internalField()); + + scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0)); + + if (unboundedness < maxUnboundedness) + { + break; + } + else + { + Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1 + << " min(" << psi.name() << ") = " << minPsi << endl; + + phiBD = psiConvectionDiffusion.flux(); + + /* + word gammaScheme("div(phi,gamma)"); + word gammarScheme("div(phirb,gamma)"); + + const surfaceScalarField& phir = + mesh.lookupObject<surfaceScalarField>("phir"); + + phiCorr = + fvc::flux + ( + phi, + psi, + gammaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - psi, gammarScheme), + psi, + gammarScheme + ) + - phiBD; + */ + } + } + + phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr; +} + + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::limiter +( + scalarField& allLambda, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiBD, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const scalarField& psiIf = psi; + const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); + + const scalarField& psi0 = psi.oldTime(); + + const fvMesh& mesh = psi.mesh(); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighb = mesh.neighbour(); + tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc(); + const scalarField& V = tVsc(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT"); + + const scalarField& phiBDIf = phiBD; + const surfaceScalarField::GeometricBoundaryField& phiBDBf = + phiBD.boundaryField(); + + const scalarField& phiCorrIf = phiCorr; + const surfaceScalarField::GeometricBoundaryField& phiCorrBf = + phiCorr.boundaryField(); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + scalarField& lambdaIf = lambda; + surfaceScalarField::GeometricBoundaryField& lambdaBf = + lambda.boundaryField(); + + scalarField psiMaxn(psiIf.size(), psiMin); + scalarField psiMinn(psiIf.size(), psiMax); + + scalarField sumPhiBD(psiIf.size(), 0.0); + + scalarField sumPhip(psiIf.size(), VSMALL); + scalarField mSumPhim(psiIf.size(), VSMALL); + + forAll(phiCorrIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); + psiMinn[own] = min(psiMinn[own], psiIf[nei]); + + psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); + psiMinn[nei] = min(psiMinn[nei], psiIf[own]); + + sumPhiBD[own] += phiBDIf[facei]; + sumPhiBD[nei] -= phiBDIf[facei]; + + scalar phiCorrf = phiCorrIf[facei]; + + if (phiCorrf > 0.0) + { + sumPhip[own] += phiCorrf; + mSumPhim[nei] += phiCorrf; + } + else + { + mSumPhim[own] -= phiCorrf; + sumPhip[nei] -= phiCorrf; + } + } + + forAll(phiCorrBf, patchi) + { + const fvPatchScalarField& psiPf = psiBf[patchi]; + const scalarField& phiBDPf = phiBDBf[patchi]; + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + if (psiPf.coupled()) + { + scalarField psiPNf = psiPf.patchNeighbourField(); + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); + } + } + else + { + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); + } + } + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + sumPhiBD[pfCelli] += phiBDPf[pFacei]; + + scalar phiCorrf = phiCorrPf[pFacei]; + + if (phiCorrf > 0.0) + { + sumPhip[pfCelli] += phiCorrf; + } + else + { + mSumPhim[pfCelli] -= phiCorrf; + } + } + } + + psiMaxn = min(psiMaxn, psiMax); + psiMinn = max(psiMinn, psiMin); + + //scalar smooth = 0.5; + //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); + //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); + + if (mesh.moving()) + { + tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0(); + + psiMaxn = + V*((rho*rDeltaT - Sp)*psiMaxn - Su) + - (V0()*rDeltaT)*rho.oldTime()*psi0 + + sumPhiBD; + + psiMinn = + V*(Su - (rho*rDeltaT - Sp)*psiMinn) + + (V0*rDeltaT)*rho.oldTime()*psi0 + - sumPhiBD; + } + else + { + psiMaxn = + V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su) + + sumPhiBD; + + psiMinn = + V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su) + - sumPhiBD; + } + + scalarField sumlPhip(psiIf.size()); + scalarField mSumlPhim(psiIf.size()); + + for(int j=0; j<nLimiterIter; j++) + { + sumlPhip = 0.0; + mSumlPhim = 0.0; + + forAll(lambdaIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[own] += lambdaPhiCorrf; + mSumlPhim[nei] += lambdaPhiCorrf; + } + else + { + mSumlPhim[own] -= lambdaPhiCorrf; + sumlPhip[nei] -= lambdaPhiCorrf; + } + } + + forAll(lambdaBf, patchi) + { + scalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[pfCelli] += lambdaPhiCorrf; + } + else + { + mSumlPhim[pfCelli] -= lambdaPhiCorrf; + } + } + } + + forAll (sumlPhip, celli) + { + sumlPhip[celli] = + max(min + ( + (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli], + 1.0), 0.0 + ); + + mSumlPhim[celli] = + max(min + ( + (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli], + 1.0), 0.0 + ); + } + + const scalarField& lambdam = sumlPhip; + const scalarField& lambdap = mSumlPhim; + + forAll(lambdaIf, facei) + { + if (phiCorrIf[facei] > 0.0) + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdap[owner[facei]], lambdam[neighb[facei]]) + ); + } + else + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdam[owner[facei]], lambdap[neighb[facei]]) + ); + } + } + + + forAll(lambdaBf, patchi) + { + fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + + syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>()); + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..205812e978dd9dfe7ace29f84753167d1186ab39 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files @@ -0,0 +1,4 @@ +LTSInterFoam.C +MULES.C + +EXE = $(FOAM_APPBIN)/LTSInterFoam diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..24349f694e016e85732d7896b7b98c57083eaa66 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options @@ -0,0 +1,15 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -ltwoPhaseInterfaceProperties \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H new file mode 100644 index 0000000000000000000000000000000000000000..a4cb6e7688c168c92d3ac1aebd5cfb1cddf93382 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H @@ -0,0 +1,36 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phic = mag(phi/mesh.magSf()); + phic = min(interface.cAlpha()*phic, max(phic)); + surfaceScalarField phir = phic*interface.nHatf(); + + for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) + { + surfaceScalarField phiAlpha = + fvc::flux + ( + phi, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + alpha1, + alpharScheme + ); + + MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0); + //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); + + rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; + } + + Info<< "Liquid phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; +} diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H new file mode 100644 index 0000000000000000000000000000000000000000..9aae37a8bef50329d7e2b4cac649219d48dcc707 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H @@ -0,0 +1,35 @@ +label nAlphaCorr +( + readLabel(piso.lookup("nAlphaCorr")) +); + +label nAlphaSubCycles +( + readLabel(piso.lookup("nAlphaSubCycles")) +); + +if (nAlphaSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + + for + ( + subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { +# include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ +# include "alphaEqn.H" +} + +interface.correct(); + +rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2; diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H new file mode 100644 index 0000000000000000000000000000000000000000..1f2f082fb36b8dc2f4538b2714990c1912482107 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H @@ -0,0 +1,30 @@ +scalar maxDeltaT +( + piso.lookupOrDefault<scalar>("maxDeltaT", GREAT) +); + +volScalarField rDeltaT +( + IOobject + ( + "rDeltaT", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT) +); + +volScalarField rSubDeltaT +( + IOobject + ( + "rSubDeltaT", + runTime.timeName(), + mesh + ), + mesh, + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT) +); diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H new file mode 100644 index 0000000000000000000000000000000000000000..f5b8450ddd9d1032784bf7f96eed44932fc54e44 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -0,0 +1,109 @@ +{ + scalar maxCo + ( + piso.lookupOrDefault<scalar>("maxCo", 0.9) + ); + + scalar maxAlphaCo + ( + piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2) + ); + + scalar rDeltaTSmoothingCoeff + ( + piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1) + ); + + label nAlphaSpreadIter + ( + piso.lookupOrDefault<label>("nAlphaSpreadIter", 1) + ); + + label nAlphaSweepIter + ( + piso.lookupOrDefault<label>("nAlphaSweepIter", 5) + ); + + scalar rDeltaTDampingCoeff + ( + piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0) + ); + + scalar maxDeltaT + ( + piso.lookupOrDefault<scalar>("maxDeltaT", GREAT) + ); + + volScalarField rDeltaT0 = rDeltaT; + + // Set the reciprocal time-step using an effective maximum Courant number + rDeltaT = max + ( + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT), + fvc::surfaceSum + ( + mag(rhoPhi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf()) + )/rho + ); + + // Limit the time-step further in the region of the interface + { + surfaceScalarField alphaf = fvc::interpolate(alpha1); + + surfaceScalarField SfUfbyDelta = + pos(alphaf - 0.01)*pos(0.99 - alphaf) + *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi); + + rDeltaT = max + ( + rDeltaT, + fvc::surfaceSum(mag(SfUfbyDelta/(maxAlphaCo*mesh.magSf()))) + ); + } + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + if (rDeltaTSmoothingCoeff < 1.0) + { + fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); + } + + if (nAlphaSpreadIter > 0) + { + fvc::spread(rDeltaT, alpha1, nAlphaSpreadIter); + } + + if (nAlphaSweepIter > 0) + { + fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter); + } + + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + // Limit rate of change of time scale + // - reduce as much as required + // - only increase at a fraction of old time scale + if + ( + rDeltaTDampingCoeff < 1.0 + && runTime.timeIndex() > runTime.startTimeIndex() + 1 + ) + { + Info<< "Damping rDeltaT" << endl; + rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff); + } + + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + label nAlphaSubCycles + ( + readLabel(piso.lookup("nAlphaSubCycles")) + ); + + rSubDeltaT = rDeltaT*nAlphaSubCycles; +} diff --git a/applications/utilities/mesh/manipulation/subsetMesh/Make/options b/applications/utilities/mesh/manipulation/subsetMesh/Make/options index d27c95d033dd5d7b1995c8ff8dc406e35ca1f586..969020c4afaf5d784299462b9e1af282040ba6b4 100644 --- a/applications/utilities/mesh/manipulation/subsetMesh/Make/options +++ b/applications/utilities/mesh/manipulation/subsetMesh/Make/options @@ -4,4 +4,5 @@ EXE_INC = \ EXE_LIBS = \ -lfiniteVolume \ - -lmeshTools + -lmeshTools \ + -lgenericPatchFields diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index 6018af2e14e272075a4ac3c727da08b2bf885edf..41f73e78a16174f543b6d985e99e95e574d935ed 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -37,8 +37,11 @@ namespace Foam const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70}; const scalar - SIBS::safe1 = 0.25, SIBS::safe2 = 0.7, - SIBS::redMax = 1.0e-5, SIBS::redMin = 0.0, SIBS::scaleMX = 0.1; + SIBS::safe1 = 0.25, + SIBS::safe2 = 0.7, + SIBS::redMax = 1.0e-5, + SIBS::redMin = 0.7, + SIBS::scaleMX = 0.1; }; diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 21b90fa409a32cf2acb74d8e419c8c6cc4133d60..dc325c1ddf89e2bab343bdf7dda360096099a462 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -337,6 +337,7 @@ $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C finiteVolume/fvc/fvcMeshPhi.C +finiteVolume/fvc/fvcSmooth/fvcSmooth.C general = cfdTools/general $(general)/findRefCell/findRefCell.C diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C new file mode 100644 index 0000000000000000000000000000000000000000..47cfb5582ac8d84adba7b67e36c2d1a2d0524752 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C @@ -0,0 +1,317 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 "fvcSmooth.H" +#include "volFields.H" +#include "FaceCellWave.H" +#include "smoothData.H" +#include "sweepData.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::scalar Foam::smoothData::maxRatio = 1.0; + +void Foam::fvc::smooth +( + volScalarField& field, + const scalar coeff +) +{ + const fvMesh& mesh = field.mesh(); + smoothData::maxRatio = 1 + coeff; + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<smoothData> changedFacesInfo(changedFaces.size()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + // Check if owner value much larger than neighbour value or vice versa + if (field[own] > smoothData::maxRatio*field[nbr]) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + else if (field[nbr] > smoothData::maxRatio*field[own]) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[nbr])); + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Set initial field on cells + List<smoothData> cellData(mesh.nCells()); + + forAll(field, celli) + { + cellData[celli] = field[celli]; + } + + // Set initial field on faces + List<smoothData> faceData(mesh.nFaces()); + + + // Propagate information over whole domain + FaceCellWave<smoothData > smoothData + ( + mesh, + changedFaces, + changedFacesInfo, + faceData, + cellData, + mesh.globalData().nTotalCells() // max iterations + ); + + forAll(field, celli) + { + field[celli] = cellData[celli].value(); + } + + field.correctBoundaryConditions(); +} + + +void Foam::fvc::spread +( + volScalarField& field, + const volScalarField& alpha, + const label nLayers +) +{ + const fvMesh& mesh = field.mesh(); + smoothData::maxRatio = 1; + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<smoothData> changedFacesInfo(changedFaces.size()); + + // Set initial field on cells + List<smoothData> cellData(mesh.nCells()); + + forAll(field, celli) + { + cellData[celli] = field[celli]; + } + + // Set initial field on faces + List<smoothData> faceData(mesh.nFaces()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + if + ( + (alpha[own] > 0.01 && alpha[own] < 0.99) + || (alpha[nbr] > 0.01 && alpha[nbr] < 0.99) + ) + { + if (mag(alpha[own] - alpha[nbr]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + smoothData(max(field[own], field[nbr])) + ); + } + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + scalarField alphapn = + alpha.boundaryField()[patchi].patchNeighbourField(); + + if + ( + (alpha[own] > 0.01 && alpha[own] < 0.99) + || (alphapn[patchFacei] > 0.01 && alphapn[patchFacei] < 0.99) + ) + { + if (mag(alpha[own] - alphapn[patchFacei]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + } + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Propagate information over whole domain + FaceCellWave<smoothData> smoothData + ( + mesh, + faceData, + cellData + ); + + smoothData.setFaceInfo(changedFaces, changedFacesInfo); + + smoothData.iterate(nLayers); + + forAll(field, celli) + { + field[celli] = cellData[celli].value(); + } + + field.correctBoundaryConditions(); +} + + +void Foam::fvc::sweep +( + volScalarField& field, + const volScalarField& alpha, + const label nLayers +) +{ + const fvMesh& mesh = field.mesh(); + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<sweepData> changedFacesInfo(changedFaces.size()); + + // Set initial field on cells + List<sweepData> cellData(mesh.nCells()); + + // Set initial field on faces + List<sweepData> faceData(mesh.nFaces()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + const vectorField& Cf = mesh.faceCentres(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + if (mag(alpha[own] - alpha[nbr]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + sweepData(max(field[own], field[nbr]), Cf[facei]) + ); + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + scalarField alphapn = + alpha.boundaryField()[patchi].patchNeighbourField(); + + if (mag(alpha[own] - alphapn[patchFacei]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + sweepData(field[own], Cf[facei]) + ); + } + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Propagate information over whole domain + FaceCellWave<sweepData> sweepData + ( + mesh, + faceData, + cellData + ); + + sweepData.setFaceInfo(changedFaces, changedFacesInfo); + + sweepData.iterate(nLayers); + + forAll(field, celli) + { + if (cellData[celli].valid()) + { + field[celli] = max(field[celli], cellData[celli].value()); + } + } + + field.correctBoundaryConditions(); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H new file mode 100644 index 0000000000000000000000000000000000000000..2ce5b1b5d321b62014820853742327adcb61ec64 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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/>. + +InNamespace + Foam::fvc + +Description + Provides functions smooth spread and sweep which use the FaceCellWave + algorithm to smooth and redistribute the first field argument. + + smooth: smooths the field by ensuring the values in neighbouring + cells are at least coeff* the cell value. + + spread: redistributes the field by spreading the maximum value within + the region defined by the value and gradient of alpha. + + sweep: redistributes the field by sweeping the maximum value where the + gradient of alpha is large away from that starting point of the + sweep. + +SourceFiles + fvcSmooth.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fvcSmooth_H +#define fvcSmooth_H + +#include "volFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fvc +{ + void smooth + ( + volScalarField& field, + const scalar coeff + ); + + void spread + ( + volScalarField& field, + const volScalarField& alpha, + const label nLayers + ); + + void sweep + ( + volScalarField& field, + const volScalarField& alpha, + const label nLayers + ); +} +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H new file mode 100644 index 0000000000000000000000000000000000000000..d027f9f2a2479e38737014d62a3ca78adfcc23c2 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H @@ -0,0 +1,204 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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::smoothData + +Description + Helper class used by the fvc::smooth and fvc::spread functions. + +SourceFiles + smoothData.H + smoothDataI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef smoothData_H +#define smoothData_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class smoothData Declaration +\*---------------------------------------------------------------------------*/ + +class smoothData +{ + scalar value_; + + // Private Member Functions + + //- Update - gets information from neighbouring face/cell and + // uses this to update itself (if necessary) and return true + inline bool update + ( + const smoothData& svf, + const scalar scale, + const scalar tol + ); + + +public: + + // Static data members + + //- Field fraction + static scalar maxRatio; + + + // Constructors + + //- Construct null + inline smoothData(); + + //- Construct from inverse field value + inline smoothData(const scalar value); + + + // Member Functions + + // Access + + //- Return value + scalar value() const + { + return value_; + } + + + // Needed by FaceCellWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value + inline bool valid() const; + + //- Check for identical geometrical data + // Used for cyclics checking + inline bool sameGeometry + ( + const polyMesh&, + const smoothData&, + const scalar + ) const; + + //- Convert any absolute coordinates into relative to + // (patch)face centre + inline void leaveDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Reverse of leaveDomain + inline void enterDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Apply rotation matrix to any coordinates + inline void transform + ( + const polyMesh&, + const tensor& + ); + + //- Influence of neighbouring face + inline bool updateCell + ( + const polyMesh&, + const label thisCellI, + const label neighbourFaceI, + const smoothData& svf, + const scalar tol + ); + + //- Influence of neighbouring cell + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const label neighbourCellI, + const smoothData& svf, + const scalar tol + ); + + //- Influence of different value on same face + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const smoothData& svf, + const scalar tol + ); + + + // Member Operators + + inline void operator=(const scalar value); + + // Needed for List IO + inline bool operator==(const smoothData&) const; + + inline bool operator!=(const smoothData&) const; + + + // IOstream Operators + + friend Ostream& operator<< + ( + Ostream& os, + const smoothData& svf + ) + { + return os << svf.value_; + } + + friend Istream& operator>>(Istream& is, smoothData& svf) + { + return is >> svf.value_; + } +}; + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "smoothDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H new file mode 100644 index 0000000000000000000000000000000000000000..3faaa694b617e5fbea9a744a2b67afde6a9992af --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H @@ -0,0 +1,192 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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/>. + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline bool Foam::smoothData::update +( + const smoothData::smoothData& svf, + const scalar scale, + const scalar tol +) +{ + if (!valid() || (value_ < VSMALL)) + { + // My value not set - take over neighbour + value_ = svf.value()/scale; + + // Something changed - let caller know + return true; + } + else if (svf.value() > (1 + tol)*scale*value_) + { + // Neighbour is too big for me - Up my value + value_ = svf.value()/scale; + + // Something changed - let caller know + return true; + } + else + { + // Neighbour is not too big for me or change is too small + // Nothing changed + return false; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::smoothData::smoothData() +: + value_(-GREAT) +{} + + +inline Foam::smoothData::smoothData(const scalar value) +: + value_(value) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline bool Foam::smoothData::valid() const +{ + return value_ > -SMALL; +} + + +inline bool Foam::smoothData::sameGeometry +( + const polyMesh&, + const smoothData&, + const scalar +) const +{ + return true; +} + + +inline void Foam::smoothData::leaveDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& +) +{} + + +inline void Foam::smoothData::transform +( + const polyMesh&, + const tensor& +) +{} + + +inline void Foam::smoothData::enterDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& +) +{} + + +inline bool Foam::smoothData::updateCell +( + const polyMesh&, + const label, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over info from face if more than deltaRatio larger + return update(svf, maxRatio, tol); +} + + +inline bool Foam::smoothData::updateFace +( + const polyMesh&, + const label, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over information from cell without any scaling (scale = 1.0) + return update(svf, 1.0, tol); +} + + +// Update this (face) with coupled face information. +inline bool Foam::smoothData::updateFace +( + const polyMesh&, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over information from coupled face without any scaling (scale = 1.0) + return update(svf, 1.0, tol); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline void Foam::smoothData::operator= +( + const scalar value +) +{ + value_ = value; +} + + +inline bool Foam::smoothData::operator== +( + const smoothData& rhs +) const +{ + return value_ == rhs.value(); +} + + +inline bool Foam::smoothData::operator!= +( + const smoothData& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H new file mode 100644 index 0000000000000000000000000000000000000000..1bfec88fdfcc67283395632746ed5dff2aa4b070 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H @@ -0,0 +1,205 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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::sweepData + +Description + Helper class used by fvc::sweep function. + +SourceFiles + sweepData.H + sweepDataI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef sweepData_H +#define sweepData_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class sweepData Declaration +\*---------------------------------------------------------------------------*/ + +class sweepData +{ + scalar value_; + point origin_; + + // Private Member Functions + + //- Update - gets information from neighbouring face/cell and + // uses this to update itself (if necessary) and return true + inline bool update + ( + const sweepData& svf, + const point& position, + const scalar tol + ); + + +public: + + // Constructors + + //- Construct null + inline sweepData(); + + //- Construct from component + inline sweepData(const scalar value, const point& origin); + + + // Member Functions + + // Access + + //- Return value + scalar value() const + { + return value_; + } + + //- Return origin + const point& origin() const + { + return origin_; + } + + + // Needed by FaceCellWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value + inline bool valid() const; + + //- Check for identical geometrical data + // Used for cyclics checking + inline bool sameGeometry + ( + const polyMesh&, + const sweepData&, + const scalar + ) const; + + //- Convert any absolute coordinates into relative to + // (patch)face centre + inline void leaveDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Reverse of leaveDomain + inline void enterDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Apply rotation matrix to any coordinates + inline void transform + ( + const polyMesh&, + const tensor& + ); + + //- Influence of neighbouring face + inline bool updateCell + ( + const polyMesh&, + const label thisCellI, + const label neighbourFaceI, + const sweepData& svf, + const scalar tol + ); + + //- Influence of neighbouring cell + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const label neighbourCellI, + const sweepData& svf, + const scalar tol + ); + + //- Influence of different value on same face + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const sweepData& svf, + const scalar tol + ); + + + // Member Operators + + inline void operator=(const scalar value); + + // Needed for List IO + inline bool operator==(const sweepData&) const; + + inline bool operator!=(const sweepData&) const; + + + // IOstream Operators + + friend Ostream& operator<< + ( + Ostream& os, + const sweepData& svf + ) + { + return os << svf.value_ << svf.origin_; + } + + friend Istream& operator>>(Istream& is, sweepData& svf) + { + return is >> svf.value_ >> svf.origin_; + } +}; + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "sweepDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H new file mode 100644 index 0000000000000000000000000000000000000000..651c5e4e5260252251da7cd445c056d419d8965c --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H @@ -0,0 +1,208 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 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 "transform.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline bool Foam::sweepData::update +( + const sweepData::sweepData& svf, + const point& position, + const scalar tol +) +{ + if (!valid()) + { + operator=(svf); + return true; + } + + scalar myDist2 = magSqr(position - origin()); + + if (myDist2 < SMALL) + { + if (svf.value() > value()) + { + operator=(svf); + return true; + } + else + { + return false; + } + } + + scalar dist2 = magSqr(position - svf.origin()); + + if (dist2 < myDist2) + { + operator=(svf); + return true; + } + + return false; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::sweepData::sweepData() +: + value_(-GREAT), + origin_(vector::max) +{} + + +inline Foam::sweepData::sweepData(const scalar value, const point& origin) +: + value_(value), + origin_(origin) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline bool Foam::sweepData::valid() const +{ + return value_ > -SMALL; +} + + +inline bool Foam::sweepData::sameGeometry +( + const polyMesh&, + const sweepData&, + const scalar +) const +{ + return true; +} + + +inline void Foam::sweepData::leaveDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& faceCentre +) +{ + origin_ -= faceCentre; +} + + +inline void Foam::sweepData::transform +( + const polyMesh&, + const tensor& rotTensor +) +{ + origin_ = Foam::transform(rotTensor, origin_); +} + + +inline void Foam::sweepData::enterDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& faceCentre +) +{ + // back to absolute form + origin_ += faceCentre; +} + + +inline bool Foam::sweepData::updateCell +( + const polyMesh& mesh, + const label thisCellI, + const label, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.cellCentres()[thisCellI], tol); +} + + +inline bool Foam::sweepData::updateFace +( + const polyMesh& mesh, + const label thisFaceI, + const label, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.faceCentres()[thisFaceI], tol); +} + + +// Update this (face) with coupled face information. +inline bool Foam::sweepData::updateFace +( + const polyMesh& mesh, + const label thisFaceI, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.faceCentres()[thisFaceI], tol); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline void Foam::sweepData::operator= +( + const scalar value +) +{ + value_ = value; +} + + +inline bool Foam::sweepData::operator== +( + const sweepData& rhs +) const +{ + return origin() == rhs.origin(); +} + + +inline bool Foam::sweepData::operator!= +( + const sweepData& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* //