src/finiteArea/faMatrices/faMatrix/faMatrixSolve.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 Description
  25     Finite-Area matrix basic solvers.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 namespace Foam
  30 {
  31
  32 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  33
  34 // Set reference level for a component of the solution
  35 // on a given patch face
  36 template<class Type>
  37 void faMatrix<Type>::setComponentReference
  38 (
  39     const label patchi,
  40     const label facei,
  41     const direction cmpt,
  42     const scalar value
  43 )
  44 {
  45     internalCoeffs_[patchi][facei].component(cmpt) +=
  46         diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
  47
  48     boundaryCoeffs_[patchi][facei].component(cmpt) +=
  49         diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]*value;
  50 }
  51
  52
  53 template<class Type>
  54 lduSolverPerformance faMatrix<Type>::solve(const dictionary& solverControls)
  55 {
  56     if (debug)
  57     {
  58         Info<< "faMatrix<Type>::solve(const dictionary&) : "
  59                "solving faMatrix<Type>"
  60             << endl;
  61     }
  62
  63     lduSolverPerformance solverPerfVec
  64     (
  65         "faMatrix<Type>::solve",
  66         psi_.name()
  67     );
  68
  69     scalarField saveDiag = diag();
  70
  71     Field<Type> source = source_;
  72     addBoundarySource(source);
  73
  74     // Make a copy of interfaces: no longer a reference
  75     // HJ, 20/Nov/2007
  76     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  77
  78     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
  79     {
  80         // copy field and source
  81
  82         scalarField psiCmpt = psi_.internalField().component(cmpt);
  83         addBoundaryDiag(diag(), solvingComponent);
  84
  85         scalarField sourceCmpt = source.component(cmpt);
  86
  87         FieldField<Field, scalar> bouCoeffsCmpt
  88         (
  89             boundaryCoeffs_.component(cmpt)
  90         );
  91
  92         FieldField<Field, scalar> intCoeffsCmpt
  93         (
  94             internalCoeffs_.component(cmpt)
  95         );
  96
  97         // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
  98         // bouCoeffsCmpt for the explicit part of the coupled boundary
  99         // conditions
 100         initMatrixInterfaces
 101         (
 102             bouCoeffsCmpt,
 103             interfaces,
 104             psiCmpt,
 105             sourceCmpt,
 106             cmpt
 107         );
 108
 109         updateMatrixInterfaces
 110         (
 111             bouCoeffsCmpt,
 112             interfaces,
 113             psiCmpt,
 114             sourceCmpt,
 115             cmpt
 116         );
 117
 118         lduSolverPerformance solverPerf;
 119
 120         // Solver call
 121         solverPerf = lduSolver::New
 122         (
 123             psi_.name() + pTraits<Type>::componentNames[cmpt],
 124             *this,
 125             bouCoeffsCmpt,
 126             intCoeffsCmpt,
 127             interfaces,
 128             solverControls
 129         )->solve(psiCmpt, sourceCmpt, cmpt);
 130
 131         solverPerf.print();
 132
 133         if
 134         (
 135             solverPerf.initialResidual() > solverPerfVec.initialResidual()
 136          && !solverPerf.singular()
 137         )
 138         {
 139             solverPerfVec = solverPerf;
 140         }
 141
 142         psi_.internalField().replace(cmpt, psiCmpt);
 143         diag() = saveDiag;
 144     }
 145
 146     psi_.correctBoundaryConditions();
 147
 148     return solverPerfVec;
 149 }
 150
 151
 152 template<class Type>
 153 autoPtr<typename faMatrix<Type>::faSolver> faMatrix<Type>::solver()
 154 {
 155     return solver(psi_.mesh().solutionDict().solverDict(psi_.name()));
 156 }
 157
 158 template<class Type>
 159 lduSolverPerformance faMatrix<Type>::faSolver::solve()
 160 {
 161     return solvei
 162     (
 163         faMat_.psi().mesh().solutionDict().solverDict
 164         (
 165             faMat_.psi().name()
 166         )
 167     );
 168 }
 169
 170
 171 template<class Type>
 172 lduSolverPerformance faMatrix<Type>::solve()
 173 {
 174     return solve
 175     (
 176         this->psi().mesh().solutionDict().solverDict
 177         (
 178             this->psi().name()
 179         )
 180     );
 181 }
 182
 183
 184 // Return the matrix residual
 185 template<class Type>
 186 tmp<Field<Type> > faMatrix<Type>::residual() const
 187 {
 188     tmp<Field<Type> > tres(source_);
 189     Field<Type>& res = tres();
 190
 191     addBoundarySource(res);
 192
 193     // Make a copy of interfaces: no longer a reference
 194     // HJ, 20/Nov/2007
 195     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 196
 197     // Loop over field components
 198     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
 199     {
 200         scalarField psiCmpt = psi_.internalField().component(cmpt);
 201
 202         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 203         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 204
 205         FieldField<Field, scalar> bouCoeffsCmpt
 206         (
 207             boundaryCoeffs_.component(cmpt)
 208         );
 209
 210         res.replace
 211         (
 212             cmpt,
 213             lduMatrix::residual
 214             (
 215                 psiCmpt,
 216                 res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
 217                 bouCoeffsCmpt,
 218                 interfaces,
 219                 cmpt
 220             )
 221         );
 222     }
 223
 224     return tres;
 225 }
 226
 227
 228 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 229
 230 } // End namespace Foam
 231
 232 // ************************************************************************* //