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