diff --git a/applications/solvers/coupled/Allwmake b/applications/solvers/coupled/Allwmake index 82e54df8e..9ce03ddae 100755 --- a/applications/solvers/coupled/Allwmake +++ b/applications/solvers/coupled/Allwmake @@ -8,4 +8,6 @@ wmake blockCoupledScalarTransportFoam wmake conjugateHeatFoam wmake conjugateHeatSimpleFoam wmake pUCoupledFoam +wmake pUCoupledFullPicardFoam +wmake pUCoupledSemiPicardFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/Make/options b/applications/solvers/coupled/blockCoupledScalarTransportFoam/Make/options index 068c46323..436ff1d82 100644 --- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/Make/options +++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/Make/options @@ -4,5 +4,4 @@ EXE_INC = \ -I$(LIB_SRC)/VectorN/lnInclude EXE_LIBS = \ - -lfiniteVolume \ - -lVectorN + -lfiniteVolume diff --git a/applications/solvers/coupled/conjugateHeatTransfer/Make/options b/applications/solvers/coupled/conjugateHeatTransfer/Make/options index 721166578..f93f4f661 100644 --- a/applications/solvers/coupled/conjugateHeatTransfer/Make/options +++ b/applications/solvers/coupled/conjugateHeatTransfer/Make/options @@ -8,6 +8,4 @@ LIB_LIBS = \ -lfiniteVolume \ -lspecie \ -lbasicThermophysicalModels \ - -lradiation \ - -lVectorN - + -lradiation diff --git a/applications/solvers/coupled/pUCoupledFoam/Make/options b/applications/solvers/coupled/pUCoupledFoam/Make/options index 5775a4d11..5d2b99e6d 100644 --- a/applications/solvers/coupled/pUCoupledFoam/Make/options +++ b/applications/solvers/coupled/pUCoupledFoam/Make/options @@ -3,13 +3,10 @@ EXE_INC = \ -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/VectorN/lnInclude + -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -lincompressibleRASModels \ -lincompressibleTransportModels \ -lfiniteVolume \ - -llduSolvers \ - -lVectorN - + -llduSolvers diff --git a/applications/solvers/coupled/pUCoupledFoam/UEqn.H b/applications/solvers/coupled/pUCoupledFoam/UEqn.H index 31693a9ab..7219eb39e 100644 --- a/applications/solvers/coupled/pUCoupledFoam/UEqn.H +++ b/applications/solvers/coupled/pUCoupledFoam/UEqn.H @@ -7,4 +7,4 @@ fvVectorMatrix UEqn UEqn.relax(); -blockMatrixTools::insertEquation(0, UEqn, A, x, b); +UpEqn.insertEquation(0, UEqn); diff --git a/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H b/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H deleted file mode 100644 index a89c6bf2a..000000000 --- a/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H +++ /dev/null @@ -1,2 +0,0 @@ -blockVectorMatrix pInU(fvm::grad(p)); -blockVectorMatrix UInp(fvm::div(U)); diff --git a/applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H b/applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H deleted file mode 100644 index 2e7d61e61..000000000 --- a/applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H +++ /dev/null @@ -1,9 +0,0 @@ -// Create block matrix -BlockLduMatrix A(mesh); - -// Block matrix - create x and b -vector4Field& x = Up.internalField(); -vector4Field b(x.size(), vector4::zero); - -// Set block interfaces properly -A.interfaces() = Up.boundaryField().blockInterfaces(); diff --git a/applications/solvers/coupled/pUCoupledFoam/pEqn.H b/applications/solvers/coupled/pUCoupledFoam/pEqn.H index 03a32161f..9029b873f 100644 --- a/applications/solvers/coupled/pUCoupledFoam/pEqn.H +++ b/applications/solvers/coupled/pUCoupledFoam/pEqn.H @@ -8,7 +8,7 @@ surfaceScalarField rUAf surfaceScalarField presSource ( "presSource", - rUAf*(fvc::interpolate(fvc::grad(p)) & mesh.Sf()) + rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf()) ); fvScalarMatrix pEqn @@ -20,4 +20,4 @@ fvScalarMatrix pEqn pEqn.setReference(pRefCell, pRefValue); -blockMatrixTools::insertEquation(3, pEqn, A, x, b); +UpEqn.insertEquation(3, pEqn); diff --git a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C index 3c30b5823..c0605549f 100644 --- a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C +++ b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C @@ -1,32 +1,33 @@ /*---------------------------------------------------------------------------*\ ========= | - \\ / F ield | foam-extend: Open Source CFD + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | For copyright notice see file Copyright + \\ / A nd | Copyright held by original author \\/ M anipulation | ------------------------------------------------------------------------------- License - This file is part of foam-extend. + This file is part of OpenFOAM. - foam-extend is free software: you can redistribute it and/or modify it + 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 + Free Software Foundation; either version 2 of the License, or (at your option) any later version. - foam-extend 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. + 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 foam-extend. If not, see . + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Application pUCoupledFoam Description Steady-state solver for incompressible, turbulent flow, with implicit - coupling between pressure and velocity achieved by BlockLduMatrix + coupling between pressure and velocity achieved by fvBlockMatrix. Turbulence is in this version solved using the existing turbulence structure. @@ -39,13 +40,7 @@ Authors #include "fvCFD.H" #include "singlePhaseTransportModel.H" #include "RASModel.H" - -#include "VectorNFieldTypes.H" -#include "volVectorNFields.H" -#include "blockLduSolvers.H" -#include "blockVectorNMatrices.H" - -#include "blockMatrixTools.H" +#include "fvBlockMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -59,12 +54,8 @@ int main(int argc, char *argv[]) # include "initContinuityErrs.H" # include "readBlockSolverControls.H" - // Calculate coupling matrices only once since the mesh doesn't change and - // implicit div and grad operators are only dependant on Sf. Actually - // coupling terms (div(U) and grad(p)) in blockMatrix do not change, - // so they could be inserted only once, resetting other parts of - // blockMatrix to zero at the end of each time step. VV, 30/April/2014 -# include "calculateCouplingMatrices.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,8 +64,8 @@ int main(int argc, char *argv[]) p.storePrevIter(); - // Initialize block matrix -# include "initializeBlockMatrix.H" + // Initialize the Up block system (matrix, source and reference to Up) + fvBlockMatrix UpEqn(Up); // Assemble and insert momentum equation # include "UEqn.H" @@ -82,27 +73,22 @@ int main(int argc, char *argv[]) // Assemble and insert pressure equation # include "pEqn.H" - // Insert coupling, updating the boundary contributions - // Last argument in insertBlockCoupling says if the first location + // 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 - blockMatrixTools::insertBlockCoupling(3, 0, UInp, U, A, b, false); - blockMatrixTools::insertBlockCoupling(0, 3, pInU, p, A, b, true); + UpEqn.insertBlockCoupling(0, 3, pInU, true); + UpEqn.insertBlockCoupling(3, 0, UInp, false); // Solve the block matrix - BlockSolverPerformance solverPerf = - BlockLduSolver::New - ( - word("Up"), - A, - mesh.solutionDict().solver("Up") - )->solve(Up, b); - - solverPerf.print(); + UpEqn.solve(); // Retrieve solution - blockMatrixTools::retrieveSolution(0, U.internalField(), Up); - blockMatrixTools::retrieveSolution(3, p.internalField(), Up); + UpEqn.retrieveSolution(0, U.internalField()); + UpEqn.retrieveSolution(3, p.internalField()); U.correctBoundaryConditions(); p.correctBoundaryConditions(); diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files new file mode 100644 index 000000000..7bda4a174 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files @@ -0,0 +1,3 @@ +pUCoupledFullPicardFoam.C + +EXE = $(FOAM_APPBIN)/pUCoupledFullPicardFoam diff --git a/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options new file mode 100644 index 000000000..5d2b99e6d --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options @@ -0,0 +1,12 @@ +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 new file mode 100644 index 000000000..aaf374e78 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H @@ -0,0 +1,12 @@ +// 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 new file mode 100644 index 000000000..1ea6ea880 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H @@ -0,0 +1,52 @@ +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 new file mode 100644 index 000000000..9029b873f --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H @@ -0,0 +1,23 @@ +// 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 new file mode 100644 index 000000000..17aa125a8 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 new file mode 100644 index 000000000..4ee98e519 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H @@ -0,0 +1,3 @@ +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 new file mode 100644 index 000000000..b7659c7d3 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files @@ -0,0 +1,3 @@ +pUCoupledSemiPicardFoam.C + +EXE = $(FOAM_APPBIN)/pUCoupledSemiPicardFoam diff --git a/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options b/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options new file mode 100644 index 000000000..5d2b99e6d --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options @@ -0,0 +1,12 @@ +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 new file mode 100644 index 000000000..f0e2d9d75 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H @@ -0,0 +1,11 @@ +// 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 new file mode 100644 index 000000000..1ea6ea880 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H @@ -0,0 +1,52 @@ +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 new file mode 100644 index 000000000..9029b873f --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H @@ -0,0 +1,23 @@ +// 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 new file mode 100644 index 000000000..160ff488e --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 new file mode 100644 index 000000000..4ee98e519 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H @@ -0,0 +1,3 @@ +label pRefCell = 0; +scalar pRefValue = 0; +setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue); diff --git a/src/Allwmake b/src/Allwmake index 59c69dd20..cdf614c4d 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -37,7 +37,6 @@ wmake libso finiteVolume wmake libso finiteArea wmake libso lduSolvers -wmake libso VectorN wmake libso tetFiniteElement diff --git a/src/VectorN/Make/files b/src/VectorN/Make/files deleted file mode 100644 index 26b8df629..000000000 --- a/src/VectorN/Make/files +++ /dev/null @@ -1,28 +0,0 @@ -foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C - -finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.C -finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.C - -finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.C -finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.C -finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.C -finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.C -finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.C -finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.C - -finiteVolume/fields/volFields/volVectorNFields.C -finiteVolume/fields/surfaceFields/surfaceVectorNFields.C - -finiteVolume/interpolation/VectorNSurfaceInterpolationSchemes.C - -LIB = $(FOAM_LIBBIN)/libVectorN diff --git a/src/VectorN/Make/options b/src/VectorN/Make/options deleted file mode 100644 index 958e50b9e..000000000 --- a/src/VectorN/Make/options +++ /dev/null @@ -1,4 +0,0 @@ -EXE_INC = \ - -I$(FOAM_SRC)/finiteVolume/lnInclude - -LIB_LIBS = diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 34759660a..941e623cf 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -172,7 +172,19 @@ $(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrosta $(derivedFvPatchFields)/pulseFixedValue/pulseFixedValueFvPatchFields.C $(derivedFvPatchFields)/waveTransmissiveInlet/waveTransmissiveInletFvPatchFields.C - +fvPatchVectorNFields = $(fvPatchFields)/fvPatchVectorNFields +$(fvPatchVectorNFields)/fvPatchVectorNFields.C +$(fvPatchVectorNFields)/genericFvPatchVectorNFields.C +$(fvPatchVectorNFields)/calculatedFvPatchVectorNFields.C +$(fvPatchVectorNFields)/emptyFvPatchVectorNFields.C +$(fvPatchVectorNFields)/transformFvPatchVectorNFields.C +$(fvPatchVectorNFields)/wedgeFvPatchVectorNFields.C +$(fvPatchVectorNFields)/coupledFvPatchVectorNFields.C +$(fvPatchVectorNFields)/processorFvPatchVectorNFields.C +$(fvPatchVectorNFields)/cyclicFvPatchVectorNFields.C +$(fvPatchVectorNFields)/fixedValueFvPatchVectorNFields.C +$(fvPatchVectorNFields)/zeroGradientFvPatchVectorNFields.C +$(fvPatchVectorNFields)/fixedGradientFvPatchVectorNFields.C fvsPatchFields = fields/fvsPatchFields $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C @@ -195,8 +207,18 @@ $(constraintFvsPatchFields)/overlapGgi/overlapGgiFvsPatchFields.C $(constraintFvsPatchFields)/mixingPlane/mixingPlaneFvsPatchFields.C $(constraintFvsPatchFields)/regionCoupling/regionCouplingFvsPatchFields.C +fvsPatchVectorNFields = $(fvsPatchFields)/fvsPatchVectorNFields +$(fvsPatchVectorNFields)/fvsPatchVectorNFields.C +$(fvsPatchVectorNFields)/calculatedFvsPatchVectorNFields.C +$(fvsPatchVectorNFields)/emptyFvsPatchVectorNFields.C +$(fvsPatchVectorNFields)/wedgeFvsPatchVectorNFields.C +$(fvsPatchVectorNFields)/coupledFvsPatchVectorNFields.C +$(fvsPatchVectorNFields)/processorFvsPatchVectorNFields.C + fields/volFields/volFields.C +fields/volFields/volVectorNFields.C fields/surfaceFields/surfaceFields.C +fields/surfaceFields/surfaceVectorNFields.C fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C @@ -224,6 +246,7 @@ $(pointVolInterpolation)/pointVolInterpolation.C surfaceInterpolation = interpolation/surfaceInterpolation $(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C $(surfaceInterpolation)/surfaceInterpolationScheme/surfaceInterpolationSchemes.C +$(surfaceInterpolation)/VectorNSurfaceInterpolationSchemes/VectorNSurfaceInterpolationSchemes.C schemes = $(surfaceInterpolation)/schemes $(schemes)/linear/linear.C @@ -315,14 +338,17 @@ $(d2dt2Schemes)/EulerD2dt2Scheme/EulerD2dt2Schemes.C $(d2dt2Schemes)/backwardD2dt2Scheme/backwardD2dt2Schemes.C divSchemes = finiteVolume/divSchemes -$(divSchemes)/divScheme/divSchemes.C +$(divSchemes)/divScheme/divSchemes.Ci +$(divSchemes)/vectorGaussDivScheme/vectorGaussDivScheme.C $(divSchemes)/gaussDivScheme/gaussDivSchemes.C gradSchemes = finiteVolume/gradSchemes $(gradSchemes)/gradScheme/gradSchemes.C +$(gradSchemes)/scalarGaussGrad/scalarGaussGrad.C $(gradSchemes)/gaussGrad/gaussGrads.C $(gradSchemes)/beGaussGrad/beGaussGrads.C $(gradSchemes)/leastSquaresGrad/leastSquaresVectors.C +$(gradSchemes)/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C $(gradSchemes)/leastSquaresGrad/leastSquaresGrads.C $(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C $(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresGrads.C diff --git a/src/finiteVolume/blockLduSystem/BlockLduSystem.C b/src/finiteVolume/blockLduSystem/BlockLduSystem.C new file mode 100644 index 000000000..e11f08196 --- /dev/null +++ b/src/finiteVolume/blockLduSystem/BlockLduSystem.C @@ -0,0 +1,119 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "BlockLduSystem.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BlockLduSystem::BlockLduSystem +( + const lduMesh& ldu +) +: + BlockLduMatrix(ldu), + source_(ldu.lduAddr().size(), pTraits::zero) +{} + + +template +Foam::BlockLduSystem::BlockLduSystem +( + const lduMesh& ldu, + const Field& s +) +: + BlockLduMatrix(ldu), + source_(s) +{ + if (ldu.lduAddr().size() != s.size()) + { + FatalErrorIn + ( + "BlockLduSystem::BlockLduSystem\n" + "(\n" + " const lduMesh& ldu," + " const Field& s," + ")\n" + ) << "Sizes of ldu addressing and source field are not the same." + << abort(FatalError); + } +} + + +template +Foam::BlockLduSystem::BlockLduSystem +( + const BlockLduMatrix& bm, + const Field& s +) +: + BlockLduMatrix(bm), + source_(s) +{ + if (this->lduAddr().size() != s.size()) + { + FatalErrorIn + ( + "BlockLduSystem::BlockLduSystem\n" + "(\n" + " const BlockLduMatrix& bm," + " const Field& s," + ")\n" + ) << "Sizes of block matrix and source field are not the same." + << abort(FatalError); + } +} + + +template +Foam::BlockLduSystem::BlockLduSystem +( + const BlockLduSystem& bs +) +: + BlockLduMatrix(bs), + source_(bs.source()) +{} + + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +template +Foam::Ostream& Foam::operator<< +( + Ostream& os, + const BlockLduSystem& bs +) +{ + os << static_cast&>(bs) << nl + << bs.source() << endl; + + return os; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/blockLduSystem/BlockLduSystem.H b/src/finiteVolume/blockLduSystem/BlockLduSystem.H new file mode 100644 index 000000000..9ce89e8d2 --- /dev/null +++ b/src/finiteVolume/blockLduSystem/BlockLduSystem.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + Foam::BlockLduSystem + +Description + BlockLduSystem is a wrapper for BlockLduMatrix with source field. This is + the return type of implicit div and grad operators needed for implicitly + coupled systems (namely p-U coupling). + +Author + Vuko Vukcevic, FMENA Zagreb. + +SourceFiles + BlockLduSystem.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BlockLduSystem_H +#define BlockLduSystem_H + +#include "blockVectorNMatrices.H" +#include "VectorNFieldTypes.H" +#include "volVectorNFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of friend functions and operators + +template +class BlockLduSystem; + +template +Ostream& operator<<(Ostream&, const BlockLduSystem&); + + +/*---------------------------------------------------------------------------*\ + Class BlockLduSystem Declaration +\*---------------------------------------------------------------------------*/ + +template +class BlockLduSystem +: + public BlockLduMatrix +{ + // Private data + + //- Source term + Field source_; + + +public: + + // Constructors + + //- Construct given addressing + explicit BlockLduSystem(const lduMesh&); + + //- Construct given addressing and source field + BlockLduSystem(const lduMesh&, const Field&); + + //- Construct from components + BlockLduSystem + ( + const BlockLduMatrix&, + const Field& + ); + + //- Construct as copy + BlockLduSystem(const BlockLduSystem&); + + + // Destructor + + virtual ~BlockLduSystem() + {} + + + // Member functions + + //- Access + + Field& source() + { + return source_; + } + + const Field& source() const + { + return source_; + } + + + // Member operators + + void operator=(const BlockLduSystem&); + + void negate(); + + void operator+=(const BlockLduSystem&); + void operator-=(const BlockLduSystem&); + + void operator*=(const scalarField&); + void operator*=(const scalar); + + + // Ostream operator + + friend Ostream& operator<< + ( + Ostream&, + const BlockLduSystem& + ); +}; + + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "BlockLduSystem.C" +# include "BlockLduSystemOperations.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/blockLduSystem/BlockLduSystemOperations.C b/src/finiteVolume/blockLduSystem/BlockLduSystemOperations.C new file mode 100644 index 000000000..e87261174 --- /dev/null +++ b/src/finiteVolume/blockLduSystem/BlockLduSystemOperations.C @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "BlockLduSystem.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::BlockLduSystem::negate() +{ + BlockLduMatrix::negate(); + source_.negate(); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +void Foam::BlockLduSystem::operator= +( + const BlockLduSystem& bs +) +{ + if (this == &bs) + { + FatalErrorIn + ( + "void BlockLduSystem::operator=" + "(const BlockLduSystem& bs)" + ) << "attempted assignment to self" + << abort(FatalError); + } + + BlockLduMatrix::operator=(bs); + source_ = bs.source(); +} + + +template +void Foam::BlockLduSystem::operator+= +( + const BlockLduSystem& bs +) +{ + BlockLduMatrix::operator+=(bs); + source_ += bs.source(); +} + +template +void Foam::BlockLduSystem::operator-= +( + const BlockLduSystem& bs +) +{ + BlockLduMatrix::operator-=(bs); + source_ -= bs.source(); +} + +template +void Foam::BlockLduSystem::operator*= +( + const scalarField& sf +) +{ + BlockLduMatrix::operator*=(sf); + source_ *= sf; +} + +template +void Foam::BlockLduSystem::operator*= +( + const scalar s +) +{ + BlockLduMatrix::operator*=(s); + source_ *= s; +} + + +// ************************************************************************* // diff --git a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C b/src/finiteVolume/blockMatrixTools/blockMatrixTools.C similarity index 78% rename from src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C rename to src/finiteVolume/blockMatrixTools/blockMatrixTools.C index 0d0cc545a..5e1e6957f 100644 --- a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C +++ b/src/finiteVolume/blockMatrixTools/blockMatrixTools.C @@ -862,12 +862,12 @@ void insertEquation } -template +template void insertBlock ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, + const BlockLduSystem& blockSystem, BlockLduMatrix& A, const bool incFirst ) @@ -875,9 +875,9 @@ void insertBlock // Sanity checks { const label blockMatrixSize = pTraits::nComponents; - const label blockSystemSize = pTraits::nComponents; + const label blockMatrixASize = pTraits::nComponents; - if (blockMatrixSize > blockSystemSize) + if (blockMatrixSize > blockMatrixASize) { FatalErrorIn ( @@ -908,7 +908,6 @@ void insertBlock << "member function." << abort(FatalError); } - } const label nCmpts = pTraits::nComponents; @@ -917,11 +916,11 @@ void insertBlock // Get references to ldu fields of blockMatrix always as linear const typename CoeffField::linearTypeField& bmd = - blockMatrix.diag().asLinear(); + blockSystem.diag().asLinear(); const typename CoeffField::linearTypeField& bmu = - blockMatrix.upper().asLinear(); + blockSystem.upper().asLinear(); const typename CoeffField::linearTypeField& bml = - blockMatrix.lower().asLinear(); + blockSystem.lower().asLinear(); // Get references to ldu fields of A matrix always as square typename CoeffField::squareTypeField& blockDiag = @@ -936,13 +935,16 @@ void insertBlock { forAll(bmd, cellI) { - blockDiag[cellI](localLoc1, localLoc2) += bmd[cellI](cmptI); + blockDiag[cellI](localLoc1, localLoc2) += + bmd[cellI].component(cmptI); } forAll(bmu, faceI) { - blockUpper[faceI](localLoc1, localLoc2) += bmu[faceI](cmptI); - blockLower[faceI](localLoc1, localLoc2) += bml[faceI](cmptI); + blockUpper[faceI](localLoc1, localLoc2) += + bmu[faceI].component(cmptI); + blockLower[faceI](localLoc1, localLoc2) += + bml[faceI].component(cmptI); } if (incFirst) @@ -957,256 +959,86 @@ void insertBlock } -template +template void insertBoundaryContributions ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, + const BlockLduSystem& blockSystem, BlockLduMatrix& A, Field& b, const bool incFirst ) { - // Interpolation scheme for pressure gradient - tmp > - tinterpScheme - ( - surfaceInterpolationScheme::New - ( - psi.mesh(), - psi.mesh().schemesDict().interpolationScheme - ( - "grad(" + psi.name() + ')' - ) - ) - ); - - surfaceScalarField weights(tinterpScheme().weights(psi)); + // Need to get reference to fvMesh instead of lduMesh + const fvBoundaryMesh& bmesh = refCast(A.mesh()).boundary(); const label nCmpts = pTraits::nComponents; + const label nSrcCmpts = pTraits::nComponents; label localLoc1 = loc1; label localLoc2 = loc2; - // Get references to ldu fields of A matrix always as square - typename CoeffField::squareTypeField& blockDiag = - A.diag().asSquare(); + const Field& source = blockSystem.source(); - // Insert boundary contributions directly into the matrix, this assumes - // that blockMatrix came from implicit grad or div operators! - psi.boundaryField().updateCoeffs(); - forAll(psi.boundaryField(), patchI) + // Insert source from block system to rhs + for (label cmptI = 0; cmptI < nSrcCmpts; cmptI++) { - const fvPatchScalarField& pf = psi.boundaryField()[patchI]; - const fvPatch& patch = pf.patch(); - const vectorField& Sf = patch.Sf(); - const fvsPatchScalarField& pw = weights.boundaryField()[patchI]; + scalarField sourceCmpt(source.component(cmptI)); - scalarField internalCoeffs(pf.valueInternalCoeffs(pw)); - - if (patch.coupled()) + forAll(b, cellI) { - typename CoeffField::squareTypeField& pcoupleUpper = - A.coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pcoupleLower = - A.coupleLower()[patchI].asSquare(); + b[cellI](localLoc1) += sourceCmpt[cellI]; + } - vectorField pcl = -pw*Sf; - vectorField pcu = pcl + Sf; - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(pf, faceI) - { - label cellI = patch.faceCells()[faceI]; - - // Diag contribution - blockDiag[cellI](localLoc1, localLoc2) += - internalCoeffs[faceI]*Sf[faceI].component(cmptI); - - // Coupling contributions - pcoupleUpper[faceI](localLoc1, localLoc2) -= - pcu[faceI].component(cmptI); - pcoupleLower[faceI](localLoc1, localLoc2) -= - pcl[faceI].component(cmptI); - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - - // Reset local locations for other patches - localLoc1 = loc1; - localLoc2 = loc2; + if (incFirst) + { + localLoc1++; } else { - scalarField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(pf, faceI) - { - label cellI = patch.faceCells()[faceI]; - - // Diag contribution - blockDiag[cellI](localLoc1, localLoc2) += - internalCoeffs[faceI]*Sf[faceI].component(cmptI); - - // Source contribution - b[cellI](localLoc1) -= - boundaryCoeffs[faceI]*Sf[faceI].component(cmptI); - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - - // Reset local locations for other patches - localLoc1 = loc1; - localLoc2 = loc2; + localLoc2++; } } -} + // Reset local locations for coupling contributions + localLoc1 = loc1; + localLoc2 = loc2; -template -void insertBoundaryContributions -( - const label loc1, - const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, - BlockLduMatrix& A, - Field& b, - const bool incFirst -) -{ - // Interpolation scheme for velocity divergence - tmp > - tinterpScheme - ( - surfaceInterpolationScheme::New - ( - psi.mesh(), - psi.mesh().schemesDict().interpolationScheme - ( - "div(" + psi.name() + ')' - ) - ) - ); - - surfaceScalarField weights(tinterpScheme().weights(mag(psi))); - - const label nCmpts = pTraits::nComponents; - label localLoc1 = loc1; - label localLoc2 = loc2; - - // Get references to ldu fields of A matrix always as square - typename CoeffField::squareTypeField& blockDiag = - A.diag().asSquare(); - - // Insert boundary contributions directly into the matrix, this assumes - // that blockMatrix came from implicit grad or div operators! - psi.boundaryField().updateCoeffs(); - forAll(psi.boundaryField(), patchI) + // Insert coupling contributions into block matrix + forAll(bmesh, patchI) { - const fvPatchVectorField& pf = psi.boundaryField()[patchI]; - const fvPatch& patch = pf.patch(); - const vectorField& Sf = patch.Sf(); - const fvsPatchScalarField& pw = weights.boundaryField()[patchI]; - - vectorField internalCoeffs(pf.valueInternalCoeffs(pw)); - - if (patch.coupled()) + if (bmesh[patchI].coupled()) { typename CoeffField::squareTypeField& pcoupleUpper = A.coupleUpper()[patchI].asSquare(); typename CoeffField::squareTypeField& pcoupleLower = A.coupleLower()[patchI].asSquare(); - vectorField pcl = -pw*Sf; - vectorField pcu = pcl + Sf; + const typename CoeffField::linearTypeField& bmcu = + blockSystem.coupleUpper()[patchI].asLinear(); + const typename CoeffField::linearTypeField& bmcl = + blockSystem.coupleLower()[patchI].asLinear(); - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(pf, faceI) + for (label cmptI = 0; cmptI < nCmpts; cmptI++) { - label cellI = patch.faceCells()[faceI]; + forAll(bmcu, faceI) + { + pcoupleUpper[faceI](localLoc1, localLoc2) += + bmcu[faceI].component(cmptI); + pcoupleLower[faceI](localLoc1, localLoc2) += + bmcl[faceI].component(cmptI); + } - // Diag contribution - blockDiag[cellI](localLoc1, localLoc2) += cmptMultiply - ( - internalCoeffs[faceI], Sf[faceI] - ).component(cmptI); - - // Coupling contributions - pcoupleLower[faceI](localLoc1, localLoc2) -= - pcl[faceI].component(cmptI); - pcoupleUpper[faceI](localLoc1, localLoc2) -= - pcu[faceI].component(cmptI); + if (incFirst) + { + localLoc1++; + } + else + { + localLoc2++; + } } - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - - // Reset local locations for other patches - localLoc1 = loc1; - localLoc2 = loc2; - } - else - { - vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(pf, faceI) - { - label cellI = patch.faceCells()[faceI]; - - // Diag contribution - blockDiag[cellI](localLoc1, localLoc2) += cmptMultiply - ( - internalCoeffs[faceI], Sf[faceI] - ).component(cmptI); - - // Source contribution - b[cellI](localLoc1) -= cmptMultiply - ( - boundaryCoeffs[faceI], Sf[faceI] - ).component(cmptI); - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - // Reset local locations for other patches localLoc1 = loc1; localLoc2 = loc2; @@ -1220,15 +1052,14 @@ void insertBlockCoupling ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, + const BlockLduSystem& blockSystem, BlockLduMatrix& A, Field& b, const bool incFirst ) { - insertBlock(loc1, loc2, blockMatrix, A, incFirst); - insertBoundaryContributions(loc1, loc2, blockMatrix, psi, A, b, incFirst); + insertBlock(loc1, loc2, blockSystem, A, incFirst); + insertBoundaryContributions(loc1, loc2, blockSystem, A, b, incFirst); } @@ -1262,7 +1093,7 @@ void insertEquationCoupling } // Get references to fvScalarMatrix fields, updating boundary contributions - const scalarField& diag = matrix.D()(); + scalarField& diag = matrix.D(); scalarField& source = matrix.source(); matrix.addBoundarySource(source, false); diff --git a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H b/src/finiteVolume/blockMatrixTools/blockMatrixTools.H similarity index 81% rename from src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H rename to src/finiteVolume/blockMatrixTools/blockMatrixTools.H index 7000c0b4e..149262ebd 100644 --- a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H +++ b/src/finiteVolume/blockMatrixTools/blockMatrixTools.H @@ -25,7 +25,11 @@ InNamespace Foam:: Description - Block matrix insertion and retrieval tools + Block matrix insertion and retrieval tools. + Note: these will be obsolete since all of the functions are now member + functions in fvBlockMatrix class. All top level coupled solvers should + be rewritten using fvBlockMatrix - then we can get rid of this. + VV, 20/July/2014. SourceFiles blockMatrixTools.C @@ -36,7 +40,7 @@ SourceFiles #ifndef blockMatrixTools_H #define blockMatrixTools_H -#include "blockLduMatrices.H" +#include "BlockLduSystem.H" #include "fvMatrices.H" #include "surfaceFieldsFwd.H" @@ -218,60 +222,37 @@ namespace blockMatrixTools Field& b ); - // Insert block coupling - template + // Insert diagonal, lower and upper of blockMatrix into A + template void insertBlock ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, + const BlockLduSystem& blockMatrix, BlockLduMatrix& A, const bool incFirst ); - // Insert boundary contributions. These functions update the matrix coeffs - // according to boundary conditions. Two functions are needed for - // specialization if the field is scalar or vector. If the field is scalar, - // it is assumed that the grad(p) was used to obtain the matrix, and the - // boundary contributions to the source are vectors. If the field is - // vector, assumed is the div(U) and the source is scalar. - // Note parameter GeometricField, this is grad specialization. - template + // Insert source and coupling coeffs of block system into A and b + template void insertBoundaryContributions ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, + const BlockLduSystem& blockSystem, BlockLduMatrix& A, Field& b, const bool incFirst ); - // Note parameter GeometricField, this is div specialization. - template - void insertBoundaryContributions - ( - const label loc1, - const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, - BlockLduMatrix& A, - Field& b, - const bool incFirst - ); - - // Insert existing block matrix (obtained by implicit grad/div operator) - // into block system, updating the matrix coefficients from boundaries. - // To be used ONLY for pressure velocity coupling terms, other scenarios - // are not considered at the moment. + // Insert existing block system (obtained by implicit grad/div operator) + // into block system. template void insertBlockCoupling ( const label loc1, const label loc2, - const BlockLduMatrix& blockMatrix, - GeometricField& psi, + const BlockLduSystem& blockSystem, BlockLduMatrix& A, Field& b, const bool incFirst diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/calculatedFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/coupledFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/cyclicFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/emptyFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedGradientFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fixedValueFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/fvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/fvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/genericFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/genericFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/genericFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/genericFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/processorFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/transformFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/wedgeFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.C b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.C rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/fvPatchVectorNFields/zeroGradientFvPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/calculatedFvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/coupledFvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/emptyFvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/fvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/processorFvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.C rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields/wedgeFvsPatchVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFields.C b/src/finiteVolume/fields/surfaceFields/surfaceVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFields.C rename to src/finiteVolume/fields/surfaceFields/surfaceVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFields.H b/src/finiteVolume/fields/surfaceFields/surfaceVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFields.H rename to src/finiteVolume/fields/surfaceFields/surfaceVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFieldsFwd.H b/src/finiteVolume/fields/surfaceFields/surfaceVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/surfaceFields/surfaceVectorNFieldsFwd.H rename to src/finiteVolume/fields/surfaceFields/surfaceVectorNFieldsFwd.H diff --git a/src/VectorN/finiteVolume/fields/volFields/volVectorNFields.C b/src/finiteVolume/fields/volFields/volVectorNFields.C similarity index 100% rename from src/VectorN/finiteVolume/fields/volFields/volVectorNFields.C rename to src/finiteVolume/fields/volFields/volVectorNFields.C diff --git a/src/VectorN/finiteVolume/fields/volFields/volVectorNFields.H b/src/finiteVolume/fields/volFields/volVectorNFields.H similarity index 100% rename from src/VectorN/finiteVolume/fields/volFields/volVectorNFields.H rename to src/finiteVolume/fields/volFields/volVectorNFields.H diff --git a/src/VectorN/finiteVolume/fields/volFields/volVectorNFieldsFwd.H b/src/finiteVolume/fields/volFields/volVectorNFieldsFwd.H similarity index 100% rename from src/VectorN/finiteVolume/fields/volFields/volVectorNFieldsFwd.H rename to src/finiteVolume/fields/volFields/volVectorNFieldsFwd.H diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C index 714fe3fac..8e3ce3640 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C @@ -98,6 +98,35 @@ divScheme::~divScheme() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +tmp +< + BlockLduSystem::type> +> +divScheme::fvmDiv +( + const GeometricField& vf +) const +{ + FatalErrorIn + ( + "tmp divScheme::fvmDiv\n" + "(\n" + " 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; +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H index dcaf936b6..06705323d 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H @@ -41,7 +41,7 @@ SourceFiles #include "linear.H" #include "typeInfo.H" #include "runTimeSelectionTables.H" -#include "blockLduMatrices.H" +#include "BlockLduSystem.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -152,12 +152,15 @@ public: const GeometricField& ) = 0; - //- Return the blockLduMatrix corresponding to the implicit - // div discretization. For block coupled system. - virtual tmp fvmDiv + //- Return the BlockLduSystem corresponding to the implicit div + // discretization. For block coupled system. + virtual tmp + < + BlockLduSystem::type> + > fvmDiv ( const GeometricField& - ) const = 0; + ) const; }; diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divSchemes.C b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divSchemes.C index 89fb0b556..2b2421260 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divSchemes.C +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divSchemes.C @@ -61,13 +61,13 @@ defineTemplateRunTimeSelectionTable defineTemplateRunTimeSelectionTable ( - divScheme, + divScheme, Istream ); defineTemplateRunTimeSelectionTable ( - divScheme, + divScheme, Istream ); diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C index 1457139fb..c2888a554 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C @@ -69,38 +69,31 @@ gaussDivScheme::fvcDiv template -tmp gaussDivScheme::fvmDiv +tmp +< + BlockLduSystem::type> +> gaussDivScheme::fvmDiv ( const GeometricField& vf ) const { - tmp tweights = this->tinterpScheme_().weights(vf); - const scalarField& wIn = tweights().internalField(); + FatalErrorIn + ( + "tmp fvmDiv\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit div operator defined only for vector." + << abort(FatalError); - const fvMesh& mesh = vf.mesh(); + typedef typename innerProduct::type DivType; - tmp tbm - ( - new blockVectorMatrix - ( - mesh - ) - ); - blockVectorMatrix& bm = tbm(); + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); - // Grab ldu parts of block matrix as linear always - typename CoeffField::linearTypeField& u = bm.upper().asLinear(); - typename CoeffField::linearTypeField& l = bm.lower().asLinear(); - - const vectorField& SfIn = mesh.Sf().internalField(); - - l = -wIn*SfIn; - u = l + SfIn; - bm.negSumDiag(); - - // Interpolation schemes with corrections not accounted for - - return tbm; + return tbs; } diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H index 85de6494a..bba569e7c 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H @@ -97,9 +97,12 @@ public: const GeometricField& ); - //- Return the blockLduMatrix corresponding to the implicit - // div discretization. For block coupled system. - tmp fvmDiv + //- Return the BlockLduSystem corresponding to the implicit div + // discretization. For block coupled system. + tmp + < + BlockLduSystem::type> + > fvmDiv ( const GeometricField& ) const; diff --git a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C new file mode 100644 index 000000000..09312b4cf --- /dev/null +++ b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C @@ -0,0 +1,130 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Description + Specialisation of gaussDivScheme for vectors. Needed for implicit fvmDiv + operator for block coupled systems. + +\*---------------------------------------------------------------------------*/ + +#include "vectorGaussDivScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > gaussDivScheme::fvmDiv +( + 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 + typename CoeffField::linearTypeField& d = bs.diag().asLinear(); + typename CoeffField::linearTypeField& u = bs.upper().asLinear(); + typename CoeffField::linearTypeField& l = bs.lower().asLinear(); + + const vectorField& SfIn = mesh.Sf().internalField(); + + l = -wIn*SfIn; + u = l + 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 vectorField internalCoeffs(pf.valueInternalCoeffs(pw)); + + // Diag contribution + forAll(pf, faceI) + { + d[fc[faceI]] += cmptMultiply(internalCoeffs[faceI], Sf[faceI]); + } + + if (patch.coupled()) + { + typename CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const vectorField pcl = -pw*Sf; + const vectorField pcu = pcl + Sf; + + // Coupling contributions + pcoupleLower -= pcl; + pcoupleUpper -= pcu; + } + else + { + const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); + + // Boundary contribution + forAll(pf, faceI) + { + source[fc[faceI]] -= boundaryCoeffs[faceI] & Sf[faceI]; + } + } + } + + // Interpolation schemes with corrections not accounted for + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H new file mode 100644 index 000000000..a09f30589 --- /dev/null +++ b/src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H @@ -0,0 +1,72 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Typedef + Foam::fv::vectorGaussDivScheme + +Description + Specialisation of gaussDivScheme for vectors. Needed for implicit fvmDiv + operator for block coupled systems. + +SourceFiles + vectorGaussDivScheme.C + +\*---------------------------------------------------------------------------*/ + +#ifndef vectorGaussDivScheme_H +#define vectorGaussDivScheme_H + +#include "gaussDivScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > gaussDivScheme::fvmDiv +( + const GeometricField& +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C index bd6ea7772..25a9028d1 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C @@ -98,8 +98,12 @@ div return Div; } + template -tmp > div +tmp +< + BlockLduSystem::type> +> div ( GeometricField& vf, const word& name @@ -114,7 +118,10 @@ tmp > div template -tmp > div +tmp +< + BlockLduSystem::type> +> div ( GeometricField& vf ) diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H index a20bff42d..465bae1d1 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H @@ -39,7 +39,7 @@ SourceFiles #include "surfaceFieldsFwd.H" #include "surfaceInterpolationScheme.H" #include "fvMatrices.H" -#include "blockLduMatrices.H" +#include "BlockLduSystem.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -85,14 +85,20 @@ namespace fvm // Implicit div operator for block systems template - tmp div + tmp + < + BlockLduSystem::type> + > div ( GeometricField&, const word& ); template - tmp div + tmp + < + BlockLduSystem::type> + > div ( GeometricField& ); diff --git a/src/finiteVolume/finiteVolume/fvm/fvmGrad.C b/src/finiteVolume/finiteVolume/fvm/fvmGrad.C index 9bd1fa27d..e2dda7b74 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmGrad.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmGrad.C @@ -1,25 +1,26 @@ /*---------------------------------------------------------------------------*\ ========= | - \\ / F ield | foam-extend: Open Source CFD + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | For copyright notice see file Copyright + \\ / A nd | Copyright held by original author \\/ M anipulation | ------------------------------------------------------------------------------- License - This file is part of foam-extend. + This file is part of OpenFOAM. - foam-extend is free software: you can redistribute it and/or modify it + 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 + Free Software Foundation; either version 2 of the License, or (at your option) any later version. - foam-extend 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. + 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 foam-extend. If not, see . + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA \*---------------------------------------------------------------------------*/ @@ -39,7 +40,10 @@ namespace fvm // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template -tmp grad +tmp +< + BlockLduSystem::type> +> grad ( GeometricField& vf, const word& name @@ -54,7 +58,10 @@ tmp grad template -tmp grad +tmp +< + BlockLduSystem::type> +> grad ( GeometricField& vf ) diff --git a/src/finiteVolume/finiteVolume/fvm/fvmGrad.H b/src/finiteVolume/finiteVolume/fvm/fvmGrad.H index b398d1ec8..ff9bb6ef1 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmGrad.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmGrad.H @@ -1,33 +1,37 @@ /*---------------------------------------------------------------------------*\ ========= | - \\ / F ield | foam-extend: Open Source CFD + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | For copyright notice see file Copyright + \\ / A nd | Copyright held by original author \\/ M anipulation | ------------------------------------------------------------------------------- License - This file is part of foam-extend. + This file is part of OpenFOAM. - foam-extend is free software: you can redistribute it and/or modify it + 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 + Free Software Foundation; either version 2 of the License, or (at your option) any later version. - foam-extend 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. + 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 foam-extend. If not, see . + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA InNamespace Foam::fvm Description - Calculate the blockLduMatrx for the gradient of the field. - BlockLduMatrix will always be of type . Intended use: block - coupled solvers. i.e. implicit grad(p) in momentum equation. + Calculate the BlockLduSystem (matrix & source) for the gradient of the field + Intended use: block coupled solvers. i.e. implicit grad(p) in momentum + equation. + +Author + Vuko Vukcevic, FMENA Zagreb. SourceFiles fvmGrad.C @@ -38,7 +42,7 @@ SourceFiles #define fvmGrad_H #include "volFieldsFwd.H" -#include "blockLduMatrices.H" +#include "BlockLduSystem.H" #include "geometricOneField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -53,14 +57,20 @@ namespace Foam namespace fvm { template - tmp grad + tmp + < + BlockLduSystem::type> + > grad ( GeometricField&, const word& ); template - tmp grad + tmp + < + BlockLduSystem::type> + > grad ( GeometricField& ); diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C index c27e5fd0a..bca98cbaf 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C @@ -144,40 +144,32 @@ gaussGrad::grad return tgGrad; } - template -tmp gaussGrad::fvmGrad +tmp +< + BlockLduSystem::type> +> gaussGrad::fvmGrad ( const GeometricField& vf ) const { - tmp tweights = this->tinterpScheme_().weights(vf); - const scalarField& wIn = tweights().internalField(); - - const fvMesh& mesh = vf.mesh(); - - tmp tbm + FatalErrorIn ( - new blockVectorMatrix - ( - mesh - ) + "tmp fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operator defined only for scalar." + << abort(FatalError); + + typedef typename outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) ); - blockVectorMatrix& bm = tbm(); - // Grab ldu parts of block matrix as linear always - typename CoeffField::linearTypeField& u = bm.upper().asLinear(); - typename CoeffField::linearTypeField& l = bm.lower().asLinear(); - - const vectorField& SfIn = mesh.Sf().internalField(); - - l = -wIn*SfIn; - u = l + SfIn; - bm.negSumDiag(); - - // Interpolation schemes with corrections not accounted for - - return tbm; + return tbs; } diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H index e85e47aa8..c9bd2ef14 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H @@ -127,7 +127,6 @@ public: const GeometricField& ); - //- Return the gradient of the given field calculated // using Gauss' theorem on the interpolated field tmp @@ -139,14 +138,16 @@ public: const GeometricField& ) const; - //- Return the blockLduMatrix corresponding to the implicit - // grad discretization. For block coupled systems. - tmp fvmGrad + //- Return the BlockLduSystem corresponding to the implicit grad + // discretization. For block coupled systems. + tmp + < + BlockLduSystem::type> + > fvmGrad ( const GeometricField& ) const; - // Correct boundary conditions moved to base class // HJ, 14/Jun/2013 }; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C index 91bebf91b..9bb85151a 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C @@ -99,29 +99,33 @@ gradScheme::~gradScheme() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -tmp gradScheme::fvmGrad +tmp +< + BlockLduSystem::type> +> +gradScheme::fvmGrad ( const GeometricField& vf ) const { - FatalIOErrorIn + FatalErrorIn ( - "tmp fvmGrad\n", + "tmp gradScheme::fvmGrad\n" "(\n" - " GeometricField&" + " GeometricField&" ")\n" - ) << "Implicit gradient operator currently defined only for Gauss grad." - << abort(FatalIOError); + ) << "Implicit gradient operator currently defined only for Gauss linear " + << "and leastSquares (cell and face limiters are optional)." + << abort(FatalError); - tmp tbm + typedef typename outerProduct::type GradType; + + tmp > tbs ( - new blockVectorMatrix - ( - vf.mesh() - ) + new BlockLduSystem(vf.mesh()) ); - return tbm; + return tbs; } diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H index fcfab3a11..50b766b84 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H @@ -40,7 +40,7 @@ SourceFiles #include "surfaceFieldsFwd.H" #include "typeInfo.H" #include "runTimeSelectionTables.H" -#include "blockLduMatrices.H" +#include "BlockLduSystem.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -138,9 +138,12 @@ public: ) const = 0; - //- Return the blockLduMatrix corresponding to the implicit - // grad discretization. For block coupled systems. - virtual tmp fvmGrad + //- Return the BlockLduSystem corresponding to the implicit grad + // discretization. For block coupled systems. + virtual tmp + < + BlockLduSystem::type> + > fvmGrad ( const GeometricField& ) const; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C index 5a1157445..51ac2caff 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C @@ -121,8 +121,8 @@ leastSquaresGrad::grad forAll(neiVsf, patchFaceI) { lsGrad[faceCells[patchFaceI]] += - patchOwnLs[patchFaceI]* - (neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]); + patchOwnLs[patchFaceI] + *(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]); } } else @@ -132,8 +132,8 @@ leastSquaresGrad::grad forAll(patchVsf, patchFaceI) { lsGrad[faceCells[patchFaceI]] += - patchOwnLs[patchFaceI]* - (patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]); + patchOwnLs[patchFaceI] + *(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]); } } } @@ -145,6 +145,34 @@ leastSquaresGrad::grad return tlsGrad; } +template +tmp +< + BlockLduSystem::type> +> leastSquaresGrad::fvmGrad +( + const GeometricField& vf +) const +{ + FatalErrorIn + ( + "tmp fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operator defined only for scalar." + << abort(FatalError); + + typedef typename outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vf.mesh()) + ); + + return tbs; +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H index 337db5cd8..4d95e1d58 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H @@ -96,6 +96,16 @@ public: ( const GeometricField& ) const; + + //- Return the BlockLduSystem corresponding to the implicit least + // quares grad discretization. For block coupled systems. + tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& + ) const; }; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C index ecbf38f43..e2a6a76fa 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C @@ -101,19 +101,15 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const // Set local references to mesh data const unallocLabelList& owner = mesh().owner(); const unallocLabelList& neighbour = mesh().neighbour(); - const volVectorField& C = mesh().C(); - const surfaceScalarField& w = mesh().weights(); -// const surfaceScalarField& magSf = mesh().magSf(); - // Set up temporary storage for the dd tensor (before inversion) symmTensorField dd(mesh().nCells(), symmTensor::zero); - forAll(owner, facei) + forAll(owner, faceI) { - label own = owner[facei]; - label nei = neighbour[facei]; + label own = owner[faceI]; + label nei = neighbour[faceI]; vector d = C[nei] - C[own]; symmTensor wdd = (1.0/magSqr(d))*sqr(d); @@ -122,113 +118,88 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const dd[nei] += wdd; } - - forAll(lsP.boundaryField(), patchi) + forAll(lsP.boundaryField(), patchI) { - const fvsPatchScalarField& pw = w.boundaryField()[patchi]; - // Note: least squares in 1.4.1 and other OpenCFD versions - // are wrong because of incorrect weighting. HJ, 23/Oct/2008 -// const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi]; - - const fvPatch& p = pw.patch(); + const fvPatch& p = mesh().boundary()[patchI]; const unallocLabelList& faceCells = p.patch().faceCells(); - // Build the d-vectors - - // Original version: closest distance to boundary - vectorField pd = - mesh().Sf().boundaryField()[patchi] - /( - mesh().magSf().boundaryField()[patchi] - *mesh().deltaCoeffs().boundaryField()[patchi] - ); - - if (!mesh().orthogonal()) - { - pd -= mesh().correctionVectors().boundaryField()[patchi] - /mesh().deltaCoeffs().boundaryField()[patchi]; - } - // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 - // Experimental: review fixed gradient condition. HJ, 30/Sep/2010 -// vectorField pd = p.delta(); + const vectorField pd = p.delta(); - if (p.coupled()) + forAll(pd, patchFaceI) { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; - - dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d); - } - } - else - { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; - - dd[faceCells[patchFacei]] += (1.0/magSqr(d))*sqr(d); - } + const vector& d = pd[patchFaceI]; + dd[faceCells[patchFaceI]] += (1.0/magSqr(d))*sqr(d); } } - - // Invert the dd tensor -// symmTensorField invDd = inv(dd); - // Fix: householder inverse. HJ, 3/Nov/2009 - symmTensorField invDd = hinv(dd); - + // For easy access of neighbour coupled patch field needed for + // lsN vectors on implicitly coupled boundaries. VV, 18/June/2014 + volSymmTensorField volInvDd + ( + IOobject + ( + "volInvDd", + mesh().pointsInstance(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh(), + dimensionedSymmTensor("zero", dimless, symmTensor::zero), + "zeroGradient" + ); + symmTensorField& invDd = volInvDd.internalField(); +// invDd = inv(dd); + invDd = hinv(dd); // Revisit all faces and calculate the lsP and lsN vectors - forAll(owner, facei) + forAll(owner, faceI) { - label own = owner[facei]; - label nei = neighbour[facei]; + label own = owner[faceI]; + label nei = neighbour[faceI]; vector d = C[nei] - C[own]; scalar magSfByMagSqrd = 1.0/magSqr(d); - lsP[facei] = magSfByMagSqrd*(invDd[own] & d); - lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d); + lsP[faceI] = magSfByMagSqrd*(invDd[own] & d); + lsN[faceI] = -magSfByMagSqrd*(invDd[nei] & d); } - forAll(lsP.boundaryField(), patchi) + forAll(lsP.boundaryField(), patchI) { - fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi]; - - const fvsPatchScalarField& pw = w.boundaryField()[patchi]; - // Note: least squares in 1.4.1 and other OpenCFD versions - // are wrong because of incorrect weighting. HJ, 23/Oct/2008 -// const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi]; - - const fvPatch& p = pw.patch(); + fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI]; + fvsPatchVectorField& patchLsN = lsN.boundaryField()[patchI]; + const fvPatch& p = mesh().boundary()[patchI]; const unallocLabelList& faceCells = p.faceCells(); - // Build the d-vectors // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 - vectorField pd = p.delta(); + const vectorField pd = p.delta(); if (p.coupled()) { - forAll(pd, patchFacei) - { - const vector& d = pd[patchFacei]; + const symmTensorField invDdNei = + volInvDd.boundaryField()[patchI].patchNeighbourField(); - patchLsP[patchFacei] = - (1.0/magSqr(d)) - *(invDd[faceCells[patchFacei]] & d); + forAll(pd, patchFaceI) + { + const vector& d = pd[patchFaceI]; + + patchLsP[patchFaceI] = (1.0/magSqr(d)) + *(invDd[faceCells[patchFaceI]] & d); + + patchLsN[patchFaceI] = - (1.0/magSqr(d)) + *(invDdNei[faceCells[patchFaceI]] & d); } } else { - forAll(pd, patchFacei) + forAll(pd, patchFaceI) { - const vector& d = pd[patchFacei]; + const vector& d = pd[patchFaceI]; - patchLsP[patchFacei] = - (1.0/magSqr(d)) - *(invDd[faceCells[patchFacei]] & d); + patchLsP[patchFaceI] = (1.0/magSqr(d)) + *(invDd[faceCells[patchFaceI]] & d); } } } diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H index dc573b7ae..c5b5a2d94 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H @@ -125,6 +125,17 @@ public: ( const GeometricField& ) const; + + + //- Return the BlockLduSystem corresponding to the implicit cell + // limted grad discretization. For block coupled systems. + tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& + ) const; }; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C index 77d0a3b24..e6392aeb6 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C @@ -404,6 +404,257 @@ tmp cellLimitedGrad::grad } +template<> +tmp > cellLimitedGrad::fvmGrad +( + const volScalarField& vsf +) const +{ + // Consider doing a calculateLimiter member function since both fvmGrad and + // grad use almost the same procedure to calculate limiter. VV, 9/June/2014. + const fvMesh& mesh = vsf.mesh(); + + tmp > tbs = basicGradScheme_().fvmGrad(vsf); + + if (k_ < SMALL) + { + return tbs; + } + + BlockLduSystem& bs = tbs(); + + // Calculate current gradient for explicit limiting + tmp tGrad = basicGradScheme_().grad(vsf); + const volVectorField& g = tGrad(); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + const volVectorField& C = mesh.C(); + const surfaceVectorField& Cf = mesh.Cf(); + + scalarField maxVsf(vsf.internalField()); + scalarField minVsf(vsf.internalField()); + + forAll(owner, facei) + { + label own = owner[facei]; + label nei = neighbour[facei]; + + scalar vsfOwn = vsf[own]; + scalar vsfNei = vsf[nei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + + maxVsf[nei] = max(maxVsf[nei], vsfOwn); + minVsf[nei] = min(minVsf[nei], vsfOwn); + } + + + const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField(); + + forAll(bsf, patchi) + { + const fvPatchScalarField& psf = bsf[patchi]; + + const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells(); + + if (psf.coupled()) + { + scalarField psfNei = psf.patchNeighbourField(); + + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + scalar vsfNei = psfNei[pFacei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + } + } + else + { + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + scalar vsfNei = psf[pFacei]; + + maxVsf[own] = max(maxVsf[own], vsfNei); + minVsf[own] = min(minVsf[own], vsfNei); + } + } + } + + maxVsf -= vsf; + minVsf -= vsf; + + if (k_ < 1.0) + { + scalarField maxMinVsf = (1.0/k_ - 1.0)*(maxVsf - minVsf); + maxVsf += maxMinVsf; + minVsf -= maxMinVsf; + + //maxVsf *= 1.0/k_; + //minVsf *= 1.0/k_; + } + + + // Create limiter as volScalarField for proper update of coupling coeffs + volScalarField limitField + ( + IOobject + ( + "limitField", + vsf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1), + "zeroGradient" + ); + scalarField& lfIn = limitField.internalField(); + + forAll(owner, facei) + { + label own = owner[facei]; + label nei = neighbour[facei]; + + // owner side + limitFace + ( + lfIn[own], + maxVsf[own], + minVsf[own], + (Cf[facei] - C[own]) & g[own] + ); + + // neighbour side + limitFace + ( + lfIn[nei], + maxVsf[nei], + minVsf[nei], + (Cf[facei] - C[nei]) & g[nei] + ); + } + + forAll(bsf, patchi) + { + const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells(); + const vectorField& pCf = Cf.boundaryField()[patchi]; + + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + + limitFace + ( + lfIn[own], + maxVsf[own], + minVsf[own], + (pCf[pFacei] - C[own]) & g[own] + ); + } + } + + limitField.correctBoundaryConditions(); + + if (fv::debug) + { + Info<< "gradient limiter for: " << vsf.name() + << " max = " << gMax(lfIn) + << " min = " << gMin(lfIn) + << " average: " << gAverage(lfIn) << endl; + } + + 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(); + + // Limit upper and lower coeffs + forAll(u, faceI) + { + label own = owner[faceI]; + label nei = neighbour[faceI]; + + u[faceI] *= lfIn[own]; + l[faceI] *= lfIn[nei]; + } + + // Limit diag and source coeffs + forAll(d, cellI) + { + d[cellI] *= lfIn[cellI]; + source[cellI] *= lfIn[cellI]; + } + + // Limit coupling coeffs + forAll(vsf.boundaryField(), patchI) + { + const fvPatchScalarField& pf = vsf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + + if (patch.coupled()) + { + typename CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const scalarField lfNei = + limitField.boundaryField()[patchI].patchNeighbourField(); + + forAll(pf, faceI) + { + label cellI = patch.faceCells()[faceI]; + + pcoupleUpper[faceI] *= lfIn[cellI]; + pcoupleLower[faceI] *= lfNei[faceI]; + } + } + } + + return tbs; +} + + +template<> +tmp +< + BlockLduSystem::type> +> +cellLimitedGrad::fvmGrad +( + const volVectorField& vsf +) const +{ + FatalErrorIn + ( + "tmp cellLimitedGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with cell limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef typename outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vsf.mesh()) + ); + + return tbs; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H index 6718583b2..ca4eee14b 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H @@ -125,6 +125,17 @@ public: ( const GeometricField& ) const; + + //- Return the BlockLduSystem corresponding to the implicit face + // limted grad discretization. For block coupled systems. + tmp + < + BlockLduSystem::type> + > fvmGrad + ( + const GeometricField& + ) const; + }; diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C index 4a5599f21..b6a2a97fc 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C @@ -362,6 +362,240 @@ tmp faceLimitedGrad::grad } +template<> +tmp > faceLimitedGrad::fvmGrad +( + const volScalarField& vsf +) const +{ + // Consider doing a calculateLimiter member function since both fvmGrad and + // grad use almost the same procedure to calculate limiter. VV, 9/June/2014. + const fvMesh& mesh = vsf.mesh(); + + tmp > tbs = basicGradScheme_().fvmGrad(vsf); + + if (k_ < SMALL) + { + return tbs; + } + + BlockLduSystem& bs = tbs(); + + // Calculate current gradient for explicit limiting + tmp tGrad = basicGradScheme_().grad(vsf); + const volVectorField& g = tGrad(); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + const volVectorField& C = mesh.C(); + const surfaceVectorField& Cf = mesh.Cf(); + + // Create limiter as volScalarField for proper update of coupling coeffs + volScalarField limitField + ( + IOobject + ( + "limitField", + vsf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1), + "zeroGradient" + ); + scalarField& lfIn = limitField.internalField(); + + scalar rk = (1.0/k_ - 1.0); + + forAll(owner, facei) + { + label own = owner[facei]; + label nei = neighbour[facei]; + + scalar vsfOwn = vsf[own]; + scalar vsfNei = vsf[nei]; + + scalar maxFace = max(vsfOwn, vsfNei); + scalar minFace = min(vsfOwn, vsfNei); + scalar maxMinFace = rk*(maxFace - minFace); + maxFace += maxMinFace; + minFace -= maxMinFace; + + // owner side + limitFace + ( + lfIn[own], + maxFace - vsfOwn, minFace - vsfOwn, + (Cf[facei] - C[own]) & g[own] + ); + + // neighbour side + limitFace + ( + lfIn[nei], + maxFace - vsfNei, minFace - vsfNei, + (Cf[facei] - C[nei]) & g[nei] + ); + } + + const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField(); + + forAll(bsf, patchi) + { + const fvPatchScalarField& psf = bsf[patchi]; + + const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells(); + const vectorField& pCf = Cf.boundaryField()[patchi]; + + if (psf.coupled()) + { + scalarField psfNei = psf.patchNeighbourField(); + + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + + scalar vsfOwn = vsf[own]; + scalar vsfNei = psfNei[pFacei]; + + scalar maxFace = max(vsfOwn, vsfNei); + scalar minFace = min(vsfOwn, vsfNei); + scalar maxMinFace = rk*(maxFace - minFace); + maxFace += maxMinFace; + minFace -= maxMinFace; + + limitFace + ( + lfIn[own], + maxFace - vsfOwn, minFace - vsfOwn, + (pCf[pFacei] - C[own]) & g[own] + ); + } + } + else if (psf.fixesValue()) + { + forAll(pOwner, pFacei) + { + label own = pOwner[pFacei]; + + scalar vsfOwn = vsf[own]; + scalar vsfNei = psf[pFacei]; + + scalar maxFace = max(vsfOwn, vsfNei); + scalar minFace = min(vsfOwn, vsfNei); + scalar maxMinFace = rk*(maxFace - minFace); + maxFace += maxMinFace; + minFace -= maxMinFace; + + limitFace + ( + lfIn[own], + maxFace - vsfOwn, minFace - vsfOwn, + (pCf[pFacei] - C[own]) & g[own] + ); + } + } + } + + if (fv::debug) + { + Info<< "gradient limiter for: " << vsf.name() + << " max = " << gMax(lfIn) + << " min = " << gMin(lfIn) + << " average: " << gAverage(lfIn) << endl; + } + + 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(); + + // Limit upper and lower coeffs + forAll(u, faceI) + { + label own = owner[faceI]; + label nei = neighbour[faceI]; + + u[faceI] *= lfIn[own]; + l[faceI] *= lfIn[nei]; + } + + // Limit diag and source coeffs + forAll(d, cellI) + { + d[cellI] *= lfIn[cellI]; + source[cellI] *= lfIn[cellI]; + } + + // Limit coupling coeffs + forAll(vsf.boundaryField(), patchI) + { + const fvPatchScalarField& pf = vsf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + + if (patch.coupled()) + { + typename CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const scalarField lfNei = + limitField.boundaryField()[patchI].patchNeighbourField(); + + forAll(pf, faceI) + { + label cellI = patch.faceCells()[faceI]; + + pcoupleUpper[faceI] *= lfIn[cellI]; + pcoupleLower[faceI] *= lfNei[faceI]; + } + } + } + + limitField.write(); + + return tbs; +} + + +template<> +tmp +< + BlockLduSystem::type> +> +faceLimitedGrad::fvmGrad +( + const volVectorField& vsf +) const +{ + FatalErrorIn + ( + "tmp faceLimitedGrad::fvmGrad\n" + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operators with face limiters defined only for " + << "scalar." + << abort(FatalError); + + typedef typename outerProduct::type GradType; + + tmp > tbs + ( + new BlockLduSystem(vsf.mesh()) + ); + + return tbs; +} + + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/.scalarGaussGrad.C.swp b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/.scalarGaussGrad.C.swp new file mode 100644 index 000000000..a62419570 Binary files /dev/null and b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/.scalarGaussGrad.C.swp differ diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C new file mode 100644 index 000000000..801185958 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C @@ -0,0 +1,130 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Description + Specialisation of gaussGrad for scalars. Needed for implicit fvmGrad + operator for block coupled systems. + +\*---------------------------------------------------------------------------*/ + +#include "scalarGaussGrad.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > gaussGrad::fvmGrad +( + 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(); + 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(); + + const vectorField& SfIn = mesh.Sf().internalField(); + + l = -wIn*SfIn; + u = l + SfIn; + bs.negSumDiag(); + + // Boundary contributions + forAll(vf.boundaryField(), patchI) + { + const fvPatchScalarField& pf = vf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + const vectorField& Sf = patch.Sf(); + const fvsPatchScalarField& pw = tweights().boundaryField()[patchI]; + const labelList& fc = patch.faceCells(); + + const scalarField internalCoeffs(pf.valueInternalCoeffs(pw)); + + // Diag contribution + forAll(pf, faceI) + { + d[fc[faceI]] += internalCoeffs[faceI]*Sf[faceI]; + } + + if (patch.coupled()) + { + typename CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + const vectorField pcl = -pw*Sf; + const vectorField pcu = pcl + Sf; + + // Coupling contributions + pcoupleLower -= pcl; + pcoupleUpper -= pcu; + } + else + { + const scalarField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); + + // Boundary contribution + forAll(pf, faceI) + { + source[fc[faceI]] -= boundaryCoeffs[faceI]*Sf[faceI]; + } + } + } + + // Interpolation schemes with corrections not accounted for + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.H new file mode 100644 index 000000000..5b137883a --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.H @@ -0,0 +1,72 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Typedef + Foam::fv::scalarGaussGrad + +Description + Specialisation of gaussGrad for scalars. Needed for implicit fvmGrad + operator for block coupled systems. + +SourceFiles + scalarGaussGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarGaussGrad_H +#define scalarGaussGrad_H + +#include "gaussGrad.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > gaussGrad::fvmGrad +( + const GeometricField& +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C new file mode 100644 index 000000000..5f74e6240 --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C @@ -0,0 +1,171 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Description + Specialisation of leastSquaresGrad for scalars. Needed for implicit fvmGrad + operator for block coupled systems. + +\*---------------------------------------------------------------------------*/ + +#include "scalarLeastSquaresGrad.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > leastSquaresGrad::fvmGrad +( + const GeometricField& vf +) const +{ + const fvMesh& mesh = vf.mesh(); + + const unallocLabelList& own = mesh.owner(); + const unallocLabelList& nei = mesh.neighbour(); + + volScalarField cellV + ( + IOobject + ( + "cellV", + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimVolume, 0), + "zeroGradient" + ); + cellV.internalField() = mesh.V(); + cellV.correctBoundaryConditions(); + const scalarField& cellVIn = cellV.internalField(); + + const surfaceScalarField& w = mesh.weights(); + + tmp > tbs + ( + new BlockLduSystem(mesh) + ); + BlockLduSystem& bs = tbs(); + 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(); + + // Get references to least square vectors + const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh); + + const surfaceVectorField& ownLs = lsv.pVectors(); + const surfaceVectorField& neiLs = lsv.nVectors(); + + forAll(nei, faceI) + { + register label owner = own[faceI]; + register label neighbour = nei[faceI]; + + u[faceI] = cellVIn[owner]*ownLs[faceI]; + l[faceI] = cellVIn[neighbour]*neiLs[faceI]; + + // Caution - this is NOT negSumDiag(). VV, 17/July/2014. + d[owner] -= u[faceI]; + d[neighbour] -= l[faceI]; + } + + // Boundary contributions + forAll(vf.boundaryField(), patchI) + { + const fvPatchScalarField& pf = vf.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + const vectorField& pownLs = ownLs.boundaryField()[patchI]; + const fvsPatchScalarField& pw = w.boundaryField()[patchI]; + const labelList& fc = patch.faceCells(); + + // Part of diagonal contribution irrespective of the patch type + forAll(pf, faceI) + { + const label cellI = fc[faceI]; + d[cellI] -= cellVIn[cellI]*pownLs[faceI]; + } + + if (patch.coupled()) + { + const vectorField& pneiLs = neiLs.boundaryField()[patchI]; + const scalarField cellVInNei = + cellV.boundaryField()[patchI].patchNeighbourField(); + + typename CoeffField::linearTypeField& pcoupleUpper = + bs.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + bs.coupleLower()[patchI].asLinear(); + + // Coupling and diagonal contributions + forAll(pf, faceI) + { + const vector upper = cellVIn[fc[faceI]]*pownLs[faceI]; + const vector lower = cellVInNei[faceI]*pneiLs[faceI]; + + pcoupleUpper[faceI] -= upper; + pcoupleLower[faceI] -= lower; + } + } + else + { + const scalarField internalCoeffs(pf.valueInternalCoeffs(pw)); + const scalarField boundaryCoeffs(pf.valueBoundaryCoeffs(pw)); + + // Diagonal and source contributions depending on the patch type + forAll(pf, faceI) + { + const label cellI = fc[faceI]; + d[cellI] += cellVIn[cellI]*pownLs[faceI]*internalCoeffs[faceI]; + source[cellI] -= cellVIn[cellI]*pownLs[faceI]* + boundaryCoeffs[faceI]; + } + } + } + + return tbs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.H new file mode 100644 index 000000000..ff111d40b --- /dev/null +++ b/src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.H @@ -0,0 +1,72 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Typedef + Foam::fv::scalarLeastSquaresGrad + +Description + Specialisation of leastSquaresGrad for scalars. Needed for implicit fvmGrad + operator for block coupled systems. + +SourceFiles + scalarLeastSquaresGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarLeastSquaresGrad_H +#define scalarLeastSquaresGrad_H + +#include "leastSquaresGrad.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > leastSquaresGrad::fvmGrad +( + const GeometricField& +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C new file mode 100644 index 000000000..f2697f212 --- /dev/null +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C @@ -0,0 +1,1082 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "fvBlockMatrix.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +template +void fvBlockMatrix::insertSolutionVector +( + const direction dir, + const Field& xSingle +) +{ + // Get number of field components and local copy of direction, for + // consistency with member functions where directions need to be reset. + const direction nCmpts = pTraits::nComponents; + direction localDir = dir; + + Field& psiIn = psi_.internalField(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField xSingleCurr(xSingle.component(cmptI)); + + forAll (xSingleCurr, cellI) + { + psiIn[cellI](localDir) = xSingleCurr[cellI]; + } + + localDir++; + } +} + + +template +template +void fvBlockMatrix::insertDiagSource +( + const direction dir, + fvMatrix& matrix +) +{ + matrix.completeAssembly(); + + // Save a copy for different components + scalarField& diag = matrix.diag(); + scalarField saveDiag(diag); + + // Add source boundary contribution + Field& source = matrix.source(); + matrix.addBoundarySource(source, false); + + const direction nCmpts = pTraits::nComponents; + direction localDir = dir; + + // Get reference to this source field of block system + Field& b = this->source(); + + if + ( + // This is needed if the matrixType is , then you need to grab + // coeffs as linear. Consider doing a matrixType check also. + // VV, 17/March/2014 + this->diag().activeType() != blockCoeffBase::SQUARE + ) + { + typename CoeffField::linearTypeField& blockDiag = + this->diag().asLinear(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + matrix.addBoundaryDiag(diag, cmptI); + scalarField sourceCmpt(source.component(cmptI)); + +// FieldField bouCoeffsCmpt +// ( +// matrix.boundaryCoeffs().component(cmptI) +// ); + + // Possible problem for coupled non-aligned boundaries. + // VV, 14/May/2014. +// matrix.correctImplicitBoundarySource +// ( +// bouCoeffsCmpt, +// sourceCmpt, +// cmptI +// ); + + forAll (diag, cellI) + { + blockDiag[cellI](localDir) = diag[cellI]; + b[cellI](localDir) += sourceCmpt[cellI]; + } + + localDir++; + + // Reset diagonal + diag = saveDiag; + } + } + else if (this->diag().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + matrix.addBoundaryDiag(diag, cmptI); + scalarField sourceCmpt(source.component(cmptI)); + +// FieldField bouCoeffsCmpt +// ( +// matrix.boundaryCoeffs().component(cmptI) +// ); + +// matrix.correctImplicitBoundarySource +// ( +// bouCoeffsCmpt, +// sourceCmpt, +// cmptI +// ); + + forAll (diag, cellI) + { + blockDiag[cellI](localDir, localDir) = diag[cellI]; + b[cellI](localDir) += sourceCmpt[cellI]; + } + + localDir++; + + // Reset diagonal + diag = saveDiag; + } + } +} + + +template +template +void fvBlockMatrix::insertUpperLower +( + const direction dir, + const fvMatrix& matrix +) +{ + if (matrix.diagonal()) + { + // Matrix for insertion is diagonal-only: nothing to do + return; + } + + const direction nCmpts = pTraits::nComponents; + direction localDir = dir; + + if (matrix.hasUpper()) + { + const scalarField& upper = matrix.upper(); + + if (this->upper().activeType() == blockCoeffBase::UNALLOCATED) + { + this->upper().asScalar() = upper; + } + else if + ( + this->upper().activeType() == blockCoeffBase::SCALAR + || this->upper().activeType() == blockCoeffBase::LINEAR + ) + { + typename CoeffField::linearTypeField& blockUpper = + this->upper().asLinear(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (upper, faceI) + { + blockUpper[faceI](localDir) = upper[faceI]; + } + + localDir++; + } + } + else if (this->upper().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockUpper = + this->upper().asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (upper, faceI) + { + blockUpper[faceI](localDir, localDir) = upper[faceI]; + } + + localDir++; + } + } + } + else + { + FatalErrorIn + ( + "void fvBlockMatrix::insertUpperLower\n" + "(\n" + " const direction dir,\n" + " const fvMatrix& matrix\n" + ")" + ) << "Error in matrix insertion: problem with block structure." + << abort(FatalError); + } + + if (matrix.symmetric() && this->symmetric()) + { + Info<< "Both matrices are symmetric: inserting only upper triangle" + << endl; + } + else + { + // Reset localDir + localDir = dir; + + // Either scalar or block matrix is asymmetric: insert lower triangle + const scalarField& lower = matrix.lower(); + + if (this->lower().activeType() == blockCoeffBase::UNALLOCATED) + { + this->lower().asScalar() = lower; + } + else if + ( + this->lower().activeType() == blockCoeffBase::SCALAR + || this->lower().activeType() == blockCoeffBase::LINEAR + ) + { + typename CoeffField::linearTypeField& blockLower = + this->lower().asLinear(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (lower, faceI) + { + blockLower[faceI](localDir) = lower[faceI]; + } + + localDir++; + } + } + else if (this->lower().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockLower = + this->lower().asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (lower, faceI) + { + blockLower[faceI](localDir, localDir) = lower[faceI]; + } + + localDir++; + } + } + } +} + + +template +template +void fvBlockMatrix::updateCouplingCoeffs +( + const direction dir, + const fvMatrix& matrix +) +{ + const direction nCmpts = pTraits::nComponents; + direction localDir = dir; + + const GeometricField& psi = matrix.psi(); + forAll(psi.boundaryField(), patchI) + { + const fvPatchField& pf = psi.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + + if (patch.coupled()) + { + const Field& icp = matrix.internalCoeffs()[patchI]; + const Field& bcp = matrix.boundaryCoeffs()[patchI]; + + if + ( + this->coupleUpper()[patchI].activeType() + != blockCoeffBase::SQUARE + ) + { + typename CoeffField::linearTypeField& pcoupleUpper = + this->coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + this->coupleLower()[patchI].asLinear(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField icpCmpt = icp.component(cmptI); + scalarField bcpCmpt = bcp.component(cmptI); + + forAll(pf, faceI) + { + pcoupleUpper[faceI](localDir) = bcpCmpt[faceI]; + pcoupleLower[faceI](localDir) = icpCmpt[faceI]; + } + + localDir++; + } + + localDir = dir; + } + else if + ( + this->coupleUpper()[patchI].activeType() + == blockCoeffBase::SQUARE + ) + { + typename CoeffField::squareTypeField& pcoupleUpper = + this->coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + this->coupleLower()[patchI].asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField icpCmpt = icp.component(cmptI); + scalarField bcpCmpt = bcp.component(cmptI); + + forAll(pf, faceI) + { + pcoupleUpper[faceI](localDir, localDir) = + bcpCmpt[faceI]; + pcoupleLower[faceI](localDir, localDir) = + icpCmpt[faceI]; + } + + localDir++; + } + + localDir = dir; + } + } + } +} + + +template +template +void fvBlockMatrix::insertBlock +( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir +) +{ + // Sanity checks + { + const direction blockMatrixSize = pTraits::nComponents; + const direction blockMatrixASize = pTraits::nComponents; + + if (blockMatrixSize > blockMatrixASize) + { + FatalErrorIn + ( + "void fvBlockMatrix::insertBlock\n" + "(\n" + " const direction dirI,\n" + " const direction dirJ,\n" + " BlockLduSystem& blockSystem,\n" + " const bool incrementColumnDir\n" + ")" + ) << "Trying to insert a block matrix from BlockLduSystem into " + << "smaller one from fvBlockMatrix." + << abort(FatalError); + } + + if (dirI == dirJ) + { + FatalErrorIn + ( + "void fvBlockMatrix::insertBlock\n" + "(\n" + " const direction dirI,\n" + " const direction dirJ,\n" + " BlockLduSystem& blockSystem,\n" + " const bool incrementColumnDir\n" + ")" + ) << "Trying to insert coupling in the position where equation " + << "should be, since dirI = dirJ. Try using insertEquation " + << "member function." + << abort(FatalError); + } + } + + const direction nCmpts = pTraits::nComponents; + direction localDirI = dirI; + direction localDirJ = dirJ; + + // Get references to ldu fields of blockMatrix always as linear + const typename CoeffField::linearTypeField& bmd = + blockSystem.diag().asLinear(); + const typename CoeffField::linearTypeField& bmu = + blockSystem.upper().asLinear(); + const typename CoeffField::linearTypeField& bml = + blockSystem.lower().asLinear(); + + // Get references to ldu fields of this block matrix always as square + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + typename CoeffField::squareTypeField& blockUpper = + this->upper().asSquare(); + typename CoeffField::squareTypeField& blockLower = + this->lower().asSquare(); + + // Insert blockMatrix that represents coupling into larger system matrix + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll(bmd, cellI) + { + blockDiag[cellI](localDirI, localDirJ) += + bmd[cellI].component(cmptI); + } + + forAll(bmu, faceI) + { + blockUpper[faceI](localDirI, localDirJ) += + bmu[faceI].component(cmptI); + blockLower[faceI](localDirI, localDirJ) += + bml[faceI].component(cmptI); + } + + if (incrementColumnDir) + { + localDirI++; + } + else + { + localDirJ++; + } + } +} + + +template +template +void fvBlockMatrix::insertBoundaryContributions +( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir +) +{ + // Need to get reference to fvMesh instead of lduMesh + const fvBoundaryMesh& bmesh = + refCast(this->mesh()).boundary(); + + const direction nCmpts = pTraits::nComponents; + const direction nSrcCmpts = pTraits::nComponents; + direction localDirI = dirI; + direction localDirJ = dirJ; + + const Field& source = blockSystem.source(); + + // Get reference to this source field of block system + Field& b = this->source(); + + // Insert source from block system to this system's rhs + for (direction cmptI = 0; cmptI < nSrcCmpts; cmptI++) + { + scalarField sourceCmpt(source.component(cmptI)); + + forAll(b, cellI) + { + b[cellI](localDirI) += sourceCmpt[cellI]; + } + + if (incrementColumnDir) + { + localDirI++; + } + else + { + localDirJ++; + } + } + + // Reset local directions for coupling contributions + localDirI = dirI; + localDirJ = dirJ; + + // Insert coupling contributions into block matrix + forAll(bmesh, patchI) + { + if (bmesh[patchI].coupled()) + { + typename CoeffField::squareTypeField& pcoupleUpper = + this->coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + this->coupleLower()[patchI].asSquare(); + + const typename CoeffField::linearTypeField& bmcu = + blockSystem.coupleUpper()[patchI].asLinear(); + const typename CoeffField::linearTypeField& bmcl = + blockSystem.coupleLower()[patchI].asLinear(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll(bmcu, faceI) + { + pcoupleUpper[faceI](localDirI, localDirJ) += + bmcu[faceI].component(cmptI); + pcoupleLower[faceI](localDirI, localDirJ) += + bmcl[faceI].component(cmptI); + } + + if (incrementColumnDir) + { + localDirI++; + } + else + { + localDirJ++; + } + } + + // Reset local directions for other patches + localDirI = dirI; + localDirJ = dirJ; + } + } +} + + +template +void fvBlockMatrix::insertCouplingDiagSource +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix +) +{ + // Prepare the diagonal and source + scalarField diag = matrix.diag(); + scalarField source = matrix.source(); + + // Add boundary source contribution + matrix.addBoundaryDiag(diag, 0); + matrix.addBoundarySource(source, false); + + // Get reference to block diagonal of the block system + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + + // Get reference to this source field of the block system + Field& b = this->source(); + + // Set off-diagonal coefficient + forAll(diag, cellI) + { + blockDiag[cellI](dirI, dirJ) += diag[cellI]; + b[cellI](dirI) += source[cellI]; + } +} + + +template +void fvBlockMatrix::insertCouplingUpperLower +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix +) +{ + if (matrix.diagonal()) + { + // Matrix for insertion is diagonal-only: nothing to do + return; + } + + if (matrix.symmetric() && this->symmetric()) + { + Info<< "Both fvScalarMatrix and block matrix are symmetric: " << nl + << "inserting only upper triangle" + << endl; + } + else + { + // Either scalar or block matrix is asymmetric: insert lower triangle + const scalarField& lower = matrix.lower(); + + typename CoeffField::squareTypeField& blockLower = + this->lower().asSquare(); + + forAll (lower, cellI) + { + blockLower[cellI](dirI, dirJ) = lower[cellI]; + } + } + + if (matrix.hasUpper()) + { + const scalarField& upper = matrix.upper(); + + typename CoeffField::squareTypeField& blockUpper = + this->upper().asSquare(); + + forAll(upper, cellI) + { + blockUpper[cellI](dirI, dirJ) = upper[cellI]; + } + } + else + { + FatalErrorIn + ( + "void fvBlockMatrix::insertCouplingUpperLower\n" + "(\n" + " const direction dirI,\n" + " const direction dirJ,\n" + " const fvScalarMatrix& m\n" + ")" + ) << "Error in matrix insertion: problem with block structure." + << abort(FatalError); + } + + // Insert block interface fields + forAll(this->interfaces(), patchI) + { + if (this->interfaces().set(patchI)) + { + // Couple upper and lower + const scalarField& cUpper = matrix.boundaryCoeffs()[patchI]; + const scalarField& cLower = matrix.internalCoeffs()[patchI]; + + typename CoeffField::squareTypeField& blockUpper = + this->coupleUpper()[patchI].asSquare(); + + typename CoeffField::squareTypeField& blockLower = + this->coupleLower()[patchI].asSquare(); + + forAll(cUpper, faceI) + { + blockUpper[faceI](dirI, dirJ) = cUpper[faceI]; + blockLower[faceI](dirI, dirJ) = cLower[faceI]; + } + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::fvBlockMatrix::fvBlockMatrix +( + GeometricField& psi +) +: + BlockLduSystem(psi.mesh()), + psi_(psi) +{ + this->interfaces() = psi.boundaryField().blockInterfaces(); +} + + +template +Foam::fvBlockMatrix::fvBlockMatrix +( + const fvBlockMatrix& bxs +) +: + BlockLduSystem(bxs), + psi_(bxs.psi()) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +template +void fvBlockMatrix::retrieveSolution +( + const direction dir, + Field& xSingle +) const +{ + const direction nCmpts = pTraits::nComponents; + direction localDir = dir; + + const Field psiIn = psi_.internalField(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField xSingleCurr(xSingle.component(cmptI)); + + forAll (xSingleCurr, cellI) + { + xSingleCurr[cellI] = psiIn[cellI](localDir); + } + + xSingle.replace(cmptI, xSingleCurr); + + localDir++; + } +} + + +template +template +void fvBlockMatrix::insertEquation +( + const direction dir, + fvMatrix& matrix +) +{ + insertSolutionVector(dir, matrix.psi().internalField()); + insertDiagSource(dir, matrix); + insertUpperLower(dir, matrix); + updateCouplingCoeffs(dir, matrix); +} + + +template +template +void fvBlockMatrix::insertBlockCoupling +( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir +) +{ + insertBlock(dirI, dirJ, blockSystem, incrementColumnDir); + insertBoundaryContributions(dirI, dirJ, blockSystem, incrementColumnDir); +} + + +template +void fvBlockMatrix::insertEquationCoupling +( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix +) +{ + insertCouplingDiagSource(dirI, dirJ, matrix); + insertCouplingUpperLower(dirI, dirJ, matrix); +} + + +template +void fvBlockMatrix::blockAdd +( + const direction dir, + const scalarField& xSingle, + Field& blockX +) +{ + forAll(xSingle, cellI) + { + blockX[cellI](dir) += xSingle[cellI]; + } +} + + +template +void fvBlockMatrix::updateSourceCoupling() +{ + // Eliminate off-diagonal block coefficients from the square diagonal + // With this change, the segregated matrix can be assembled with complete + // source terms and linearisation can be provided independently. + // Once the linearisation coefficients are set (off-diagonal entries + // in the square block matrix, they are multiplied by the current value + // of the field and subtracted from the source term + + if (this->diag().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + + typename CoeffField::linearTypeField lf(blockDiag.size()); + typename CoeffField::squareTypeField sf(blockDiag.size()); + + // Expand and contract + + // Take out the diagonal entries from the square coefficient + contractLinear(lf, blockDiag); + + // Expand the diagonal for full square, with zeroes in the off-diagonal + expandLinear(sf, lf); + + // Subtract from the source the difference between the full block + // diagonal and the diagonal terms only + // Sign is the same as in the derivative + this->source() += (blockDiag - sf) & psi_.internalField(); + } +} + + +template +void fvBlockMatrix::insertPicardTensor +( + const direction UEqnDir, + const volVectorField& U, + const surfaceScalarField& phi +) +{ + const direction nCmpts = pTraits::nComponents; + direction localDirI = UEqnDir; + direction localDirJ = UEqnDir; + + // Sanity check + { + const direction blockMatrixSize = pTraits::nComponents; + + if (nCmpts > blockMatrixSize) + { + FatalErrorIn + ( + "void fvBlockMatrix::insertPicardTensor\n" + "(\n" + " const direction UEqnDir,\n" + " const volVectorField& U,\n" + " const surfaceScalarField& phi\n" + ")" + ) << "Trying to insert Picard tensor term into smaller " + << "fvBlockMatrix. Do you have momentum equation?" + << abort(FatalError); + } + } + + // Get weights for U which makes the implicit flux part + const fvMesh& mesh = U.mesh(); + + tmp > tinterpScheme = + fvc::scheme(mesh, U.name()); + + tmp tweights = tinterpScheme().weights(U); + const scalarField& wIn = tweights().internalField(); + + // Calculate the Pi tensor. Consider hard coding the interpolation scheme to + // correspond to the div(U, phi) interpolation scheme for consistency. + // VV, 21/July/2014. + const surfaceTensorField pi(mesh.Sf()*fvc::interpolate(U, phi, "(phi,U)")); +// const surfaceTensorField pi(fvc::interpolate(U, phi, "(phi,U)")*mesh.Sf()); + const tensorField& piIn = pi.internalField(); + + BlockLduSystem piSystem(mesh); + + // Get references to ldu fields of pi block system always as square + typename CoeffField::squareTypeField& piDiag = + piSystem.diag().asSquare(); + typename CoeffField::squareTypeField& piUpper = + piSystem.upper().asSquare(); + typename CoeffField::squareTypeField& piLower = + piSystem.lower().asSquare(); + + vectorField& piSource = piSystem.source(); + + piLower = -wIn*piIn; + piUpper = piLower + piIn; + piSystem.negSumDiag(); + + Info << "Diag coeff: " << piDiag[125] << nl + << "Lower coeff: " << piLower[125] << nl + << "Upper coeff: " << piUpper[125] << endl; + + // Boundary contributions - treated explicitly because of the problem with + // inconsistent return type of internalCoeffs. VV, 21/July/2014. + forAll(U.boundaryField(), patchI) + { + const fvPatchVectorField& Ub = U.boundaryField()[patchI]; + const fvPatch& patch = Ub.patch(); + const tensorField& pib = pi.boundaryField()[patchI]; + const fvsPatchScalarField& wb = tweights().boundaryField()[patchI]; + const unallocLabelList& fc = patch.faceCells(); + +// const vectorField internalCoeffs(Ub.valueInternalCoeffs(wb)); + + // Explicit diag boundary contribution +// forAll(Ub, faceI) +// { +// piSource[fc[faceI]] -= pib[faceI] & Ub[faceI]; +// piSource[fc[faceI]] -= Ub[faceI] & pib[faceI]; +// } + + // Hard coded zeroGradient if patch does not fix value + if (!U.boundaryField()[patchI].fixesValue()) + { + forAll(Ub, faceI) + { + piDiag[fc[faceI]] += pib[faceI]; + } + } + // Coupled patches treated implicitly + else if (patch.coupled()) + { + typename CoeffField::squareTypeField& pipCoupleUpper = + piSystem.coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pipCoupleLower = + piSystem.coupleLower()[patchI].asSquare(); + + const tensorField pcl = -wb*pib; + const tensorField pcu = pcl + pib; + + // Coupling contributions + pipCoupleLower -= pcl; + pipCoupleUpper -= pcu; + } + else + { + const vectorField boundaryCoeffs(Ub.valueBoundaryCoeffs(wb)); + + // Boundary contribution + forAll(Ub, faceI) + { + piSource[fc[faceI]] -= pib[faceI] & boundaryCoeffs[faceI]; + } + } + } + + // Consider chucking the above code into fvm::picardTensor operator. + // VV, 21/July/2014. + + // Finally add Picard piSystem into this fvBlockMatrix + typename CoeffField::squareTypeField& blockDiag = + this->diag().asSquare(); + typename CoeffField::squareTypeField& blockUpper = + this->upper().asSquare(); + typename CoeffField::squareTypeField& blockLower = + this->lower().asSquare(); + + Field& blockSource = this->source(); + + // Add diagonal, source, lower and upper + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) + { + forAll(piDiag, cellI) + { + blockDiag[cellI](localDirI, localDirJ) += + piDiag[cellI](cmptI, cmptJ); +// blockSource[cellI](localDirI, localDirJ) += +// piSource[cellI](cmptI, cmptJ); + } + + forAll(piUpper, faceI) + { + blockUpper[faceI](localDirI, localDirJ) += + piUpper[faceI](cmptI, cmptJ); + blockLower[faceI](localDirI, localDirJ) += + piLower[faceI](cmptI, cmptJ); + } + + localDirJ++; + } + + forAll(piSource, cellI) + { + blockSource[cellI](localDirI) += piSource[cellI](cmptI); + } + + localDirI++; + } + + // Reset local direction for coupling contributions + localDirI = UEqnDir; + localDirJ = UEqnDir; + + // Add coupling contributions + forAll(U.boundaryField(), patchI) + { + if (U.boundaryField()[patchI].patch().coupled()) + { + typename CoeffField::squareTypeField& pcoupleUpper = + this->coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + this->coupleLower()[patchI].asSquare(); + + const typename CoeffField::squareTypeField& pipcu = + piSystem.coupleUpper()[patchI].asSquare(); + const typename CoeffField::squareTypeField& pipcl = + piSystem.coupleLower()[patchI].asSquare(); + + for (direction cmptI = 0; cmptI < nCmpts; cmptI++) + { + for (direction cmptJ = 0; cmptJ < nCmpts; cmptJ++) + { + forAll(pipcu, faceI) + { + pcoupleUpper[faceI](localDirI, localDirJ) += + pipcu[faceI](cmptI, cmptJ); + pcoupleLower[faceI](localDirI, localDirJ) += + pipcl[faceI](cmptI, cmptJ); + } + } + } + + // Reset local directions for other patches + localDirI = UEqnDir; + localDirJ = UEqnDir; + } + } +} + + +template +Foam::BlockSolverPerformance fvBlockMatrix::solve +( + const dictionary& solverControls +) +{ + // Solver call + BlockSolverPerformance solverPerf = + BlockLduSolver::New + ( + psi_.name(), + *this, + solverControls + )->solve(psi_.internalField(), this->source()); + + // Print performance + solverPerf.print(); + + return solverPerf; +} + + +template +Foam::BlockSolverPerformance fvBlockMatrix::solve() +{ + return solve(psi_.mesh().solutionDict().solverDict(psi_.name())); +} + + +// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // + +template +Foam::Ostream& Foam::operator<< +( + Ostream& os, + const fvBlockMatrix& bxs +) +{ + os << static_cast&>(bxs) << nl + << bxs.psi() << endl; + + return os; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H new file mode 100644 index 000000000..d212a7253 --- /dev/null +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.H @@ -0,0 +1,302 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + Foam::fvBlockMatrix + +Description + fvBlockMatrix is an extension of fvMatrix for block coupled types. It holds + general insertion and retrieval tools for block coupling along with specific + functions for p-U coupling. + +Author + Vuko Vukcevic, FMENA Zagreb. + +SourceFiles + fvBlockMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fvBlockMatrix_H +#define fvBlockMatrix_H + +#include "BlockLduSystem.H" +#include "fvMatrices.H" +#include "blockLduSolvers.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of friend functions and operators + +template +class fvBlockMatrix; + +template +Ostream& operator<<(Ostream&, const fvBlockMatrix&); + + +/*---------------------------------------------------------------------------*\ + Class fvBlockMatrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class fvBlockMatrix +: + public BlockLduSystem +{ + // Private data + + //- Internal field reference + GeometricField& psi_; + + + // Private member functions + + // Insertion functions for fvMatrix into diagonal positions + + //- Insert internal field into this fvBlockMatrix + template + void insertSolutionVector + ( + const direction dir, + const Field& xSingle + ); + + //- Insert matrix diagonal and source into this fvBlockMatrix + template + void insertDiagSource + ( + const direction dir, + fvMatrix& matrix + ); + + //- Insert upper and lower part into this fvBlockMatrix + template + void insertUpperLower + ( + const direction dir, + const fvMatrix& matrix + ); + + //- Update coupling coefficients in this fvBlockMatrix + template + void updateCouplingCoeffs + ( + const direction dir, + const fvMatrix& matrix + ); + + + // Insertion functions for fvScalarMatrix into off-diagonal positions + // (coupling matrices) + // These two functions are dubious. Reconsider. HJ, 17/July/2014. + + //- Insert coupling matrix diag and source into this fvBlockMatrix + void insertCouplingDiagSource + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix + ); + + //- Insert coupling matrix lower and upper into this fvBlockMatrix + void insertCouplingUpperLower + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix + ); + + + // Pressure - velocity (p-U) coupling specific functions + + //- Insert BlockLduSystem (obtained by implicit grad/div operator) + // into this fvBlockMatrix + template + void insertBlock + ( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir + ); + + //- Insert source and coupling coeffs of BlockLduSystem (obtained by + // implicit grad/div operator) into this fvBlockMatrix + template + void insertBoundaryContributions + ( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir + ); + + +public: + + // Constructors + + //- Construct given a field to solve for + fvBlockMatrix(GeometricField&); + + //- Construct as copy + fvBlockMatrix(const fvBlockMatrix&); + + + // Destructor + + virtual ~fvBlockMatrix() + {} + + + // Member functions + + // Access + + //- Access to GeometricField + GeometricField& psi() + { + return psi_; + } + + //- const reference to GeometricField + const GeometricField& psi() const + { + return psi_; + } + + + // Insertion and retrieval public tools + + //- Retrieve part of internal (solved) field from this fvBlockMatrix + template + void retrieveSolution + ( + const direction dir, + Field& xSingle + ) const; + + //- Insert matrix into this fvBlockMatrix + template + void insertEquation + ( + const direction dir, + fvMatrix& matrix + ); + + //- Insert existing block system (obtained by implicit grad/div + // operator) into this fvBlockMatrix + template + void insertBlockCoupling + ( + const direction dirI, + const direction dirJ, + const BlockLduSystem& blockSystem, + const bool incrementColumnDir + ); + + //- Insert scalar equation coupling into this fvBlockMatrix + // Not tested: VV, 10/July/2014 + void insertEquationCoupling + ( + const direction dirI, + const direction dirJ, + const fvScalarMatrix& matrix + ); + + //- Add field into block field + void blockAdd + ( + const direction dir, + const scalarField& xSingle, + Field& blockX + ); + + //- Update coupling of block system. + // Subtracts the block-coefficient coupling as specified by the user + // from the source, leaving the implicit update given by + // linearisation + void updateSourceCoupling(); + + //- Insert Picard tensor term that comes from Picard linearisation + // of convection term in momentum equation. VV, 21/July/2014. + void insertPicardTensor + ( + const direction UEqnDir, + const volVectorField& U, + const surfaceScalarField& phi + ); + + + // Solver calls for fvBlockMatrix + + //- Solve returning the solution statistics. + // Use the given solver controls + BlockSolverPerformance solve(const dictionary&); + + //- Solve returning the solution statistics. + // Solver controls read from fvSolution + BlockSolverPerformance solve(); + + + // Member operators + + void operator=(const fvBlockMatrix&); + + void negate(); + + void operator+=(const fvBlockMatrix&); + void operator-=(const fvBlockMatrix&); + + void operator*=(const scalarField&); + void operator*=(const scalar); + + + // Ostream operator + + friend Ostream& operator<< + ( + Ostream&, + const fvBlockMatrix& + ); +}; + + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "fvBlockMatrix.C" +# include "fvBlockMatrixOperations.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrixOperations.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrixOperations.C new file mode 100644 index 000000000..e6d6c1a4b --- /dev/null +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrixOperations.C @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "fvBlockMatrix.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::fvBlockMatrix::negate() +{ + BlockLduSystem::negate(); + psi_.negate(); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +void Foam::fvBlockMatrix::operator= +( + const fvBlockMatrix& bxs +) +{ + if (this == &bxs) + { + FatalErrorIn + ( + "void fvBlockMatrix::operator=" + "(const fvBlockMatrix& bs)" + ) << "attempted assignment to self" + << abort(FatalError); + } + + BlockLduSystem::operator=(bxs); + psi_ = bxs.psi(); +} + + +template +void Foam::fvBlockMatrix::operator+= +( + const fvBlockMatrix& bxs +) +{ + BlockLduSystem::operator+=(bxs); + psi_ += bxs.psi(); +} + +template +void Foam::fvBlockMatrix::operator-= +( + const fvBlockMatrix& bxs +) +{ + BlockLduSystem::operator-=(bxs); + psi_ -= bxs.psi(); +} + +template +void Foam::fvBlockMatrix::operator*= +( + const scalarField& sf +) +{ + BlockLduSystem::operator*=(sf); + psi_ *= sf; +} + +template +void Foam::fvBlockMatrix::operator*= +( + const scalar s +) +{ + BlockLduSystem::operator*=(s); + psi_ *= s; +} + + +// ************************************************************************* // diff --git a/src/VectorN/finiteVolume/interpolation/VectorNSurfaceInterpolationSchemes.C b/src/finiteVolume/interpolation/surfaceInterpolation/VectorNSurfaceInterpolationSchemes/VectorNSurfaceInterpolationSchemes.C similarity index 100% rename from src/VectorN/finiteVolume/interpolation/VectorNSurfaceInterpolationSchemes.C rename to src/finiteVolume/interpolation/surfaceInterpolation/VectorNSurfaceInterpolationSchemes/VectorNSurfaceInterpolationSchemes.C diff --git a/src/foam/Make/files b/src/foam/Make/files index c436c0bdd..7a39691bd 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -212,6 +212,7 @@ dimensionedTypes/dimensionedDiagTensor/dimensionedDiagTensor.C dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.C dimensionedTypes/dimensionedSymmTensor4thOrder/dimensionedSymmTensor4thOrder.C dimensionedTypes/dimensionedTensor/dimensionedTensor.C +dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C matrices/solution/solution.C matrices/constraint/scalarConstraint.C diff --git a/src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C b/src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C similarity index 100% rename from src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C rename to src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C diff --git a/src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.H b/src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.H similarity index 100% rename from src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.H rename to src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.H diff --git a/src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorNI.H b/src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorNI.H similarity index 100% rename from src/VectorN/foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorNI.H rename to src/foam/dimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorNI.H diff --git a/src/VectorN/foam/GeometricFields/GeometricDiagTensorNFields.C b/src/foam/fields/GeometricFields/GeometricDiagTensorNFields/GeometricDiagTensorNFields.C similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricDiagTensorNFields.C rename to src/foam/fields/GeometricFields/GeometricDiagTensorNFields/GeometricDiagTensorNFields.C diff --git a/src/VectorN/foam/GeometricFields/GeometricDiagTensorNFields.H b/src/foam/fields/GeometricFields/GeometricDiagTensorNFields/GeometricDiagTensorNFields.H similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricDiagTensorNFields.H rename to src/foam/fields/GeometricFields/GeometricDiagTensorNFields/GeometricDiagTensorNFields.H diff --git a/src/VectorN/foam/GeometricFields/GeometricSphericalTensorNFields.C b/src/foam/fields/GeometricFields/GeometricSphericalTensorNFields/GeometricSphericalTensorNFields.C similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricSphericalTensorNFields.C rename to src/foam/fields/GeometricFields/GeometricSphericalTensorNFields/GeometricSphericalTensorNFields.C diff --git a/src/VectorN/foam/GeometricFields/GeometricSphericalTensorNFields.H b/src/foam/fields/GeometricFields/GeometricSphericalTensorNFields/GeometricSphericalTensorNFields.H similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricSphericalTensorNFields.H rename to src/foam/fields/GeometricFields/GeometricSphericalTensorNFields/GeometricSphericalTensorNFields.H diff --git a/src/VectorN/foam/GeometricFields/GeometricTensorNFields.C b/src/foam/fields/GeometricFields/GeometricTensorNFields/GeometricTensorNFields.C similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricTensorNFields.C rename to src/foam/fields/GeometricFields/GeometricTensorNFields/GeometricTensorNFields.C diff --git a/src/VectorN/foam/GeometricFields/GeometricTensorNFields.H b/src/foam/fields/GeometricFields/GeometricTensorNFields/GeometricTensorNFields.H similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricTensorNFields.H rename to src/foam/fields/GeometricFields/GeometricTensorNFields/GeometricTensorNFields.H diff --git a/src/VectorN/foam/GeometricFields/GeometricVectorNFields.C b/src/foam/fields/GeometricFields/GeometricVectorNFields/GeometricVectorNFields.C similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricVectorNFields.C rename to src/foam/fields/GeometricFields/GeometricVectorNFields/GeometricVectorNFields.C diff --git a/src/VectorN/foam/GeometricFields/GeometricVectorNFields.H b/src/foam/fields/GeometricFields/GeometricVectorNFields/GeometricVectorNFields.H similarity index 100% rename from src/VectorN/foam/GeometricFields/GeometricVectorNFields.H rename to src/foam/fields/GeometricFields/GeometricVectorNFields/GeometricVectorNFields.H diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockGAMGInterfaceField/processorBlockGAMGInterfaceField.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockGAMGInterfaceField/processorBlockGAMGInterfaceField.C index 41667d40e..9775a5e0b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockGAMGInterfaceField/processorBlockGAMGInterfaceField.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField/processorBlockGAMGInterfaceField/processorBlockGAMGInterfaceField.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "processorBlockGAMGInterfaceField.H" +#include "processorLduInterfaceField.H" #include "addToRunTimeSelectionTable.H" #include "lduMatrix.H" @@ -56,6 +57,14 @@ Foam::processorBlockGAMGInterfaceField::processorBlockGAMGInterfaceField doTransform_ = p.doTransform(); rank_ = p.rank(); } + else if (isA(fineInterfaceField)) + { + const processorLduInterfaceField& p = + refCast(fineInterfaceField); + + doTransform_ = p.doTransform(); + rank_ = p.rank(); + } else { FatalErrorIn("processorBlockGAMGInterfaceField Constructor")