src/cudaSolvers/cudaBiCGStab/cudaBiCGStab.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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     Preconditioned Conjugate Gradient solver with run-time selectable
  26     preconditioning
  27
  28 Author
  29     Dominik Christ, Wikki Ltd.
  30     Based on Cufflink library by Daniel P. Combest
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "cudaBiCGStab.H"
  35
  36 // Preconditioner and solver are hardwired due to
  37 // code structure of Cusp library. Below are
  38 // external functions with solver/preconditioner combinations
  39
  40 // BiCG with diagonal preconditioning
  41 extern "C" void bicgDiag
  42 (
  43     cuspEquationSystem* ces,
  44     cudaSolverPerformance* solverParam
  45 );
  46
  47 // BiCG with Ainv preconditioning
  48 extern "C" void bicgAinv
  49 (
  50     cuspEquationSystem* ces,
  51     cudaSolverPerformance* solverPerf,
  52     ValueType drop_tolerance,
  53     int nonzero_per_row,
  54     bool lin_dropping,
  55     int lin_param
  56 );
  57
  58
  59 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  60
  61 namespace Foam
  62 {
  63     defineTypeNameAndDebug(cudaBiCGStab, 0);
  64
  65     lduSolver::addasymMatrixConstructorToTable<cudaBiCGStab>
  66         addcudaBiCGStabAsymMatrixConstructorToTable_;
  67 }
  68
  69
  70 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  71
  72 // - Construct from matrix and solver data stream
  73 Foam::cudaBiCGStab::cudaBiCGStab
  74 (
  75     const word& fieldName,
  76     const lduMatrix& matrix,
  77     const FieldField<Field, scalar>& coupleBouCoeffs,
  78     const FieldField<Field, scalar>& coupleIntCoeffs,
  79     const lduInterfaceFieldPtrsList& interfaces,
  80     const dictionary& dict
  81 )
  82 :
  83     cudaSolver
  84     (
  85         fieldName,
  86         matrix,
  87         coupleBouCoeffs,
  88         coupleIntCoeffs,
  89         interfaces,
  90         dict
  91     )
  92 {}
  93
  94
  95 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  96
  97 Foam::lduSolverPerformance Foam::cudaBiCGStab::solve
  98 (
  99     scalarField& x,
 100     const scalarField& b,
 101     const direction cmpt
 102 ) const
 103 {
 104     // Initialize Cusp solver perfomance
 105     cudaSolverPerformance solverPerf = cudaSolverPerformanceDefault();
 106     solverPerf.minIter = minIter(); // Minimum iterations
 107     solverPerf.maxIter = maxIter(); // Maximum iterations
 108     solverPerf.relTol = relTolerance(); // Relative tolerance
 109     solverPerf.tol = tolerance(); // Tolerance
 110
 111     if (lduMatrix::debug >= 2)
 112     {
 113          solverPerf.debugCusp = true;
 114     }
 115
 116     // Initialize and copy matrix data to GPU
 117     cuspEquationSystem ces = createAsymCuspMatrix(matrix(), x, b);
 118
 119     // Call solver externally
 120     word preconName(dict().lookup("preconditioner"));
 121     if (preconName == "diagonal")
 122     {
 123         bicgDiag(&ces, &solverPerf);
 124     }
 125     else if(preconName == "Ainv")
 126     {
 127         bicgAinv
 128         (
 129             &ces,
 130             &solverPerf,
 131             dict().lookupOrDefault<scalar>("dropTolerance", 0.1),
 132             dict().lookupOrDefault<label>("nonzeroPerRow", -1),
 133             dict().lookupOrDefault<bool>("LinDropping", false),
 134             dict().lookupOrDefault<label>("LinParameter", 1)
 135         );
 136     }
 137     else
 138     {
 139         FatalErrorIn("cudaBiCGStab::solver()")
 140             << "Unknown preconditioner name. "
 141             << "Options are:" << nl
 142             << "(" << nl
 143             << "diagonal" << nl
 144             << "Ainv" << nl
 145             << ")" << nl
 146             << abort(FatalError);
 147     }
 148
 149     // copy the x vector back to Openfoam
 150         thrust::copy(ces.X.begin(), ces.X.end(), x.begin());
 151
 152     // Return solver output
 153     return lduSolverPerformance
 154     (
 155         typeName,
 156         fieldName(),
 157         solverPerf.iRes,
 158         solverPerf.fRes,
 159         solverPerf.nIterations,
 160         solverPerf.converged,
 161         solverPerf.singular
 162     );
 163
 164 }
 165
 166
 167 // ************************************************************************* //