src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  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.
  15
  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
  19     for more details.
  20
  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
  24
  25 Description
  26      Finite-Area scalar matrix member functions and operators
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "faScalarMatrix.H"
  31 #include "zeroGradientFaPatchFields.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  39
  40 // Set reference level for a component of the solution
  41 // on a given patch face
  42 template<>
  43 void faMatrix<scalar>::setComponentReference
  44 (
  45     const label patchI,
  46     const label edgeI,
  47     const direction,
  48     const scalar value
  49 )
  50 {
  51     const unallocLabelList& faceLabels =
  52         psi_.mesh().boundary()[patchI].edgeFaces();
  53
  54     internalCoeffs_[patchI][edgeI] +=
  55         diag()[faceLabels[edgeI]];
  56
  57     boundaryCoeffs_[patchI][edgeI] = value;
  58 }
  59
  60
  61 template<>
  62 autoPtr<faMatrix<scalar>::faSolver> faMatrix<scalar>::solver
  63 (
  64     const dictionary& solverControls
  65 )
  66 {
  67     if (debug)
  68     {
  69         Info<< "faMatrix<scalar>::solver(const dictionary&) : "
  70                "solver for faMatrix<scalar>"
  71             << endl;
  72     }
  73
  74     scalarField saveDiag = diag();
  75     addBoundaryDiag(diag(), 0);
  76
  77     // Make a copy of interfaces: no longer a reference
  78     // HJ, 20/Nov/2007
  79     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  80
  81     autoPtr<faMatrix<scalar>::faSolver> solverPtr
  82     (
  83         new faMatrix<scalar>::faSolver
  84         (
  85             *this,
  86             lduSolver::New
  87             (
  88                 psi_.name(),
  89                 *this,
  90                 boundaryCoeffs_,
  91                 internalCoeffs_,
  92                 interfaces,
  93                 solverControls
  94             )
  95         )
  96     );
  97
  98     diag() = saveDiag;
  99
 100     return solverPtr;
 101 }
 102
 103
 104 template<>
 105 lduSolverPerformance faMatrix<scalar>::faSolver::solve
 106 (
 107     const dictionary& solverControls
 108 )
 109 {
 110     scalarField saveDiag = faMat_.diag();
 111     faMat_.addBoundaryDiag(faMat_.diag(), 0);
 112
 113     scalarField totalSource = faMat_.source();
 114     faMat_.addBoundarySource(totalSource, false);
 115
 116     solver_->read(solverControls);
 117     lduSolverPerformance solverPerf =
 118         solver_->solve(faMat_.psi().internalField(), totalSource);
 119
 120     solverPerf.print();
 121
 122     faMat_.diag() = saveDiag;
 123
 124     faMat_.psi().correctBoundaryConditions();
 125
 126     return solverPerf;
 127 }
 128
 129
 130 template<>
 131 lduSolverPerformance faMatrix<scalar>::solve
 132 (
 133     const dictionary& solverControls
 134 )
 135 {
 136     if (debug)
 137     {
 138         Info<< "faMatrix<scalar>::solve(const dictionary&) : "
 139                "solving faMatrix<scalar>"
 140             << endl;
 141     }
 142
 143     scalarField saveDiag = diag();
 144     addBoundaryDiag(diag(), 0);
 145
 146     scalarField totalSource = source_;
 147     addBoundarySource(totalSource, 0);
 148
 149     // Make a copy of interfaces: no longer a reference
 150     // HJ, 20/Nov/2007
 151     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 152
 153     // Solver call
 154     lduSolverPerformance solverPerf = lduSolver::New
 155     (
 156         psi_.name(),
 157         *this,
 158         boundaryCoeffs_,
 159         internalCoeffs_,
 160         interfaces,
 161         solverControls
 162     )->solve(psi_.internalField(), totalSource);
 163
 164     solverPerf.print();
 165
 166     diag() = saveDiag;
 167
 168     psi_.correctBoundaryConditions();
 169
 170     return solverPerf;
 171 }
 172
 173
 174 // Return the matrix residual
 175 template<>
 176 tmp<scalarField> faMatrix<scalar>::residual() const
 177 {
 178     scalarField boundaryDiag(psi_.size(), 0.0);
 179     addBoundaryDiag(boundaryDiag, 0);
 180
 181     // Make a copy of interfaces: no longer a reference
 182     // HJ, 20/Nov/2007
 183     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 184
 185     tmp<scalarField> tres
 186     (
 187         lduMatrix::residual
 188         (
 189             psi_.internalField(),
 190             source_ - boundaryDiag*psi_.internalField(),
 191             boundaryCoeffs_,
 192             interfaces,
 193             0
 194         )
 195     );
 196
 197     addBoundarySource(tres());
 198
 199     return tres;
 200 }
 201
 202
 203 // H operator
 204 template<>
 205 tmp<areaScalarField> faMatrix<scalar>::H() const
 206 {
 207     tmp<areaScalarField> tHphi
 208     (
 209         new areaScalarField
 210         (
 211             IOobject
 212             (
 213                 "H("+psi_.name()+')',
 214                 psi_.instance(),
 215                 psi_.db(),
 216                 IOobject::NO_READ,
 217                 IOobject::NO_WRITE
 218             ),
 219             psi_.mesh(),
 220             dimensions_/dimArea,
 221             zeroGradientFaPatchScalarField::typeName
 222         )
 223     );
 224     areaScalarField Hphi = tHphi();
 225
 226     Hphi.internalField() = (lduMatrix::H(psi_.internalField()) + source_);
 227     addBoundarySource(Hphi.internalField());
 228
 229     Hphi.internalField() /= psi_.mesh().S();
 230     Hphi.correctBoundaryConditions();
 231
 232     return tHphi;
 233 }
 234
 235
 236 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 237
 238 } // End namespace Foam
 239
 240 // ************************************************************************* //