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