src/coupledMatrix/coupledFvMatrices/coupledFvMatrix/coupledFvMatrix.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     Coupled Finite Volume matrix
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "coupledLduSolver.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 template<class Type>
  34 void Foam::coupledFvMatrix<Type>::checkSize() const
  35 {
  36     if (this->size() < 1)
  37     {
  38         FatalErrorIn("void Foam::coupledFvMatrix<Type>::checkSize() const")
  39             << "No matrices added to coupled matrix"
  40             << abort(FatalError);
  41     }
  42 }
  43
  44
  45 template<class Type>
  46 Foam::word Foam::coupledFvMatrix<Type>::coupledPsiName() const
  47 {
  48     word cpn;
  49
  50     const PtrList<lduMatrix>& matrices = *this;
  51
  52     forAll (matrices, rowI)
  53     {
  54         const fvMatrix<Type>& curMatrix =
  55             static_cast<const fvMatrix<Type>& >(matrices[rowI]);
  56
  57         cpn += curMatrix.psi().name();
  58
  59         if (rowI < matrices.size() - 1)
  60         {
  61             cpn += "+";
  62         }
  63     }
  64
  65     return cpn;
  66 }
  67
  68 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  69
  70 template<class Type>
  71 Foam::coupledSolverPerformance
  72 Foam::coupledFvMatrix<Type>::solve(const dictionary& solverControls)
  73 {
  74     if (debug)
  75     {
  76         InfoIn("coupledFvMatrix<Type>::solve(const dictionary&)")
  77             << "solving coupledFvMatrix<Type>" << endl;
  78     }
  79
  80     coupledSolverPerformance solverPerfVec
  81     (
  82         "coupledFvMatrix<Type>::solve",
  83         this->coupledPsiName()
  84     );
  85
  86     typedef FieldField<Field, scalar> scalarFieldField;
  87
  88     PtrList<lduMatrix>& matrices = *this;
  89
  90     // Make a copy of the diagonal and complete the source
  91     scalarFieldField saveDiag(this->size());
  92     scalarFieldField psiCmpt(this->size());
  93     FieldField<Field, Type> source(this->size());
  94     scalarFieldField sourceCmpt(this->size());
  95     lduInterfaceFieldPtrsListList interfaces(this->size());
  96
  97     // Prepare block solution
  98     forAll (matrices, rowI)
  99     {
 100         fvMatrix<Type>& curMatrix =
 101             static_cast<fvMatrix<Type>& >(matrices[rowI]);
 102
 103         saveDiag.set(rowI, new scalarField(curMatrix.diag()));
 104         psiCmpt.set(rowI, new scalarField(curMatrix.psi().size()));
 105         source.set(rowI, new Field<Type>(curMatrix.source()));
 106         sourceCmpt.set(rowI, new scalarField(curMatrix.psi().size()));
 107
 108         curMatrix.addBoundarySource(source[rowI]);
 109
 110         interfaces[rowI] = curMatrix.psi().boundaryField().interfaces();
 111     }
 112
 113     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
 114     {
 115         // Copy field and source
 116
 117         PtrList<FieldField<Field, scalar> > bouCoeffsCmpt(this->size());
 118         PtrList<FieldField<Field, scalar> > intCoeffsCmpt(this->size());
 119
 120         forAll (matrices, rowI)
 121         {
 122
 123             fvMatrix<Type>& curMatrix =
 124                 static_cast<fvMatrix<Type>& >(matrices[rowI]);
 125
 126             psiCmpt[rowI] = curMatrix.psi().internalField().component(cmpt);
 127             curMatrix.addBoundaryDiag(curMatrix.diag(), cmpt);
 128
 129             sourceCmpt[rowI] = source[rowI].component(cmpt);
 130
 131             bouCoeffsCmpt.set
 132             (
 133                 rowI,
 134                 new FieldField<Field, scalar>
 135                 (
 136                     curMatrix.boundaryCoeffs().component(cmpt)
 137                 )
 138             );
 139
 140             intCoeffsCmpt.set
 141             (
 142                 rowI,
 143                 new FieldField<Field, scalar>
 144                 (
 145                     curMatrix.internalCoeffs().component(cmpt)
 146                 )
 147             );
 148         }
 149
 150         this->initMatrixInterfaces
 151         (
 152             bouCoeffsCmpt,
 153             interfaces,
 154             psiCmpt,
 155             sourceCmpt,
 156             cmpt
 157         );
 158
 159         this->updateMatrixInterfaces
 160         (
 161             bouCoeffsCmpt,
 162             interfaces,
 163             psiCmpt,
 164             sourceCmpt,
 165             cmpt
 166         );
 167
 168         coupledSolverPerformance solverPerf
 169         (
 170             "coupledFvMatrix<Type>::solve",
 171             this->coupledPsiName()
 172         );
 173
 174         // Solver call
 175         solverPerf = coupledLduSolver::New
 176         (
 177             this->coupledPsiName() + pTraits<Type>::componentNames[cmpt],
 178             *this,
 179             bouCoeffsCmpt,
 180             intCoeffsCmpt,
 181             interfaces,
 182             solverControls
 183         )->solve(psiCmpt, sourceCmpt, cmpt);
 184
 185         solverPerf.print();
 186
 187         if
 188         (
 189             solverPerf.initialResidual() > solverPerfVec.initialResidual()
 190          && !solverPerf.singular()
 191         )
 192         {
 193             solverPerfVec = solverPerf;
 194         }
 195
 196         // Update solution
 197         forAll (matrices, rowI)
 198         {
 199             fvMatrix<Type>& curMatrix =
 200                 static_cast<fvMatrix<Type>& >(matrices[rowI]);
 201             curMatrix.psi().internalField().replace(cmpt, psiCmpt[rowI]);
 202             curMatrix.diag() = saveDiag[rowI];
 203         }
 204     }
 205
 206     // Correct boundary conditions
 207     forAll (matrices, rowI)
 208     {
 209         fvMatrix<Type>& curMatrix =
 210             static_cast<fvMatrix<Type>& >(matrices[rowI]);
 211         curMatrix.psi().correctBoundaryConditions();
 212     }
 213
 214     return solverPerfVec;
 215 }
 216
 217
 218 template<class Type>
 219 Foam::coupledSolverPerformance Foam::coupledFvMatrix<Type>::solve()
 220 {
 221     // Check matrices are present.  Using first matrix for controls
 222     checkSize();
 223
 224     const fvMatrix<Type>& m =
 225         static_cast<const fvMatrix<Type>& >(this->operator[](0));
 226
 227     return solve(m.psi().mesh().solutionDict().solver(coupledPsiName()));
 228 }
 229
 230
 231 // ************************************************************************* //