src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "fvScalarMatrix.H"
  27 #include "zeroGradientFvPatchFields.H"
  28
  29 #include "profiling.H"
  30
  31 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  32
  33 template<>
  34 void Foam::fvMatrix<Foam::scalar>::setComponentReference
  35 (
  36     const label patchi,
  37     const label facei,
  38     const direction,
  39     const scalar value
  40 )
  41 {
  42     if (psi_.needReference())
  43     {
  44         if (Pstream::master())
  45         {
  46             internalCoeffs_[patchi][facei] +=
  47                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
  48
  49             boundaryCoeffs_[patchi][facei] +=
  50                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
  51                *value;
  52         }
  53     }
  54 }
  55
  56
  57 template<>
  58 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
  59 (
  60     const dictionary& solverControls
  61 )
  62 {
  63     profilingTrigger profSolve("fvMatrix::solve_"+psi_.name());
  64
  65     if (debug)
  66     {
  67         Info<< "fvMatrix<scalar>::solve(const dictionary& solverControls) : "
  68                "solving fvMatrix<scalar>"
  69             << endl;
  70     }
  71
  72     // Complete matrix assembly.  HJ, 17/Apr/2012
  73     completeAssembly();
  74
  75     scalarField saveDiag = diag();
  76     addBoundaryDiag(diag(), 0);
  77
  78     scalarField totalSource = source_;
  79     addBoundarySource(totalSource, false);
  80
  81     // Make a copy of interfaces: no longer a reference
  82     // HJ, 20/Nov/2007
  83     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  84
  85     // Solver call
  86     lduSolverPerformance solverPerf = lduSolver::New
  87     (
  88         psi_.name(),
  89         *this,
  90         boundaryCoeffs_,
  91         internalCoeffs_,
  92         interfaces,
  93         solverControls
  94     )->solve(psi_.internalField(), totalSource);
  95
  96     solverPerf.print();
  97
  98     // Diagonal has been restored, clear complete assembly flag?
  99     diag() = saveDiag;
 100
 101     psi_.correctBoundaryConditions();
 102
 103     return solverPerf;
 104 }
 105
 106
 107 template<>
 108 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
 109 {
 110     // Complete matrix assembly.  HJ, 17/Apr/2012
 111     fvMatrix<scalar>& m = const_cast<fvMatrix<scalar>&>(*this);
 112     m.completeAssembly();
 113
 114     // To avoid copying the diagonal, execute boundary diag multiplication
 115     // separately.  See source provided to lduMatrix::residual
 116     // HJ, 26/Jun/2015
 117     scalarField boundaryDiag(psi_.size(), 0.0);
 118     addBoundaryDiag(boundaryDiag, 0);
 119
 120     // Bug fix: boundary source must be added on rhs
 121     // HJ, 26/Jun/2015
 122     scalarField totalSource = source_;
 123     addBoundarySource(totalSource, false);
 124
 125    // Make a copy of interfaces: no longer a reference
 126     // HJ, 20/Nov/2007
 127     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 128
 129     tmp<scalarField> tres
 130     (
 131         lduMatrix::residual
 132         (
 133             psi_.internalField(),
 134             totalSource - boundaryDiag*psi_.internalField(),
 135             boundaryCoeffs_,
 136             interfaces,
 137             0
 138         )
 139     );
 140
 141     return tres;
 142 }
 143
 144
 145 template<>
 146 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H() const
 147 {
 148     tmp<volScalarField> tHphi
 149     (
 150         new volScalarField
 151         (
 152             IOobject
 153             (
 154                 "H("+psi_.name()+')',
 155                 psi_.instance(),
 156                 psi_.mesh(),
 157                 IOobject::NO_READ,
 158                 IOobject::NO_WRITE
 159             ),
 160             psi_.mesh(),
 161             dimensions_/dimVol,
 162             zeroGradientFvPatchScalarField::typeName
 163         )
 164     );
 165     volScalarField& Hphi = tHphi();
 166
 167     Hphi.internalField() = (lduMatrix::H(psi_.internalField()) + source_);
 168     addBoundarySource(Hphi.internalField());
 169
 170     Hphi.internalField() /= psi_.mesh().V();
 171     Hphi.correctBoundaryConditions();
 172
 173     return tHphi;
 174 }
 175
 176
 177 template<>
 178 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Foam::scalar>::H1() const
 179 {
 180     tmp<volScalarField> tH1
 181     (
 182         new volScalarField
 183         (
 184             IOobject
 185             (
 186                 "H(1)",
 187                 psi_.instance(),
 188                 psi_.mesh(),
 189                 IOobject::NO_READ,
 190                 IOobject::NO_WRITE
 191             ),
 192             psi_.mesh(),
 193             dimensions_/(dimVol*psi_.dimensions()),
 194             zeroGradientFvPatchScalarField::typeName
 195         )
 196     );
 197     volScalarField& H1_ = tH1();
 198
 199     H1_.internalField() = lduMatrix::H1();
 200     //addBoundarySource(Hphi.internalField());
 201
 202     H1_.internalField() /= psi_.mesh().V();
 203     H1_.correctBoundaryConditions();
 204
 205     return tH1;
 206 }
 207
 208
 209 // ************************************************************************* //