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