diff --git a/applications/solvers/coupled/pUCoupledFoam/couplingTerms.H b/applications/solvers/coupled/pUCoupledFoam/couplingTerms.H new file mode 100644 index 000000000..9aa6f9510 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/couplingTerms.H @@ -0,0 +1,15 @@ +{ + // Calculate grad p coupling matrix. Needs to be here if one uses + // gradient schemes with limiters. VV, 9/June/2014 + BlockLduSystem pInU(fvm::grad(p)); + + // Calculate div U coupling. Could be calculated only once since + // it is only geometry dependent. VV, 9/June/2014 + BlockLduSystem UInp(fvm::UDiv(U)); + + // Last argument in insertBlockCoupling says if the column direction + // should be incremented. This is needed for arbitrary positioning + // of U and p in the system. This could be better. VV, 30/April/2014 + UpEqn.insertBlockCoupling(0, 3, pInU, true); + UpEqn.insertBlockCoupling(3, 0, UInp, false); +} diff --git a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C index c0605549f..088d41759 100644 --- a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C +++ b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C @@ -38,9 +38,9 @@ Authors \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "fvBlockMatrix.H" #include "singlePhaseTransportModel.H" #include "RASModel.H" -#include "fvBlockMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -54,9 +54,6 @@ int main(int argc, char *argv[]) # include "initContinuityErrs.H" # include "readBlockSolverControls.H" - // Calculate div U coupling only once since it is only geometry dependant. - BlockLduSystem UInp(fvm::div(U)); - Info<< "\nStarting time loop\n" << endl; while (runTime.loop()) { @@ -73,15 +70,8 @@ int main(int argc, char *argv[]) // Assemble and insert pressure equation # include "pEqn.H" - // Calculate grad p coupling matrix. Needs to be here if one uses - // gradient schemes with limiters. VV, 9/June/2014 - BlockLduSystem pInU(fvm::grad(p)); - - // Last argument in insertBlockCoupling says if the column direction - // should be incremented. This is needed for arbitrary positioning - // of U and p in the system. This could be better. VV, 30/April/2014 - UpEqn.insertBlockCoupling(0, 3, pInU, true); - UpEqn.insertBlockCoupling(3, 0, UInp, false); + // Assemble and insert coupling terms +# include "couplingTerms.H" // Solve the block matrix UpEqn.solve(); diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files deleted file mode 100644 index 7bda4a174..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -pUCoupledFullPicardFoam.C - -EXE = $(FOAM_APPBIN)/pUCoupledFullPicardFoam diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options deleted file mode 100644 index 5d2b99e6d..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options +++ /dev/null @@ -1,12 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/turbulenceModels \ - -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/finiteVolume/lnInclude - -EXE_LIBS = \ - -lincompressibleRASModels \ - -lincompressibleTransportModels \ - -lfiniteVolume \ - -llduSolvers diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H b/applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H deleted file mode 100644 index aaf374e78..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H +++ /dev/null @@ -1,12 +0,0 @@ -// Momentum equation -fvVectorMatrix UEqn -( - fvm::div(phi, U) - + fvc::div(-phi, U, "div(phi,U)") - + turbulence->divDevReff(U) -); - -UEqn.relax(); - -UpEqn.insertEquation(0, UEqn); -UpEqn.insertPicardTensor(0, U, phi); diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H b/applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H deleted file mode 100644 index 1ea6ea880..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H +++ /dev/null @@ -1,52 +0,0 @@ -Info << "Reading field p\n" << endl; -volScalarField p -( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -Info << "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "createPhi.H" - -singlePhaseTransportModel laminarTransport(U, phi); -autoPtr turbulence -( - incompressible::RASModel::New(U, phi, laminarTransport) -); - -// Block vector field for velocity (first entry) and pressure (second -// entry). -Info << "Creating field Up\n" << endl; -volVector4Field Up -( - IOobject - ( - "Up", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedVector4("zero", dimless, vector4::zero) -); diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H b/applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H deleted file mode 100644 index 9029b873f..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H +++ /dev/null @@ -1,23 +0,0 @@ -// Pressure parts of the continuity equation -surfaceScalarField rUAf -( - "rUAf", - fvc::interpolate(1.0/UEqn.A()) -); - -surfaceScalarField presSource -( - "presSource", - rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf()) -); - -fvScalarMatrix pEqn -( - - fvm::laplacian(rUAf, p) - == - - fvc::div(presSource) -); - -pEqn.setReference(pRefCell, pRefValue); - -UpEqn.insertEquation(3, pEqn); diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C b/applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C deleted file mode 100644 index 17aa125a8..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C +++ /dev/null @@ -1,116 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Application - pUCoupledFullPicardFoam - -Description - Steady-state solver for incompressible, turbulent flow, with implicit - coupling between pressure and velocity achieved by fvBlockMatrix. - Turbulence is in this version solved using the existing turbulence - structure. - Convective term is properly linearised using Picard linearisation and - additional block coefficients. This is full, proper variant! - VV, 23/July/2014. - -Authors - Klas Jareteg, Chalmers University of Technology, - Vuko Vukcevic, FMENA Zagreb. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "singlePhaseTransportModel.H" -#include "RASModel.H" -#include "fvBlockMatrix.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "createFields.H" -# include "initContinuityErrs.H" -# include "readBlockSolverControls.H" - - // Calculate div U coupling only once since it is only geometry dependant. - BlockLduSystem UInp(fvm::div(U)); - - Info<< "\nStarting time loop\n" << endl; - while (runTime.loop()) - { - Info<< "Time = " << runTime.timeName() << nl << endl; - - p.storePrevIter(); - - // Initialize the Up block system (matrix, source and reference to Up) - fvBlockMatrix UpEqn(Up); - - // Assemble and insert momentum equation -# include "UEqn.H" - - // Assemble and insert pressure equation -# include "pEqn.H" - - // Calculate grad p coupling matrix. Needs to be here if one uses - // gradient schemes with limiters. VV, 9/June/2014 - BlockLduSystem pInU(fvm::grad(p)); - - // Last argument in insertBlockCoupling says if the column direction - // should be incremented. This is needed for arbitrary positioning - // of U and p in the system. This could be better. VV, 30/April/2014 - UpEqn.insertBlockCoupling(0, 3, pInU, true); - UpEqn.insertBlockCoupling(3, 0, UInp, false); - - // Solve the block matrix - UpEqn.solve(); - - // Retrieve solution - UpEqn.retrieveSolution(0, U.internalField()); - UpEqn.retrieveSolution(3, p.internalField()); - - U.correctBoundaryConditions(); - p.correctBoundaryConditions(); - - phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; - -# include "continuityErrs.H" - - p.relax(); - - 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/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H b/applications/solvers/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H deleted file mode 100644 index 4ee98e519..000000000 --- a/applications/solvers/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H +++ /dev/null @@ -1,3 +0,0 @@ -label pRefCell = 0; -scalar pRefValue = 0; -setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue); diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files b/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files deleted file mode 100644 index b7659c7d3..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -pUCoupledSemiPicardFoam.C - -EXE = $(FOAM_APPBIN)/pUCoupledSemiPicardFoam diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options b/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options deleted file mode 100644 index 5d2b99e6d..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options +++ /dev/null @@ -1,12 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/turbulenceModels \ - -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ - -I$(LIB_SRC)/transportModels \ - -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/finiteVolume/lnInclude - -EXE_LIBS = \ - -lincompressibleRASModels \ - -lincompressibleTransportModels \ - -lfiniteVolume \ - -llduSolvers diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H b/applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H deleted file mode 100644 index f0e2d9d75..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H +++ /dev/null @@ -1,11 +0,0 @@ -// Momentum equation -fvVectorMatrix UEqn -( - fvm::div(2*phi, U,"div(phi,U)") - + fvc::div(-phi, U,"div(phi,U)") - + turbulence->divDevReff(U) -); - -UEqn.relax(); - -UpEqn.insertEquation(0, UEqn); diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H b/applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H deleted file mode 100644 index 1ea6ea880..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H +++ /dev/null @@ -1,52 +0,0 @@ -Info << "Reading field p\n" << endl; -volScalarField p -( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -Info << "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "createPhi.H" - -singlePhaseTransportModel laminarTransport(U, phi); -autoPtr turbulence -( - incompressible::RASModel::New(U, phi, laminarTransport) -); - -// Block vector field for velocity (first entry) and pressure (second -// entry). -Info << "Creating field Up\n" << endl; -volVector4Field Up -( - IOobject - ( - "Up", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedVector4("zero", dimless, vector4::zero) -); diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H b/applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H deleted file mode 100644 index 9029b873f..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H +++ /dev/null @@ -1,23 +0,0 @@ -// Pressure parts of the continuity equation -surfaceScalarField rUAf -( - "rUAf", - fvc::interpolate(1.0/UEqn.A()) -); - -surfaceScalarField presSource -( - "presSource", - rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf()) -); - -fvScalarMatrix pEqn -( - - fvm::laplacian(rUAf, p) - == - - fvc::div(presSource) -); - -pEqn.setReference(pRefCell, pRefValue); - -UpEqn.insertEquation(3, pEqn); diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C b/applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C deleted file mode 100644 index 160ff488e..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C +++ /dev/null @@ -1,116 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Application - pUCoupledSemiPicardFoam - -Description - Steady-state solver for incompressible, turbulent flow, with implicit - coupling between pressure and velocity achieved by fvBlockMatrix. - Turbulence is in this version solved using the existing turbulence - structure. - Convective term is properly linearised using Picard linearisation and - additional block coefficients. This is semi variant where the difference - between succesive velocity fields is neglected! VV, 23/July/2014. - -Authors - Klas Jareteg, Chalmers University of Technology, - Vuko Vukcevic, FMENA Zagreb. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "singlePhaseTransportModel.H" -#include "RASModel.H" -#include "fvBlockMatrix.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - -# include "setRootCase.H" -# include "createTime.H" -# include "createMesh.H" -# include "createFields.H" -# include "initContinuityErrs.H" -# include "readBlockSolverControls.H" - - // Calculate div U coupling only once since it is only geometry dependant. - BlockLduSystem UInp(fvm::div(U)); - - Info<< "\nStarting time loop\n" << endl; - while (runTime.loop()) - { - Info<< "Time = " << runTime.timeName() << nl << endl; - - p.storePrevIter(); - - // Initialize the Up block system (matrix, source and reference to Up) - fvBlockMatrix UpEqn(Up); - - // Assemble and insert momentum equation -# include "UEqn.H" - - // Assemble and insert pressure equation -# include "pEqn.H" - - // Calculate grad p coupling matrix. Needs to be here if one uses - // gradient schemes with limiters. VV, 9/June/2014 - BlockLduSystem pInU(fvm::grad(p)); - - // Last argument in insertBlockCoupling says if the column direction - // should be incremented. This is needed for arbitrary positioning - // of U and p in the system. This could be better. VV, 30/April/2014 - UpEqn.insertBlockCoupling(0, 3, pInU, true); - UpEqn.insertBlockCoupling(3, 0, UInp, false); - - // Solve the block matrix - UpEqn.solve(); - - // Retrieve solution - UpEqn.retrieveSolution(0, U.internalField()); - UpEqn.retrieveSolution(3, p.internalField()); - - U.correctBoundaryConditions(); - p.correctBoundaryConditions(); - - phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; - -# include "continuityErrs.H" - - p.relax(); - - 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/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H b/applications/solvers/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H deleted file mode 100644 index 4ee98e519..000000000 --- a/applications/solvers/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H +++ /dev/null @@ -1,3 +0,0 @@ -label pRefCell = 0; -scalar pRefValue = 0; -setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue); diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C index 8e3ce3640..2f4d56c86 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C @@ -103,7 +103,7 @@ tmp < BlockLduSystem::type> > -divScheme::fvmDiv +divScheme::fvmUDiv ( const GeometricField& vf ) const @@ -128,6 +128,38 @@ divScheme::fvmDiv } +template +tmp +< + BlockLduSystem::type> +> +divScheme::fvmUDiv +( + const surfaceScalarField& flux, + const GeometricField& vf +) const +{ + FatalErrorIn + ( + "tmp divScheme::fvmDiv\n" + "(\n" + " surfaceScalarField&" + " GeometricField&" + ")\n" + ) << "Implicit div operator currently defined only for Gauss linear. " + << abort(FatalError); + + typedef typename innerProduct::type DivType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H index 06705323d..d22159bcc 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H @@ -157,10 +157,19 @@ public: virtual tmp < BlockLduSystem::type> - > fvmDiv + > fvmUDiv ( const GeometricField& ) const; + + virtual tmp + < + BlockLduSystem::type> + > fvmUDiv + ( + const surfaceScalarField&, + const GeometricField& + ) const; }; diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C index c2888a554..1815412b9 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C @@ -72,14 +72,14 @@ template tmp < BlockLduSystem::type> -> gaussDivScheme::fvmDiv +> gaussDivScheme::fvmUDiv ( const GeometricField& vf ) const { FatalErrorIn ( - "tmp fvmDiv\n" + "tmp fvmUDiv\n" "(\n" " GeometricField&" ")\n" @@ -97,6 +97,37 @@ tmp } +template +tmp +< + BlockLduSystem::type> +> gaussDivScheme::fvmUDiv +( + const surfaceScalarField& flux, + const GeometricField& vf +) const +{ + FatalErrorIn + ( + "tmp fvmUDiv\n" + "(\n" + " const surfaceScalarField& flux" + " const GeometricField&" + ")\n" + ) << "Implicit div operator defined only for vector." + << abort(FatalError); + + typedef typename innerProduct::type DivType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H index bba569e7c..68e812a10 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H @@ -102,10 +102,19 @@ public: tmp < BlockLduSystem::type> - > fvmDiv + > fvmUDiv ( const GeometricField& ) const; + + tmp + < + BlockLduSystem::type> + > fvmUDiv + ( + const surfaceScalarField& flux, + const GeometricField& + ) const; }; diff --git a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C index 09312b4cf..93250f398 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C @@ -42,7 +42,7 @@ namespace fv // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template<> -tmp > gaussDivScheme::fvmDiv +tmp > gaussDivScheme::fvmUDiv ( const GeometricField& vf ) const @@ -60,9 +60,9 @@ tmp > gaussDivScheme::fvmDiv scalarField& source = bs.source(); // Grab ldu parts of block matrix as linear always - typename CoeffField::linearTypeField& d = bs.diag().asLinear(); - typename CoeffField::linearTypeField& u = bs.upper().asLinear(); - typename CoeffField::linearTypeField& l = bs.lower().asLinear(); + CoeffField::linearTypeField& d = bs.diag().asLinear(); + CoeffField::linearTypeField& u = bs.upper().asLinear(); + CoeffField::linearTypeField& l = bs.lower().asLinear(); const vectorField& SfIn = mesh.Sf().internalField(); @@ -89,9 +89,9 @@ tmp > gaussDivScheme::fvmDiv if (patch.coupled()) { - typename CoeffField::linearTypeField& pcoupleUpper = + CoeffField::linearTypeField& pcoupleUpper = bs.coupleUpper()[patchI].asLinear(); - typename CoeffField::linearTypeField& pcoupleLower = + CoeffField::linearTypeField& pcoupleLower = bs.coupleLower()[patchI].asLinear(); const vectorField pcl = -pw*Sf; @@ -119,6 +119,95 @@ tmp > gaussDivScheme::fvmDiv } +template<> +tmp > gaussDivScheme::fvmUDiv +( + const surfaceScalarField& flux, + const GeometricField& vf +) const +{ + tmp tweights = this->tinterpScheme_().weights(vf); + const scalarField& wIn = tweights().internalField(); + + const fvMesh& mesh = vf.mesh(); + + tmp > tbs + ( + new BlockLduSystem(mesh) + ); + BlockLduSystem& bs = tbs(); + scalarField& source = bs.source(); + + // Grab ldu parts of block matrix as linear always + CoeffField::linearTypeField& d = bs.diag().asLinear(); + CoeffField::linearTypeField& u = bs.upper().asLinear(); + CoeffField::linearTypeField& l = bs.lower().asLinear(); + + const vectorField& SfIn = mesh.Sf().internalField(); + + const scalarField& fluxIn = flux.internalField(); + + l = -wIn*fluxIn*SfIn; + u = l + fluxIn*SfIn; + bs.negSumDiag(); + + // Boundary contributions + forAll(vf.boundaryField(), patchI) + { + const fvPatchVectorField& pf = vf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + const vectorField& Sf = patch.Sf(); + const fvsPatchScalarField& pw = tweights().boundaryField()[patchI]; + const unallocLabelList& fc = patch.faceCells(); + + const scalarField& pFlux = flux.boundaryField()[patchI]; + + const vectorField internalCoeffs(pf.valueInternalCoeffs(pw)); + + // Diag contribution + forAll(pf, faceI) + { + d[fc[faceI]] += cmptMultiply + ( + internalCoeffs[faceI], + pFlux[faceI]*Sf[faceI] + ); + } + + if (patch.coupled()) + { + CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const vectorField pcl = -pw*pFlux*Sf; + const vectorField pcu = pcl + pFlux*Sf; + + // Coupling contributions + pcoupleLower -= pcl; + pcoupleUpper -= pcu; + } + else + { + const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); + + // Boundary contribution + forAll(pf, faceI) + { + source[fc[faceI]] -= + ( + boundaryCoeffs[faceI] & (pFlux[faceI]*Sf[faceI]) + ); + } + } + } + + // Interpolation schemes with corrections not accounted for + + return tbs; +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H index a09f30589..0fe0f58e9 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H @@ -51,11 +51,18 @@ namespace fv // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template<> -tmp > gaussDivScheme::fvmDiv +tmp > gaussDivScheme::fvmUDiv ( const GeometricField& ) const; +template<> +tmp > gaussDivScheme::fvmUDiv +( + const surfaceScalarField& flux, + const GeometricField& +) const; + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C index 25a9028d1..e57a61af3 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C @@ -103,7 +103,7 @@ template tmp < BlockLduSystem::type> -> div +> UDiv ( GeometricField& vf, const word& name @@ -113,7 +113,7 @@ tmp ( vf.mesh(), vf.mesh().schemesDict().divScheme(name) - )().fvmDiv(vf); + )().fvmUDiv(vf); } @@ -121,12 +121,52 @@ template tmp < BlockLduSystem::type> -> div +> UDiv +( + const surfaceScalarField& flux, + GeometricField& vf, + const word& name +) +{ + return fv::divScheme::New + ( + vf.mesh(), + vf.mesh().schemesDict().divScheme(name) + )().fvmUDiv(flux, vf); +} + + +template +tmp +< + BlockLduSystem::type> +> UDiv +( + const tmp& tflux, + GeometricField& vf, + const word& name +) +{ + tmp + < + BlockLduSystem::type> + > + Div(fvm::UDiv(tflux(), vf, name)); + tflux.clear(); + return Div; +} + + +template +tmp +< + BlockLduSystem::type> +> UDiv ( GeometricField& vf ) { - return fvm::div + return fvm::UDiv ( vf, "div(" + vf.name() + ')' @@ -134,6 +174,45 @@ tmp } +template +tmp +< + BlockLduSystem::type> +> UDiv +( + const surfaceScalarField& flux, + GeometricField& vf +) +{ + return fvm::UDiv + ( + flux, + vf, + "div(" + vf.name() + ')' + ); +} + + +template +tmp +< + BlockLduSystem::type> +> UDiv +( + const tmp& tflux, + GeometricField& vf +) +{ + tmp + < + BlockLduSystem::type> + > + Div(fvm::UDiv(tflux(), vf)); + tflux.clear(); + return Div; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fvm diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H index 465bae1d1..a386a83db 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H @@ -68,7 +68,6 @@ namespace fvm const word& name ); - template tmp > div ( @@ -83,12 +82,12 @@ namespace fvm GeometricField& ); - // Implicit div operator for block systems + // Implicit div operators for block systems template tmp < BlockLduSystem::type> - > div + > UDiv ( GeometricField&, const word& @@ -98,8 +97,50 @@ namespace fvm tmp < BlockLduSystem::type> - > div + > UDiv ( + const surfaceScalarField&, + GeometricField&, + const word& + ); + + template + tmp + < + BlockLduSystem::type> + > UDiv + ( + const tmp&, + GeometricField&, + const word& + ); + + template + tmp + < + BlockLduSystem::type> + > UDiv + ( + GeometricField& + ); + + template + tmp + < + BlockLduSystem::type> + > UDiv + ( + const surfaceScalarField&, + GeometricField& + ); + + template + tmp + < + BlockLduSystem::type> + > UDiv + ( + const tmp&, GeometricField& ); } diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C index 801185958..40d3e1e5d 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C @@ -60,9 +60,9 @@ tmp > gaussGrad::fvmGrad vectorField& source = bs.source(); // Grab ldu parts of block matrix as linear always - typename CoeffField::linearTypeField& d = bs.diag().asLinear(); - typename CoeffField::linearTypeField& u = bs.upper().asLinear(); - typename CoeffField::linearTypeField& l = bs.lower().asLinear(); + CoeffField::linearTypeField& d = bs.diag().asLinear(); + CoeffField::linearTypeField& u = bs.upper().asLinear(); + CoeffField::linearTypeField& l = bs.lower().asLinear(); const vectorField& SfIn = mesh.Sf().internalField(); @@ -89,9 +89,9 @@ tmp > gaussGrad::fvmGrad if (patch.coupled()) { - typename CoeffField::linearTypeField& pcoupleUpper = + CoeffField::linearTypeField& pcoupleUpper = bs.coupleUpper()[patchI].asLinear(); - typename CoeffField::linearTypeField& pcoupleLower = + CoeffField::linearTypeField& pcoupleLower = bs.coupleLower()[patchI].asLinear(); const vectorField pcl = -pw*Sf; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C index 5f74e6240..1b3691312 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C @@ -80,9 +80,9 @@ tmp > leastSquaresGrad::fvmGrad vectorField& source = bs.source(); // Grab ldu parts of block matrix as linear always - typename CoeffField::linearTypeField& d = bs.diag().asLinear(); - typename CoeffField::linearTypeField& u = bs.upper().asLinear(); - typename CoeffField::linearTypeField& l = bs.lower().asLinear(); + CoeffField::linearTypeField& d = bs.diag().asLinear(); + CoeffField::linearTypeField& u = bs.upper().asLinear(); + CoeffField::linearTypeField& l = bs.lower().asLinear(); // Get references to least square vectors const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh); @@ -125,9 +125,9 @@ tmp > leastSquaresGrad::fvmGrad const scalarField cellVInNei = cellV.boundaryField()[patchI].patchNeighbourField(); - typename CoeffField::linearTypeField& pcoupleUpper = + CoeffField::linearTypeField& pcoupleUpper = bs.coupleUpper()[patchI].asLinear(); - typename CoeffField::linearTypeField& pcoupleLower = + CoeffField::linearTypeField& pcoupleLower = bs.coupleLower()[patchI].asLinear(); // Coupling and diagonal contributions