src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
   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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Description
  26     Segregated solver for symmetric and asymmetric matrices.
  27     Calls scalar solver component-by-component
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "SegregatedSolver.H"
  32
  33 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  34
  35 // Construct from matrix and solver data stream
  36 template<class Type>
  37 Foam::SegregatedSolver<Type>::SegregatedSolver
  38 (
  39     const word& fieldName,
  40     const BlockLduMatrix<Type>& matrix,
  41     const dictionary& dict
  42 )
  43 :
  44     BlockLduSolver<Type>(fieldName, matrix, dict),
  45     scalarX_(matrix.lduAddr().size()),
  46     scalarMatrix_(matrix.mesh()),
  47     scalarB_(matrix.lduAddr().size())
  48 {}
  49
  50
  51 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  52
  53 template<class Type>
  54 Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
  55 (
  56     Field<Type>& x,
  57     const Field<Type>& b
  58 )
  59 {
  60     // Get reference to matrix, x and b
  61     const BlockLduMatrix<Type>& blockMatrix = this->matrix_;
  62
  63     // Determine if the diagonal or off-diagonal is expanded
  64
  65     // Check switching diagonal
  66     bool switchingDiag = false;
  67
  68     if (blockMatrix.diag().activeType() == blockCoeffBase::SCALAR)
  69     {
  70         scalarMatrix_.diag() = blockMatrix.diag().asScalar();
  71     }
  72     else
  73     {
  74         switchingDiag = true;
  75     }
  76
  77     // Check switching upper
  78     bool switchingUpper = false;
  79
  80     if (blockMatrix.thereIsUpper())
  81     {
  82         if (blockMatrix.upper().activeType() == blockCoeffBase::SCALAR)
  83         {
  84             scalarMatrix_.upper() = blockMatrix.upper().asScalar();
  85         }
  86         else
  87         {
  88             switchingUpper = true;
  89         }
  90     }
  91
  92     // Check switching lower
  93     bool switchingLower = false;
  94
  95     if (blockMatrix.thereIsLower())
  96     {
  97         if (blockMatrix.lower().activeType() == blockCoeffBase::SCALAR)
  98         {
  99             scalarMatrix_.lower() = blockMatrix.lower().asScalar();
 100         }
 101         else
 102         {
 103             switchingLower = true;
 104         }
 105     }
 106
 107     // Decouple quadratic coupling by multiplying out the square coefficient
 108     // coupling
 109     Field<Type>* dBPtr = NULL;
 110
 111     if (blockMatrix.componentCoupled())
 112     {
 113         if (BlockLduMatrix<Type>::debug >= 2)
 114         {
 115             Info << " Component coupled segregation" << endl;
 116         }
 117
 118         dBPtr = new Field<Type>(b);
 119         blockMatrix.segregateB(*dBPtr, x);
 120     }
 121
 122     // Prepare solver performance
 123     word segSolverName(this->dict().lookup("solver"));
 124
 125     BlockSolverPerformance<Type> solverPerf
 126     (
 127         typeName + "_" + segSolverName,
 128         this->fieldName()
 129     );
 130
 131     // Loop through the form component by component and call the
 132     // scalar solver
 133     for
 134     (
 135         direction cmpt = 0;
 136         cmpt < pTraits<Type>::nComponents;
 137         cmpt++
 138     )
 139     {
 140         // Grab the x and b for the current component
 141         scalarX_ = x.component(cmpt);
 142
 143         // Handle b decoupling
 144         if (dBPtr)
 145         {
 146             scalarB_ = dBPtr->component(cmpt);
 147         }
 148         else
 149         {
 150             scalarB_ = b.component(cmpt);
 151         }
 152
 153         // Switch diagonal, upper and lower
 154         if (switchingDiag)
 155         {
 156             scalarMatrix_.diag() = blockMatrix.diag().component(cmpt);
 157         }
 158
 159         if (switchingUpper)
 160         {
 161             scalarMatrix_.upper() = blockMatrix.upper().component(cmpt);
 162         }
 163
 164         if (switchingLower)
 165         {
 166             scalarMatrix_.lower() = blockMatrix.lower().component(cmpt);
 167         }
 168
 169         // Call the scalar solver
 170         BlockSolverPerformance<scalar> scalarPerf =
 171             blockScalarSolver::New
 172             (
 173                 this->fieldName(),
 174                 scalarMatrix_,
 175                 this->dict()
 176             )->solve(scalarX_, scalarB_);
 177
 178         // Replace the solution component
 179         x.replace(cmpt, scalarX_);
 180
 181         // Grab solver performance and combine
 182         solverPerf.initialResidual().replace
 183         (
 184             cmpt,
 185             scalarPerf.initialResidual()
 186         );
 187
 188         solverPerf.finalResidual().replace
 189         (
 190             cmpt,
 191             scalarPerf.finalResidual()
 192         );
 193
 194         solverPerf.nIterations() =
 195             max
 196             (
 197                 solverPerf.nIterations(),
 198                 scalarPerf.nIterations()
 199             );
 200
 201         solverPerf.converged() =
 202             (solverPerf.converged() && scalarPerf.converged());
 203
 204         solverPerf.singular() =
 205             solverPerf.singular() && scalarPerf.singular();
 206     }
 207
 208     // Clear decoupled b if allocated
 209     if (dBPtr)
 210     {
 211         delete dBPtr;
 212         dBPtr = NULL;
 213     }
 214
 215     return solverPerf;
 216 }
 217
 218
 219 // ************************************************************************* //