src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.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 \*---------------------------------------------------------------------------*/
  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     // Make a copy of interfaces: no longer a reference
  74     // HJ, 20/Nov/2007
  75     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  76
  77     autoPtr<fvMatrix<scalar>::fvSolver> solverPtr
  78     (
  79         new fvMatrix<scalar>::fvSolver
  80         (
  81             *this,
  82             lduSolver::New
  83             (
  84                 psi_.name(),
  85                 *this,
  86                 boundaryCoeffs_,
  87                 internalCoeffs_,
  88                 interfaces,
  89                 solverControls
  90             )
  91         )
  92     );
  93
  94     diag() = saveDiag;
  95
  96     return solverPtr;
  97 }
  98
  99
 100 template<>
 101 Foam::lduMatrix::solverPerformance
 102 Foam::fvMatrix<Foam::scalar>::fvSolver::solve
 103 (
 104     const dictionary& solverControls
 105 )
 106 {
 107     scalarField saveDiag = fvMat_.diag();
 108     fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
 109
 110     scalarField totalSource = fvMat_.source();
 111     fvMat_.addBoundarySource(totalSource, false);
 112
 113     // assign new solver controls
 114     solver_->read(solverControls);
 115     lduSolverPerformance solverPerf =
 116         solver_->solve(fvMat_.psi().internalField(), totalSource);
 117
 118     solverPerf.print();
 119
 120     fvMat_.diag() = saveDiag;
 121
 122     fvMat_.psi().correctBoundaryConditions();
 123
 124     return solverPerf;
 125 }
 126
 127
 128 template<>
 129 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
 130 (
 131     const dictionary& solverControls
 132 )
 133 {
 134     if (debug)
 135     {
 136         Info<< "fvMatrix<scalar>::solve(const dictionary& solverControls) : "
 137                "solving fvMatrix<scalar>"
 138             << endl;
 139     }
 140
 141     scalarField saveDiag = diag();
 142     addBoundaryDiag(diag(), 0);
 143
 144     scalarField totalSource = source_;
 145     addBoundarySource(totalSource, false);
 146
 147     // Make a copy of interfaces: no longer a reference
 148     // HJ, 20/Nov/2007
 149     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 150
 151     // Solver call
 152     lduSolverPerformance solverPerf = lduSolver::New
 153     (
 154         psi_.name(),
 155         *this,
 156         boundaryCoeffs_,
 157         internalCoeffs_,
 158         interfaces,
 159         solverControls
 160     )->solve(psi_.internalField(), totalSource);
 161
 162     solverPerf.print();
 163
 164     diag() = saveDiag;
 165
 166     psi_.correctBoundaryConditions();
 167
 168     return solverPerf;
 169 }
 170
 171
 172 template<>
 173 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
 174 {
 175     scalarField boundaryDiag(psi_.size(), 0.0);
 176     addBoundaryDiag(boundaryDiag, 0);
 177
 178     // Make a copy of interfaces: no longer a reference
 179     // HJ, 20/Nov/2007
 180     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 181
 182     tmp<scalarField> tres
 183     (
 184         lduMatrix::residual
 185         (
 186             psi_.internalField(),
 187             source_ - boundaryDiag*psi_.internalField(),
 188             boundaryCoeffs_,
 189             interfaces,
 190             0
 191         )
 192     );
 193
 194     addBoundarySource(tres());
 195
 196     return tres;
 197 }
 198
 199
 200 template<>
 201 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H() const
 202 {
 203     tmp<volScalarField> tHphi
 204     (
 205         new volScalarField
 206         (
 207             IOobject
 208             (
 209                 "H("+psi_.name()+')',
 210                 psi_.instance(),
 211                 psi_.mesh(),
 212                 IOobject::NO_READ,
 213                 IOobject::NO_WRITE
 214             ),
 215             psi_.mesh(),
 216             dimensions_/dimVol,
 217             zeroGradientFvPatchScalarField::typeName
 218         )
 219     );
 220     volScalarField& Hphi = tHphi();
 221
 222     Hphi.internalField() = (lduMatrix::H(psi_.internalField()) + source_);
 223     addBoundarySource(Hphi.internalField());
 224
 225     Hphi.internalField() /= psi_.mesh().V();
 226     Hphi.correctBoundaryConditions();
 227
 228     return tHphi;
 229 }
 230
 231
 232 template<>
 233 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H1() const
 234 {
 235     tmp<volScalarField> tH1
 236     (
 237         new volScalarField
 238         (
 239             IOobject
 240             (
 241                 "H(1)",
 242                 psi_.instance(),
 243                 psi_.mesh(),
 244                 IOobject::NO_READ,
 245                 IOobject::NO_WRITE
 246             ),
 247             psi_.mesh(),
 248             dimensions_/(dimVol*psi_.dimensions()),
 249             zeroGradientFvPatchScalarField::typeName
 250         )
 251     );
 252     volScalarField& H1_ = tH1();
 253
 254     H1_.internalField() = lduMatrix::H1();
 255     //addBoundarySource(Hphi.internalField());
 256
 257     H1_.internalField() /= psi_.mesh().V();
 258     H1_.correctBoundaryConditions();
 259
 260     return tH1;
 261 }
 262
 263
 264 // ************************************************************************* //