1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
26 blockCoupledScalarTransportFoam
29 Solves two coupled transport equations in a block-coupled manner
31 1) transport equation for a passive scalar
34 This resembles heat exchanging flow through a porous medium
36 \*---------------------------------------------------------------------------*/
39 #include "fieldTypes.H"
40 #include "blockLduMatrices.H"
41 #include "blockLduSolvers.H"
45 #include "blockVector2Matrix.H"
47 #include "vector2Field.H"
48 #include "tensor2Field.H"
49 #include "blockMatrixTools.H"
51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53 int main(int argc, char *argv[])
56 # include "setRootCase.H"
57 # include "createTime.H"
58 # include "createMesh.H"
59 # include "createFields.H"
61 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
63 Info<< "\nCalculating scalar transport\n" << endl;
65 # include "CourantNo.H"
67 for (runTime++; !runTime.end(); runTime++)
69 Info<< "Time = " << runTime.timeName() << nl << endl;
71 # include "readSIMPLEControls.H"
73 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
78 - fvm::laplacian(DT, T)
88 - fvm::laplacian(DTs, Ts)
96 // Prepare block system
97 BlockLduMatrix<vector2> blockM(mesh);
99 // Grab block diagonal and set it to zero
100 Field<tensor2>& d = blockM.diag().asSquare();
103 // Grab linear off-diagonal and set it to zero
104 Field<vector2>& u = blockM.upper().asLinear();
105 Field<vector2>& l = blockM.lower().asLinear();
109 vector2Field blockX(mesh.nCells(), vector2::zero);
110 vector2Field blockB(mesh.nCells(), vector2::zero);
112 //- Inset equations into block Matrix
113 blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB);
114 blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB);
116 //- Add off-diagonal terms and remove from Block source
119 d[i](0,1) = -alpha.value()*mesh.V()[i];
120 d[i](1,0) = -alpha.value()*mesh.V()[i];
122 blockB[i][0] -= alpha.value()*blockX[i][1]*mesh.V()[i];
123 blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i];
127 BlockSolverPerformance<vector2> solverPerf =
128 BlockLduSolver<vector2>::New
132 mesh.solver("blockVar")
133 )->solve(blockX, blockB);
138 blockMatrixTools::blockRetrieve(0, T.internalField(), blockX);
139 blockMatrixTools::blockRetrieve(1, Ts.internalField(), blockX);
141 T.correctBoundaryConditions();
142 Ts.correctBoundaryConditions();
148 Info<< "End\n" << endl;
154 // ************************************************************************* //