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         lduSolverPerformance solverPerf;
  99
 100         // Solver call
 101         solverPerf = lduSolver::New
 102         (
 103             psi_.name() + pTraits<Type>::componentNames[cmpt],
 104             *this,
 105             bouCoeffsCmpt,
 106             intCoeffsCmpt,
 107             interfaces,
 108             solverControls
 109         )->solve(psiCmpt, sourceCmpt, cmpt);
 110
 111         solverPerf.print();
 112
 113         if
 114         (
 115             solverPerf.initialResidual() > solverPerfVec.initialResidual()
 116          && !solverPerf.singular()
 117         )
 118         {
 119             solverPerfVec = solverPerf;
 120         }
 121
 122         psi_.internalField().replace(cmpt, psiCmpt);
 123         diag() = saveDiag;
 124     }
 125
 126     psi_.correctBoundaryConditions();
 127
 128     return solverPerfVec;
 129 }
 130
 131
 132 template<class Type>
 133 autoPtr<typename faMatrix<Type>::faSolver> faMatrix<Type>::solver()
 134 {
 135     return solver(psi_.mesh().solverDict(psi_.name()));
 136 }
 137
 138 template<class Type>
 139 lduSolverPerformance faMatrix<Type>::faSolver::solve()
 140 {
 141     return solve(psi_.mesh().solverDict(psi_.name()));
 142 }
 143
 144
 145 template<class Type>
 146 lduSolverPerformance faMatrix<Type>::solve()
 147 {
 148     return solve(psi_.mesh().solverDict(psi_.name()));
 149 }
 150
 151
 152 // Return the matrix residual
 153 template<class Type>
 154 tmp<Field<Type> > faMatrix<Type>::residual() const
 155 {
 156     tmp<Field<Type> > tres(source_);
 157     Field<Type>& res = tres();
 158
 159     addBoundarySource(res);
 160
 161     // Make a copy of interfaces: no longer a reference
 162     // HJ, 20/Nov/2007
 163     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 164
 165     // Loop over field components
 166     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
 167     {
 168         scalarField psiCmpt = psi_.internalField().component(cmpt);
 169
 170         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 171         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 172
 173         FieldField<Field, scalar> bouCoeffsCmpt
 174         (
 175             boundaryCoeffs_.component(cmpt)
 176         );
 177
 178         res.replace
 179         (
 180             cmpt,
 181             lduMatrix::residual
 182             (
 183                 psiCmpt,
 184                 res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
 185                 bouCoeffsCmpt,
 186                 interfaces,
 187                 cmpt
 188             )
 189         );
 190     }
 191
 192     return tres;
 193 }
 194
 195
 196 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 197
 198 } // End namespace Foam
 199
 200 // ************************************************************************* //