From 8848498e5db59015d3ee66ce9ddfe08d5d3c9078 Mon Sep 17 00:00:00 2001 From: Dominik Christ Date: Fri, 23 May 2014 13:47:26 +0100 Subject: [PATCH] Added pUCoupledSolver with core library changes and tutorial case --- applications/solvers/coupled/Allwmake | 3 +- .../solvers/coupled/pUCoupledFoam/Make/files | 3 + .../coupled/pUCoupledFoam/Make/options | 15 + .../solvers/coupled/pUCoupledFoam/UEqn.H | 10 + .../pUCoupledFoam/calculateCouplingMatrices.H | 2 + .../coupled/pUCoupledFoam/createFields.H | 52 ++ .../pUCoupledFoam/initializeBlockMatrix.H | 9 + .../solvers/coupled/pUCoupledFoam/pEqn.H | 23 + .../coupled/pUCoupledFoam/pUCoupledFoam.C | 128 +++ .../pUCoupledFoam/readBlockSolverControls.H | 3 + .../blockMatrixTools/blockMatrixTools.C | 832 ++++++++++++++++++ .../blockMatrixTools/blockMatrixTools.H | 129 +++ .../divSchemes/divScheme/divScheme.H | 8 + .../gaussDivScheme/gaussDivScheme.C | 36 + .../gaussDivScheme/gaussDivScheme.H | 7 + src/finiteVolume/finiteVolume/fvm/fvm.H | 1 + src/finiteVolume/finiteVolume/fvm/fvmDiv.C | 30 + src/finiteVolume/finiteVolume/fvm/fvmDiv.H | 15 + src/finiteVolume/finiteVolume/fvm/fvmGrad.C | 79 ++ src/finiteVolume/finiteVolume/fvm/fvmGrad.H | 85 ++ .../gradSchemes/gaussGrad/gaussGrad.C | 36 + .../gradSchemes/gaussGrad/gaussGrad.H | 8 + .../gradSchemes/gradScheme/gradScheme.C | 27 + .../gradSchemes/gradScheme/gradScheme.H | 9 + .../incompressible/pUCoupledFoam/cavity/0/U | 41 + .../incompressible/pUCoupledFoam/cavity/0/p | 39 + .../pUCoupledFoam/cavity/Allclean | 11 + .../pUCoupledFoam/cavity/Allrun | 15 + .../cavity/constant/RASProperties | 200 +++++ .../cavity/constant/polyMesh/blockMeshDict | 63 ++ .../cavity/constant/polyMesh/boundary | 40 + .../cavity/constant/transportProperties | 21 + .../pUCoupledFoam/cavity/system/controlDict | 47 + .../pUCoupledFoam/cavity/system/fvSchemes | 59 ++ .../pUCoupledFoam/cavity/system/fvSolution | 38 + 35 files changed, 2123 insertions(+), 1 deletion(-) create mode 100644 applications/solvers/coupled/pUCoupledFoam/Make/files create mode 100644 applications/solvers/coupled/pUCoupledFoam/Make/options create mode 100644 applications/solvers/coupled/pUCoupledFoam/UEqn.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/createFields.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/pEqn.H create mode 100644 applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C create mode 100644 applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H create mode 100644 src/finiteVolume/finiteVolume/fvm/fvmGrad.C create mode 100644 src/finiteVolume/finiteVolume/fvm/fvmGrad.H create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/0/U create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/0/p create mode 100755 tutorials/incompressible/pUCoupledFoam/cavity/Allclean create mode 100755 tutorials/incompressible/pUCoupledFoam/cavity/Allrun create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/constant/RASProperties create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/blockMeshDict create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/constant/transportProperties create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/system/controlDict create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes create mode 100644 tutorials/incompressible/pUCoupledFoam/cavity/system/fvSolution diff --git a/applications/solvers/coupled/Allwmake b/applications/solvers/coupled/Allwmake index 3292954ac..eb198d84b 100755 --- a/applications/solvers/coupled/Allwmake +++ b/applications/solvers/coupled/Allwmake @@ -6,5 +6,6 @@ wmake libso conjugateHeatTransfer wmake blockCoupledScalarTransportFoam wmake conjugateHeatFoam -wmake conjugateHeatSimpleFoam +wmake conjugateHeatSimpleFoam +wmake pUCoupledFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/coupled/pUCoupledFoam/Make/files b/applications/solvers/coupled/pUCoupledFoam/Make/files new file mode 100644 index 000000000..c2200751c --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/Make/files @@ -0,0 +1,3 @@ +pUCoupledFoam.C + +EXE = $(FOAM_APPBIN)/pUCoupledFoam diff --git a/applications/solvers/coupled/pUCoupledFoam/Make/options b/applications/solvers/coupled/pUCoupledFoam/Make/options new file mode 100644 index 000000000..63db9d04b --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/Make/options @@ -0,0 +1,15 @@ +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 \ + -I$(LIB_SRC)/VectorN/lnInclude + +EXE_LIBS = \ + -lincompressibleRASModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -llduSolvers \ + -lVectorN + diff --git a/applications/solvers/coupled/pUCoupledFoam/UEqn.H b/applications/solvers/coupled/pUCoupledFoam/UEqn.H new file mode 100644 index 000000000..31693a9ab --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/UEqn.H @@ -0,0 +1,10 @@ +// Momentum equation +fvVectorMatrix UEqn +( + fvm::div(phi, U) + + turbulence->divDevReff(U) +); + +UEqn.relax(); + +blockMatrixTools::insertEquation(0, UEqn, A, x, b); diff --git a/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H b/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H new file mode 100644 index 000000000..a89c6bf2a --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H @@ -0,0 +1,2 @@ +blockVectorMatrix pInU(fvm::grad(p)); +blockVectorMatrix UInp(fvm::div(U)); diff --git a/applications/solvers/coupled/pUCoupledFoam/createFields.H b/applications/solvers/coupled/pUCoupledFoam/createFields.H new file mode 100644 index 000000000..1ea6ea880 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/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/pUCoupledFoam/initializeBlockMatrix.H b/applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H new file mode 100644 index 000000000..2e7d61e61 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H @@ -0,0 +1,9 @@ +// 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 new file mode 100644 index 000000000..03a32161f --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/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)) & mesh.Sf()) +); + +fvScalarMatrix pEqn +( + - fvm::laplacian(rUAf, p) + == + - fvc::div(presSource) +); + +pEqn.setReference(pRefCell, pRefValue); + +blockMatrixTools::insertEquation(3, pEqn, A, x, b); diff --git a/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C new file mode 100644 index 000000000..8aaef8d7a --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + pUCoupledFoam + +Description + Steady-state solver for incompressible, turbulent flow, with implicit + coupling between pressure and velocity achieved by BlockLduMatrix + Turbulence is in this version solved using the existing turbulence + structure. + +Authors + Klas Jareteg, Chalmers University of Technology, + Vuko Vukcevic, FMENA Zagreb. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "RASModel.H" + +#include "VectorNFieldTypes.H" +#include "volVectorNFields.H" +#include "blockLduSolvers.H" +#include "blockVectorNMatrices.H" + +#include "blockMatrixTools.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 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" + + Info<< "\nStarting time loop\n" << endl; + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + p.storePrevIter(); + + // Initialize block matrix +# include "initializeBlockMatrix.H" + + // Assemble and insert momentum equation +# include "UEqn.H" + + // Assemble and insert pressure equation +# include "pEqn.H" + + // Insert coupling, updating the boundary contributions + // Last argument in insertBlockCoupling says if the first location + // 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); + + // Solve the block matrix + BlockSolverPerformance solverPerf = + BlockLduSolver::New + ( + word("Up"), + A, + mesh.solutionDict().solver("Up") + )->solve(Up, b); + + solverPerf.print(); + + // Retrieve solution + blockMatrixTools::retrieveSolution(0, U.internalField(), Up); + blockMatrixTools::retrieveSolution(3, p.internalField(), Up); + + U.correctBoundaryConditions(); + p.correctBoundaryConditions(); + + phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource; + + p.relax(); + + turbulence->correct(); + runTime.write(); + +# include "continuityErrs.H" + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} diff --git a/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H b/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H new file mode 100644 index 000000000..4ee98e519 --- /dev/null +++ b/applications/solvers/coupled/pUCoupledFoam/readBlockSolverControls.H @@ -0,0 +1,3 @@ +label pRefCell = 0; +scalar pRefValue = 0; +setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue); diff --git a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C index 04f8f1ef1..a337af777 100644 --- a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C +++ b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C @@ -483,6 +483,838 @@ void updateSourceCoupling } } +template +void insertSolutionVector +( + const label loc, + const Field& xSingle, + Field& xBlock +) +{ + // Get number of field components and local copy of location, for + // consistency with member functions where locations need to be reset. + const label nCmpts = pTraits::nComponents; + label localLoc = loc; + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField xSingleCurr(xSingle.component(cmptI)); + + forAll (xSingleCurr, cellI) + { + xBlock[cellI](localLoc) = xSingleCurr[cellI]; + } + + localLoc++; + } +} + + +template +void retrieveSolution +( + const label loc, + Field& xSingle, + const Field& xBlock +) +{ + const label nCmpts = pTraits::nComponents; + label localLoc = loc; + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField xSingleCurr(xSingle.component(cmptI)); + + forAll (xSingleCurr, cellI) + { + xSingleCurr[cellI] = xBlock[cellI](localLoc); + } + + xSingle.replace(cmptI, xSingleCurr); + + localLoc++; + } +} + + +template +void insertDiagSource +( + const label loc, + fvMatrix& matrix, + BlockLduMatrix& A, + Field& b +) +{ + 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 label nCmpts = pTraits::nComponents; + label localLoc = loc; + + 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 + A.diag().activeType() != blockCoeffBase::SQUARE + ) + { + typename CoeffField::linearTypeField& blockDiag = + A.diag().asLinear(); + + for (label 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](localLoc) = diag[cellI]; + b[cellI](localLoc) += sourceCmpt[cellI]; + } + + localLoc++; + + // Reset diagonal + diag = saveDiag; + } + } + else if (A.diag().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockDiag = + A.diag().asSquare(); + + for (label 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](localLoc, localLoc) = diag[cellI]; + b[cellI](localLoc) += sourceCmpt[cellI]; + } + + localLoc++; + + // Reset diagonal + diag = saveDiag; + } + } +} + + +template +void insertUpperLower +( + const label loc, + const fvMatrix& matrix, + BlockLduMatrix& A +) +{ + if (matrix.diagonal()) + { + // Matrix for insertion is diagonal-only: nothing to do + return; + } + + const label nCmpts = pTraits::nComponents; + label localLoc = loc; + + if (matrix.hasUpper()) + { + const scalarField& upper = matrix.upper(); + + if (A.upper().activeType() == blockCoeffBase::UNALLOCATED) + { + A.upper().asScalar() = upper; + } + else if + ( + A.upper().activeType() == blockCoeffBase::SCALAR + || A.upper().activeType() == blockCoeffBase::LINEAR + ) + { + typename CoeffField::linearTypeField& blockUpper = + A.upper().asLinear(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (upper, faceI) + { + blockUpper[faceI](localLoc) = upper[faceI]; + } + + localLoc++; + } + } + else if (A.upper().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockUpper = + A.upper().asSquare(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (upper, faceI) + { + blockUpper[faceI](localLoc, localLoc) = upper[faceI]; + } + + localLoc++; + } + } + } + else + { + FatalErrorIn + ( + "void insertUpperLower\n" + "(\n" + " const location loc,\n" + " const fvMatrix& matrix,\n" + " BlockLduMatrix& A\n" + ")" + ) << "Error in matrix insertion: problem with block structure." + << abort(FatalError); + } + + if (matrix.symmetric() && A.symmetric()) + { + Info<< "Both matrices are symmetric: inserting only upper triangle" + << endl; + } + else + { + // Reset localLoc + localLoc = loc; + + // Either scalar or block matrix is asymmetric: insert lower triangle + const scalarField& lower = matrix.lower(); + + if (A.lower().activeType() == blockCoeffBase::UNALLOCATED) + { + A.lower().asScalar() = lower; + } + else if + ( + A.lower().activeType() == blockCoeffBase::SCALAR + || A.lower().activeType() == blockCoeffBase::LINEAR + ) + { + typename CoeffField::linearTypeField& blockLower = + A.lower().asLinear(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (lower, faceI) + { + blockLower[faceI](localLoc) = lower[faceI]; + } + + localLoc++; + } + } + else if (A.lower().activeType() == blockCoeffBase::SQUARE) + { + typename CoeffField::squareTypeField& blockLower = + A.lower().asSquare(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll (lower, faceI) + { + blockLower[faceI](localLoc, localLoc) = lower[faceI]; + } + + localLoc++; + } + } + } +} + + +template +void updateCouplingCoeffs +( + const label loc, + const fvMatrix& matrix, + BlockLduMatrix& A +) +{ + const label nCmpts = pTraits::nComponents; + label localLoc = loc; + + 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 (A.coupleUpper()[patchI].activeType() != blockCoeffBase::SQUARE) + { + typename CoeffField::linearTypeField& pcoupleUpper = + A.coupleUpper()[patchI].asLinear(); + typename CoeffField::linearTypeField& pcoupleLower = + A.coupleLower()[patchI].asLinear(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField icpCmpt = icp.component(cmptI); + scalarField bcpCmpt = bcp.component(cmptI); + + forAll(pf, faceI) + { + pcoupleUpper[faceI](localLoc) = bcpCmpt[faceI]; + pcoupleLower[faceI](localLoc) = icpCmpt[faceI]; + } + + localLoc++; + } + + localLoc = loc; + } + else if + ( + A.coupleUpper()[patchI].activeType() == blockCoeffBase::SQUARE + ) + { + typename CoeffField::squareTypeField& pcoupleUpper = + A.coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + A.coupleLower()[patchI].asSquare(); + + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + scalarField icpCmpt = icp.component(cmptI); + scalarField bcpCmpt = bcp.component(cmptI); + + forAll(pf, faceI) + { + pcoupleUpper[faceI](localLoc, localLoc) = + bcpCmpt[faceI]; + pcoupleLower[faceI](localLoc, localLoc) = + icpCmpt[faceI]; + } + + localLoc++; + } + + localLoc = loc; + } + } + } +} + + +template +void insertEquation +( + const label loc, + fvMatrix& matrix, + BlockLduMatrix& A, + Field& x, + Field& b +) +{ + insertDiagSource(loc, matrix, A, b); + insertUpperLower(loc, matrix, A); + updateCouplingCoeffs(loc, matrix, A); + insertSolutionVector(loc, matrix.psi().internalField(), x); +} + + +template +void insertBlock +( + const label loc1, + const label loc2, + const BlockLduMatrix& blockMatrix, + BlockLduMatrix& A, + const bool incFirst +) +{ + // Sanity checks + { + const label blockMatrixSize = pTraits::nComponents; + const label blockSystemSize = pTraits::nComponents; + + if (blockMatrixSize > blockSystemSize) + { + FatalErrorIn + ( + "void insertBlock\n" + "(\n" + " const label loc1,\n" + " const label loc2,\n" + " BlockLduMatrix& blockMatrix,\n" + " BlockLduMatrix& A\n" + ")" + ) << "Trying to insert a block matrix into smaller one." + << abort(FatalError); + } + + if (loc1 == loc2) + { + FatalErrorIn + ( + "void insertCoupling\n" + "(\n" + " const label loc1,\n" + " const label loc2,\n" + " BlockLduMatrix& blockMatrix,\n" + " BlockLduMatrix& A\n" + ")" + ) << "Trying to insert coupling in the position where equation " + << "should be, since loc1 = loc2. Try using insertEquatiion " + << "member function." + << abort(FatalError); + } + + } + + const label nCmpts = pTraits::nComponents; + label localLoc1 = loc1; + label localLoc2 = loc2; + + // Get references to ldu fields of blockMatrix always as linear + const typename CoeffField::linearTypeField& bmd = + blockMatrix.diag().asLinear(); + const typename CoeffField::linearTypeField& bmu = + blockMatrix.upper().asLinear(); + const typename CoeffField::linearTypeField& bml = + blockMatrix.lower().asLinear(); + + // Get references to ldu fields of A matrix always as square + typename CoeffField::squareTypeField& blockDiag = + A.diag().asSquare(); + typename CoeffField::squareTypeField& blockUpper = + A.upper().asSquare(); + typename CoeffField::squareTypeField& blockLower = + A.lower().asSquare(); + + // Insert blockMatrix that represents coupling into larger system matrix + for (label cmptI = 0; cmptI < nCmpts; cmptI++) + { + forAll(bmd, cellI) + { + blockDiag[cellI](localLoc1, localLoc2) += bmd[cellI](cmptI); + } + + forAll(bmu, faceI) + { + blockUpper[faceI](localLoc1, localLoc2) += bmu[faceI](cmptI); + blockLower[faceI](localLoc1, localLoc2) += bml[faceI](cmptI); + } + + if (incFirst) + { + localLoc1++; + } + else + { + localLoc2++; + } + } +} + + +template +void insertBoundaryContributions +( + const label loc1, + const label loc2, + const BlockLduMatrix& blockMatrix, + GeometricField& psi, + 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)); + + 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) + { + const fvPatchScalarField& pf = psi.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + const vectorField& Sf = patch.Sf(); + const fvsPatchScalarField& pw = weights.boundaryField()[patchI]; + + scalarField internalCoeffs(pf.valueInternalCoeffs(pw)); + + if (patch.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; + + 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; + } + 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; + } + } +} + + +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) + { + 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()) + { + typename CoeffField::squareTypeField& pcoupleUpper = + A.coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + A.coupleLower()[patchI].asSquare(); + + 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) += 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++; + } + } + + // 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; + } + } +} + + +template +void insertBlockCoupling +( + const label loc1, + const label loc2, + const BlockLduMatrix& blockMatrix, + GeometricField& psi, + BlockLduMatrix& A, + Field& b, + const bool incFirst +) +{ + insertBlock(loc1, loc2, blockMatrix, A, incFirst); + insertBoundaryContributions(loc1, loc2, blockMatrix, psi, A, b, incFirst); +} + + +template +void insertEquationCoupling +( + const label loc1, + const label loc2, + fvScalarMatrix& matrix, + BlockLduMatrix& A, + Field& b +) +{ + // Sanity check + if (loc1 == loc2) + { + FatalErrorIn + ( + "void insertEquationCoupling\n" + "(\n" + " const label loc1,\n" + " const label loc2,\n" + " fvScalarMatrix& matrix\n" + " BlockLduMatrix& A\n" + " Field& b\n" + ")" + ) << "Trying to insert coupling in the position where equation " + << "should be since loc1 = loc2. Try using insertEquatiion " + << "member function." + << abort(FatalError); + } + + // Get references to fvScalarMatrix fields, updating boundary contributions + scalarField& diag = matrix.D(); + scalarField& source = matrix.source(); + matrix.addBoundarySource(source, false); + + const scalarField& upper = matrix.upper(); + const scalarField& lower = matrix.lower(); + + // Get references to ldu fields of A matrix always as square + typename CoeffField::squareTypeField& blockDiag = + A.diag().asSquare(); + typename CoeffField::squareTypeField& blockUpper = + A.upper().asSquare(); + typename CoeffField::squareTypeField& blockLower = + A.lower().asSquare(); + + forAll(diag, cellI) + { + blockDiag[cellI](loc1, loc2) += diag[cellI]; + b[cellI](loc1) += source[cellI]; + } + + forAll(upper, faceI) + { + blockUpper[faceI](loc1, loc2) += upper[faceI]; + blockLower[faceI](loc1, loc2) += lower[faceI]; + } + + // Update coupling contributions + const volScalarField& psi = matrix.psi(); + forAll(psi.boundaryField(), patchI) + { + const fvPatchScalarField& pf = psi.boundaryField()[patchI]; + const fvPatch& patch = pf.patch(); + + if (patch.coupled()) + { + const scalarField& icp = matrix.internalCoeffs()[patchI]; + const scalarField& bcp = matrix.boundaryCoeffs()[patchI]; + + typename CoeffField::squareTypeField& pcoupleUpper = + A.coupleUpper()[patchI].asSquare(); + typename CoeffField::squareTypeField& pcoupleLower = + A.coupleLower()[patchI].asSquare(); + + forAll(pf, faceI) + { + pcoupleUpper[faceI](loc1, loc2) = bcp[faceI]; + pcoupleLower[faceI](loc1, loc2) = icp[faceI]; + } + } + } +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H index a53b099f3..7000c0b4e 100644 --- a/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H +++ b/src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H @@ -158,6 +158,135 @@ namespace blockMatrixTools Field& x, Field& b ); + + //- Tools specific to pU coupled block matrix + + //- Insert field into block field to be solved for + template + void insertSolutionVector + ( + const label loc, + const Field& xSingle, + Field& xBlock + ); + + //- Retrieve solved field from block field + template + void retrieveSolution + ( + const label loc, + Field& xSingle, + const Field& xBlock + ); + + //- Insert matrix diagonal and source into the block system + template + void insertDiagSource + ( + const label loc, + fvMatrix& matrix, + BlockLduMatrix& A, + Field& b + ); + + // Insert upper and lower part into the block system + template + void insertUpperLower + ( + const label loc, + const fvMatrix& matrix, + BlockLduMatrix& A + ); + + // Update coupling coefficients in the block matrix + template + void updateCouplingCoeffs + ( + const label loc, + const fvMatrix& matrix, + BlockLduMatrix& A + ); + + // Insert matrix into the block system + template + void insertEquation + ( + const label loc, + fvMatrix& matrix, + BlockLduMatrix& A, + Field& x, + Field& b + ); + + // Insert block coupling + template + void insertBlock + ( + const label loc1, + const label loc2, + const BlockLduMatrix& 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 + void insertBoundaryContributions + ( + const label loc1, + const label loc2, + const BlockLduMatrix& blockMatrix, + GeometricField& psi, + 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. + template + void insertBlockCoupling + ( + const label loc1, + const label loc2, + const BlockLduMatrix& blockMatrix, + GeometricField& psi, + BlockLduMatrix& A, + Field& b, + const bool incFirst + ); + + // Insert scalar equation coupling. Not tested: VV, 9/May/2014 + template + void insertEquationCoupling + ( + const label loc1, + const label loc2, + fvScalarMatrix& matrix, + BlockLduMatrix& A, + Field b + ); } diff --git a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H index dd5cf5d1b..de18f5763 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H @@ -41,6 +41,7 @@ SourceFiles #include "linear.H" #include "typeInfo.H" #include "runTimeSelectionTables.H" +#include "blockLduMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -150,6 +151,13 @@ public: ( const GeometricField& ) = 0; + + //- Return the blockLduMatrix corresponding to the implicit + // div discretization. For block coupled system. + virtual tmp fvmDiv + ( + const GeometricField& + ) const = 0; }; diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C index 13910e25a..1457139fb 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C @@ -68,6 +68,42 @@ gaussDivScheme::fvcDiv } +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 tbm + ( + new blockVectorMatrix + ( + 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; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H index 2c22a6494..8ebd4e227 100644 --- a/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H +++ b/src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H @@ -96,6 +96,13 @@ public: ( const GeometricField& ); + + //- Return the blockLduMatrix corresponding to the implicit + // div discretization. For block coupled system. + tmp fvmDiv + ( + const GeometricField& + ) const; }; diff --git a/src/finiteVolume/finiteVolume/fvm/fvm.H b/src/finiteVolume/finiteVolume/fvm/fvm.H index 7a140afe9..47b8c0f68 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvm.H +++ b/src/finiteVolume/finiteVolume/fvm/fvm.H @@ -43,6 +43,7 @@ Description #include "fvmDdt.H" #include "fvmD2dt2.H" #include "fvmDiv.H" +#include "fvmGrad.H" #include "fvmAdjDiv.H" #include "fvmLaplacian.H" #include "fvmSup.H" diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C index 5be5d778c..bd6ea7772 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C @@ -26,6 +26,7 @@ License #include "fvmDiv.H" #include "fvMesh.H" #include "convectionScheme.H" +#include "divScheme.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -83,6 +84,7 @@ div return fvm::div(flux, vf, "div("+flux.name()+','+vf.name()+')'); } + template tmp > div @@ -96,6 +98,34 @@ div return Div; } +template +tmp > div +( + GeometricField& vf, + const word& name +) +{ + return fv::divScheme::New + ( + vf.mesh(), + vf.mesh().schemesDict().divScheme(name) + )().fvmDiv(vf); +} + + +template +tmp > div +( + GeometricField& vf +) +{ + return fvm::div + ( + vf, + "div(" + vf.name() + ')' + ); +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H index 6b23c3389..a20bff42d 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H @@ -39,6 +39,7 @@ SourceFiles #include "surfaceFieldsFwd.H" #include "surfaceInterpolationScheme.H" #include "fvMatrices.H" +#include "blockLduMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -81,6 +82,20 @@ namespace fvm const tmp&, GeometricField& ); + + // Implicit div operator for block systems + template + tmp div + ( + GeometricField&, + const word& + ); + + template + tmp div + ( + GeometricField& + ); } diff --git a/src/finiteVolume/finiteVolume/fvm/fvmGrad.C b/src/finiteVolume/finiteVolume/fvm/fvmGrad.C new file mode 100644 index 000000000..b17168368 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvm/fvmGrad.C @@ -0,0 +1,79 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +\*---------------------------------------------------------------------------*/ + +#include "fvmGrad.H" +#include "gradScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fvm +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp grad +( + GeometricField& vf, + const word& name +) +{ + return fv::gradScheme::New + ( + vf.mesh(), + vf.mesh().schemesDict().gradScheme(name) + )().fvmGrad(vf); +} + + +template +tmp grad +( + GeometricField& vf +) +{ + return fvm::grad + ( + vf, + "grad(" + vf.name() + ')' + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fvm + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvm/fvmGrad.H b/src/finiteVolume/finiteVolume/fvm/fvmGrad.H new file mode 100644 index 000000000..52ff1ab62 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvm/fvmGrad.H @@ -0,0 +1,85 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +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. + +SourceFiles + fvmGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fvmGrad_H +#define fvmGrad_H + +#include "volFieldsFwd.H" +#include "blockLduMatrices.H" +#include "geometricOneField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace fvm functions Declaration +\*---------------------------------------------------------------------------*/ + +namespace fvm +{ + template + tmp grad + ( + GeometricField&, + const word& + ); + + template + tmp grad + ( + GeometricField& + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "fvmGrad.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C index aa926b0f7..c27e5fd0a 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C @@ -145,6 +145,42 @@ gaussGrad::grad } +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 tbm + ( + new blockVectorMatrix + ( + 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; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace fv diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H index 13253ce55..3010011fb 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H @@ -139,6 +139,14 @@ public: const GeometricField& ) const; + //- Return the blockLduMatrix corresponding to the implicit + // grad discretization. For block coupled systems. + tmp 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 1237a2695..bb273364b 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C @@ -98,6 +98,33 @@ gradScheme::~gradScheme() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +tmp gradScheme::fvmGrad +( + const GeometricField& vf +) const +{ + FatalIOErrorIn + ( + "tmp fvmGrad\n", + "(\n" + " GeometricField&" + ")\n" + ) << "Implicit gradient operator currently defined only for Gauss grad." + << abort(FatalIOError); + + tmp tbm + ( + new blockVectorMatrix + ( + vf.mesh() + ) + ); + + return tbm; +} + + template void gradScheme::correctBoundaryConditions ( diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H index 959cfe129..b7cf31e66 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H +++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H @@ -40,6 +40,7 @@ SourceFiles #include "surfaceFieldsFwd.H" #include "typeInfo.H" #include "runTimeSelectionTables.H" +#include "blockLduMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -137,6 +138,14 @@ public: ) const = 0; + //- Return the blockLduMatrix corresponding to the implicit + // grad discretization. For block coupled systems. + virtual tmp fvmGrad + ( + const GeometricField& + ) const; + + // Moved from gaussGrad into base class. HJ, 14/Jun/2013 //- Correct the boundary values of the gradient using the patchField // snGrad functions diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/0/U b/tutorials/incompressible/pUCoupledFoam/cavity/0/U new file mode 100644 index 000000000..385eb2b65 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/0/U @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + movingWall + { + type fixedValue; + value uniform (1 0 0); + } + + fixedWalls + { + type fixedValue; + value uniform (0 0 0); + } + + frontAndBack + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/0/p b/tutorials/incompressible/pUCoupledFoam/cavity/0/p new file mode 100644 index 000000000..24af6884d --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/0/p @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + movingWall + { + type zeroGradient; + } + + fixedWalls + { + type zeroGradient; + } + + frontAndBack + { + type empty; + } +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/Allclean b/tutorials/incompressible/pUCoupledFoam/cavity/Allclean new file mode 100755 index 000000000..0627e8600 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/Allclean @@ -0,0 +1,11 @@ +#!/bin/bash + +# Klas Jareteg, 2012-10-13 +# Description: +# Script to run comparative case between coupled and segregated solvers. + + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +cleanCase diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/Allrun b/tutorials/incompressible/pUCoupledFoam/cavity/Allrun new file mode 100755 index 000000000..0f4b2840c --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/Allrun @@ -0,0 +1,15 @@ +#!/bin/bash + +# Klas Jareteg, 2012-10-13 +# Description: +# Script to run comparative case between coupled and segregated solvers. + + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +application=pUCoupledFoam + +runApplication blockMesh +runApplication $application + diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/constant/RASProperties b/tutorials/incompressible/pUCoupledFoam/cavity/constant/RASProperties new file mode 100644 index 000000000..626f9da45 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/constant/RASProperties @@ -0,0 +1,200 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object RASProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +RASModel laminar; + +turbulence on; + +printCoeffs on; + +laminarCoeffs +{ +} + +kEpsilonCoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphaEps 0.76923; +} + +RNGkEpsilonCoeffs +{ + Cmu 0.0845; + C1 1.42; + C2 1.68; + alphak 1.39; + alphaEps 1.39; + eta0 4.38; + beta 0.012; +} + +realizableKECoeffs +{ + Cmu 0.09; + A0 4.0; + C2 1.9; + alphak 1; + alphaEps 0.833333; +} + +kOmegaSSTCoeffs +{ + alphaK1 0.85034; + alphaK2 1.0; + alphaOmega1 0.5; + alphaOmega2 0.85616; + gamma1 0.5532; + gamma2 0.4403; + beta1 0.0750; + beta2 0.0828; + betaStar 0.09; + a1 0.31; + c1 10; + + Cmu 0.09; +} + +NonlinearKEShihCoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphak 1; + alphaEps 0.76932; + A1 1.25; + A2 1000; + Ctau1 -4; + Ctau2 13; + Ctau3 -2; + alphaKsi 0.9; +} + +LienCubicKECoeffs +{ + C1 1.44; + C2 1.92; + alphak 1; + alphaEps 0.76923; + A1 1.25; + A2 1000; + Ctau1 -4; + Ctau2 13; + Ctau3 -2; + alphaKsi 0.9; +} + +QZetaCoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphaZeta 0.76923; + anisotropic no; +} + +LaunderSharmaKECoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphaEps 0.76923; +} + +LamBremhorstKECoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphaEps 0.76923; +} + +LienCubicKELowReCoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphak 1; + alphaEps 0.76923; + A1 1.25; + A2 1000; + Ctau1 -4; + Ctau2 13; + Ctau3 -2; + alphaKsi 0.9; + Am 0.016; + Aepsilon 0.263; + Amu 0.00222; +} + +LienLeschzinerLowReCoeffs +{ + Cmu 0.09; + C1 1.44; + C2 1.92; + alphak 1; + alphaEps 0.76923; + Am 0.016; + Aepsilon 0.263; + Amu 0.00222; +} + +LRRCoeffs +{ + Cmu 0.09; + Clrr1 1.8; + Clrr2 0.6; + C1 1.44; + C2 1.92; + Cs 0.25; + Ceps 0.15; + alphaEps 0.76923; +} + +LaunderGibsonRSTMCoeffs +{ + Cmu 0.09; + Clg1 1.8; + Clg2 0.6; + C1 1.44; + C2 1.92; + C1Ref 0.5; + C2Ref 0.3; + Cs 0.25; + Ceps 0.15; + alphaEps 0.76923; + alphaR 1.22; +} + +SpalartAllmarasCoeffs +{ + alphaNut 1.5; + Cb1 0.1355; + Cb2 0.622; + Cw2 0.3; + Cw3 2; + Cv1 7.1; + Cv2 5.0; +} + +wallFunctionCoeffs +{ + kappa 0.4187; + E 9; +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/blockMeshDict b/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/blockMeshDict new file mode 100644 index 000000000..4656ad967 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/blockMeshDict @@ -0,0 +1,63 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 0.1; + +vertices +( + (0 0 0) + (1 0 0) + (1 1 0) + (0 1 0) + (0 0 0.1) + (1 0 0.1) + (1 1 0.1) + (0 1 0.1) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (50 50 1) simpleGrading (1 1 1) +); + +edges +( +); + +patches +( + wall movingWall + ( + (3 7 6 2) + ) + wall fixedWalls + ( + (0 4 7 3) + (2 6 5 1) + (1 5 4 0) + ) + empty frontAndBack + ( + (0 3 2 1) + (4 5 6 7) + ) +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary b/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary new file mode 100644 index 000000000..dd979bf9a --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/constant/polyMesh/boundary @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | foam-extend: Open Source CFD | +| \\ / O peration | Version: 3.0 | +| \\ / A nd | Web: http://www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +3 +( + movingWall + { + type wall; + nFaces 50; + startFace 4900; + } + fixedWalls + { + type wall; + nFaces 150; + startFace 4950; + } + frontAndBack + { + type empty; + nFaces 5000; + startFace 5100; + } +) + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/constant/transportProperties b/tutorials/incompressible/pUCoupledFoam/cavity/constant/transportProperties new file mode 100644 index 000000000..46022c2b4 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/constant/transportProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +transportModel Newtonian; + +nu nu [0 2 -1 0 0 0 0] 0.01; + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/system/controlDict b/tutorials/incompressible/pUCoupledFoam/cavity/system/controlDict new file mode 100644 index 000000000..a82ec943e --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/system/controlDict @@ -0,0 +1,47 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application pUCoupledFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 0.5; + +deltaT 0.005; + +writeControl timeStep; + +writeInterval 20; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression uncompressed; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes new file mode 100644 index 000000000..136a79ddc --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSchemes @@ -0,0 +1,59 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default steadyState; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default none; + div(phi,U) Gauss linearUpwind Gauss linear; + div((nuEff*dev(grad(U).T()))) Gauss linear; + + div(U) Gauss linear; +} + +laplacianSchemes +{ + default none; + laplacian(nuEff,U) Gauss linear corrected; + laplacian(rUAf,p) Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} + +fluxRequired +{ + default no; + p; +} + +// ************************************************************************* // diff --git a/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSolution b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSolution new file mode 100644 index 000000000..d75f853b4 --- /dev/null +++ b/tutorials/incompressible/pUCoupledFoam/cavity/system/fvSolution @@ -0,0 +1,38 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM Extend Project: Open Source CFD | +| \\ / O peration | Version: 1.6-ext | +| \\ / A nd | Web: www.extend-project.de | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + Up + { + solver BiCGStab; + preconditioner Cholesky; + + tolerance 1e-09; + relTol 0.0; + + minIter 1; + maxIter 500; + } +} + +blockSolver +{ + pRefCell 0; + pRefValue 0; +} + +// ************************************************************************* //