src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 \*---------------------------------------------------------------------------*/
  26
  27 #include "fvScalarMatrix.H"
  28 #include "zeroGradientFvPatchFields.H"
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 template<>
  33 void Foam::fvMatrix<Foam::scalar>::setComponentReference
  34 (
  35     const label patchi,
  36     const label facei,
  37     const direction,
  38     const scalar value
  39 )
  40 {
  41     if (psi_.needReference())
  42     {
  43         if (Pstream::master())
  44         {
  45             internalCoeffs_[patchi][facei] +=
  46                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
  47
  48             boundaryCoeffs_[patchi][facei] +=
  49                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
  50                *value;
  51         }
  52     }
  53 }
  54
  55
  56 template<>
  57 Foam::autoPtr<Foam::fvMatrix<Foam::scalar>::fvSolver>
  58 Foam::fvMatrix<Foam::scalar>::solver
  59 (
  60     const dictionary& solverControls
  61 )
  62 {
  63     if (debug)
  64     {
  65         Info<< "fvMatrix<scalar>::solver(const dictionary& solverControls) : "
  66                "solver for fvMatrix<scalar>"
  67             << endl;
  68     }
  69
  70     scalarField saveDiag = diag();
  71     addBoundaryDiag(diag(), 0);
  72
  73     autoPtr<fvMatrix<scalar>::fvSolver> solverPtr
  74     (
  75         new fvMatrix<scalar>::fvSolver
  76         (
  77             *this,
  78             lduMatrix::solver::New
  79             (
  80                 psi_.name(),
  81                 *this,
  82                 boundaryCoeffs_,
  83                 internalCoeffs_,
  84                 psi_.boundaryField().interfaces(),
  85                 solverControls
  86             )
  87         )
  88     );
  89
  90     diag() = saveDiag;
  91
  92     return solverPtr;
  93 }
  94
  95
  96 template<>
  97 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve
  98 (
  99     const dictionary& solverControls
 100 )
 101 {
 102     scalarField saveDiag = fvMat_.diag();
 103     fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
 104
 105     scalarField totalSource = fvMat_.source();
 106     fvMat_.addBoundarySource(totalSource, false);
 107
 108     // assign new solver controls
 109     solver_->read(solverControls);
 110
 111     lduMatrix::solverPerformance solverPerf =
 112         solver_->solve(fvMat_.psi().internalField(), totalSource);
 113
 114     solverPerf.print();
 115
 116     fvMat_.diag() = saveDiag;
 117
 118     fvMat_.psi().correctBoundaryConditions();
 119
 120     return solverPerf;
 121 }
 122
 123
 124 template<>
 125 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
 126 (
 127     const dictionary& solverControls
 128 )
 129 {
 130     if (debug)
 131     {
 132         Info<< "fvMatrix<scalar>::solve(const dictionary& solverControls) : "
 133                "solving fvMatrix<scalar>"
 134             << endl;
 135     }
 136
 137     scalarField saveDiag = diag();
 138     addBoundaryDiag(diag(), 0);
 139
 140     scalarField totalSource = source_;
 141     addBoundarySource(totalSource, false);
 142
 143     // Solver call
 144     lduMatrix::solverPerformance solverPerf = lduMatrix::solver::New
 145     (
 146         psi_.name(),
 147         *this,
 148         boundaryCoeffs_,
 149         internalCoeffs_,
 150         psi_.boundaryField().interfaces(),
 151         solverControls
 152     )->solve(psi_.internalField(), totalSource);
 153
 154     solverPerf.print();
 155
 156     diag() = saveDiag;
 157
 158     psi_.correctBoundaryConditions();
 159
 160     return solverPerf;
 161 }
 162
 163
 164 template<>
 165 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
 166 {
 167     scalarField boundaryDiag(psi_.size(), 0.0);
 168     addBoundaryDiag(boundaryDiag, 0);
 169
 170     tmp<scalarField> tres
 171     (
 172         lduMatrix::residual
 173         (
 174             psi_.internalField(),
 175             source_ - boundaryDiag*psi_.internalField(),
 176             boundaryCoeffs_,
 177             psi_.boundaryField().interfaces(),
 178             0
 179         )
 180     );
 181
 182     addBoundarySource(tres());
 183
 184     return tres;
 185 }
 186
 187
 188 template<>
 189 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H() const
 190 {
 191     tmp<volScalarField> tHphi
 192     (
 193         new volScalarField
 194         (
 195             IOobject
 196             (
 197                 "H("+psi_.name()+')',
 198                 psi_.instance(),
 199                 psi_.mesh(),
 200                 IOobject::NO_READ,
 201                 IOobject::NO_WRITE
 202             ),
 203             psi_.mesh(),
 204             dimensions_/dimVol,
 205             zeroGradientFvPatchScalarField::typeName
 206         )
 207     );
 208     volScalarField& Hphi = tHphi();
 209
 210     Hphi.internalField() = (lduMatrix::H(psi_.internalField()) + source_);
 211     addBoundarySource(Hphi.internalField());
 212
 213     Hphi.internalField() /= psi_.mesh().V();
 214     Hphi.correctBoundaryConditions();
 215
 216     return tHphi;
 217 }
 218
 219
 220 template<>
 221 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H1() const
 222 {
 223     tmp<volScalarField> tH1
 224     (
 225         new volScalarField
 226         (
 227             IOobject
 228             (
 229                 "H(1)",
 230                 psi_.instance(),
 231                 psi_.mesh(),
 232                 IOobject::NO_READ,
 233                 IOobject::NO_WRITE
 234             ),
 235             psi_.mesh(),
 236             dimensions_/(dimVol*psi_.dimensions()),
 237             zeroGradientFvPatchScalarField::typeName
 238         )
 239     );
 240     volScalarField& H1_ = tH1();
 241
 242     H1_.internalField() = lduMatrix::H1();
 243     //addBoundarySource(Hphi.internalField());
 244
 245     H1_.internalField() /= psi_.mesh().V();
 246     H1_.correctBoundaryConditions();
 247
 248     return tH1;
 249 }
 250
 251
 252 // ************************************************************************* //