src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  28
  29 template<class Type>
  30 void Foam::fvMatrix<Type>::setComponentReference
  31 (
  32     const label patchi,
  33     const label facei,
  34     const direction cmpt,
  35     const scalar value
  36 )
  37 {
  38     if (psi_.needReference())
  39     {
  40         if (Pstream::master())
  41         {
  42             internalCoeffs_[patchi][facei].component(cmpt) +=
  43                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
  44
  45             boundaryCoeffs_[patchi][facei].component(cmpt) +=
  46                 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
  47                *value;
  48         }
  49     }
  50 }
  51
  52
  53 template<class Type>
  54 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
  55 (
  56     Istream& solverControls
  57 )
  58 {
  59     if (debug)
  60     {
  61         Info<< "fvMatrix<Type>::solve(Istream& solverControls) : "
  62                "solving fvMatrix<Type>"
  63             << endl;
  64     }
  65
  66     lduMatrix::solverPerformance solverPerfVec
  67     (
  68         "fvMatrix<Type>::solve",
  69         psi_.name()
  70     );
  71
  72     scalarField saveDiag = diag();
  73
  74     Field<Type> source = source_;
  75
  76     // At this point include the boundary source from the coupled boundaries.
  77     // This is corrected for the implict part by updateMatrixInterfaces within
  78     // the component loop.
  79     addBoundarySource(source);
  80
  81     typename Type::labelType validComponents
  82     (
  83         pow
  84         (
  85             psi_.mesh().directions(),
  86             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
  87         )
  88     );
  89
  90     for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
  91     {
  92         if (validComponents[cmpt] == -1) continue;
  93
  94         // copy field and source
  95
  96         scalarField psiCmpt = psi_.internalField().component(cmpt);
  97         addBoundaryDiag(diag(), cmpt);
  98
  99         scalarField sourceCmpt = source.component(cmpt);
 100
 101         FieldField<Field, scalar> bouCoeffsCmpt
 102         (
 103             boundaryCoeffs_.component(cmpt)
 104         );
 105
 106         FieldField<Field, scalar> intCoeffsCmpt
 107         (
 108             internalCoeffs_.component(cmpt)
 109         );
 110
 111         lduInterfaceFieldPtrsList interfaces =
 112             psi_.boundaryField().interfaces();
 113
 114         // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
 115         // bouCoeffsCmpt for the explicit part of the coupled boundary
 116         // conditions
 117         initMatrixInterfaces
 118         (
 119             bouCoeffsCmpt,
 120             interfaces,
 121             psiCmpt,
 122             sourceCmpt,
 123             cmpt
 124         );
 125
 126         updateMatrixInterfaces
 127         (
 128             bouCoeffsCmpt,
 129             interfaces,
 130             psiCmpt,
 131             sourceCmpt,
 132             cmpt
 133         );
 134
 135         lduMatrix::solverPerformance solverPerf;
 136
 137         // Solver call
 138         solverPerf = lduMatrix::solver::New
 139         (
 140             psi_.name() + pTraits<Type>::componentNames[cmpt],
 141             *this,
 142             bouCoeffsCmpt,
 143             intCoeffsCmpt,
 144             interfaces,
 145             solverControls.rewind()
 146         )->solve(psiCmpt, sourceCmpt, cmpt);
 147
 148         solverPerf.print();
 149
 150         if
 151         (
 152             solverPerf.initialResidual() > solverPerfVec.initialResidual()
 153          && !solverPerf.singular()
 154         )
 155         {
 156             solverPerfVec = solverPerf;
 157         }
 158
 159         psi_.internalField().replace(cmpt, psiCmpt);
 160         diag() = saveDiag;
 161     }
 162
 163     psi_.correctBoundaryConditions();
 164
 165     return solverPerfVec;
 166 }
 167
 168
 169 template<class Type>
 170 Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
 171 Foam::fvMatrix<Type>::solver()
 172 {
 173     return solver(psi_.mesh().solver(psi_.name()));
 174 }
 175
 176 template<class Type>
 177 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
 178 {
 179     return solve(psi_.mesh().solver(psi_.name()));
 180 }
 181
 182
 183 template<class Type>
 184 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
 185 {
 186     return solve(psi_.mesh().solver(psi_.name()));
 187 }
 188
 189
 190 template<class Type>
 191 Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
 192 {
 193     tmp<Field<Type> > tres(source_);
 194     Field<Type>& res = tres();
 195
 196     addBoundarySource(res);
 197
 198     // Loop over field components
 199     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 200     {
 201         scalarField psiCmpt = psi_.internalField().component(cmpt);
 202
 203         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 204         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 205
 206         FieldField<Field, scalar> bouCoeffsCmpt
 207         (
 208             boundaryCoeffs_.component(cmpt)
 209         );
 210
 211         res.replace
 212         (
 213             cmpt,
 214             lduMatrix::residual
 215             (
 216                 psiCmpt,
 217                 res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
 218                 bouCoeffsCmpt,
 219                 psi_.boundaryField().interfaces(),
 220                 cmpt
 221             )
 222         );
 223     }
 224
 225     return tres;
 226 }
 227
 228
 229 // ************************************************************************* //