GAMGSolver.H 11.9 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8 9
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2019 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29

Class
    Foam::GAMGSolver

30 31 32
Group
    grpLduMatrixSolvers

33 34 35
Description
    Geometric agglomerated algebraic multigrid solver.

36
    Characteristics:
37 38 39 40 41 42 43 44 45 46 47
      - Requires positive definite, diagonally dominant matrix.
      - Agglomeration algorithm: selectable and optionally cached.
      - Restriction operator: summation.
      - Prolongation operator: injection.
      - Smoother: Gauss-Seidel.
      - Coarse matrix creation: central coefficient: summation of fine grid
        central coefficients with the removal of intra-cluster face;
        off-diagonal coefficient: summation of off-diagonal faces.
      - Coarse matrix scaling: performed by correction scaling, using steepest
        descent optimisation.
      - Type of cycle: V-cycle with optional pre-smoothing.
48
      - Coarsest-level matrix solved using PCG or PBiCGStab.
49 50 51

SourceFiles
    GAMGSolver.C
52
    GAMGSolverAgglomerateMatrix.C
53 54
    GAMGSolverInterpolate.C
    GAMGSolverScale.C
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
    GAMGSolverSolve.C

\*---------------------------------------------------------------------------*/

#ifndef GAMGSolver_H
#define GAMGSolver_H

#include "GAMGAgglomeration.H"
#include "lduMatrix.H"
#include "labelField.H"
#include "primitiveFields.H"
#include "LUscalarMatrix.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
74
                         Class GAMGSolver Declaration
75 76 77 78 79 80 81 82
\*---------------------------------------------------------------------------*/

class GAMGSolver
:
    public lduMatrix::solver
{
    // Private data

83
        bool cacheAgglomeration_;
84 85 86 87

        //- Number of pre-smoothing sweeps
        label nPreSweeps_;

88 89 90 91 92 93
        //- Lever multiplier for the number of pre-smoothing sweeps
        label preSweepsLevelMultiplier_;

        //- Maximum number of pre-smoothing sweeps
        label maxPreSweeps_;

94 95 96
        //- Number of post-smoothing sweeps
        label nPostSweeps_;

97 98 99 100 101 102
        //- Lever multiplier for the number of post-smoothing sweeps
        label postSweepsLevelMultiplier_;

        //- Maximum number of post-smoothing sweeps
        label maxPostSweeps_;

103 104 105
        //- Number of smoothing sweeps on finest mesh
        label nFinestSweeps_;

106 107 108 109
        //- Choose if the corrections should be interpolated after injection.
        //  By default corrections are not interpolated.
        bool interpolateCorrection_;

110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
        //- Choose if the corrections should be scaled.
        //  By default corrections for symmetric matrices are scaled
        //  but not for asymmetric matrices.
        bool scaleCorrection_;

        //- Direct or iteratively solve the coarsest level
        bool directSolveCoarsest_;

        //- The agglomeration
        const GAMGAgglomeration& agglomeration_;

        //- Hierarchy of matrix levels
        PtrList<lduMatrix> matrixLevels_;

        //- Hierarchy of interfaces.
125
        PtrList<PtrList<lduInterfaceField>> primitiveInterfaceLevels_;
126

127 128 129
        //- Hierarchy of interfaces in lduInterfaceFieldPtrs form
        PtrList<lduInterfaceFieldPtrsList> interfaceLevels_;

130
        //- Hierarchy of interface boundary coefficients
131
        PtrList<FieldField<Field, scalar>> interfaceLevelsBouCoeffs_;
132 133

        //- Hierarchy of interface internal coefficients
134
        PtrList<FieldField<Field, scalar>> interfaceLevelsIntCoeffs_;
135

Andrew Heather's avatar
Andrew Heather committed
136
        //- LU decomposed coarsest matrix
137 138
        autoPtr<LUscalarMatrix> coarsestLUMatrixPtr_;

139 140 141
        //- Sparse coarsest matrix solver
        autoPtr<lduMatrix::solver> coarsestSolverPtr_;

142

143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
    // Private Member Functions

        //- Read control parameters from the control dictionary
        virtual void readControls();

        //- Simplified access to interface level
        const lduInterfaceFieldPtrsList& interfaceLevel
        (
            const label i
        ) const;

        //- Simplified access to matrix level
        const lduMatrix& matrixLevel(const label i) const;

        //- Simplified access to interface boundary coeffs level
        const FieldField<Field, scalar>& interfaceBouCoeffsLevel
        (
            const label i
        ) const;

        //- Simplified access to interface internal coeffs level
        const FieldField<Field, scalar>& interfaceIntCoeffsLevel
        (
            const label i
        ) const;

169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
        //- Agglomerate coarse matrix. Supply mesh to use - so we can
        //  construct temporary matrix on the fine mesh (instead of the coarse
        //  mesh)
        void agglomerateMatrix
        (
            const label fineLevelIndex,
            const lduMesh& coarseMesh,
            const lduInterfacePtrsList& coarseMeshInterfaces
        );

        //- Agglomerate coarse interface coefficients
        void agglomerateInterfaceCoefficients
        (
            const label fineLevelIndex,
            const lduInterfacePtrsList& coarseMeshInterfaces,
184
            PtrList<lduInterfaceField>& coarsePrimInterfaces,
185 186 187 188
            lduInterfaceFieldPtrsList& coarseInterfaces,
            FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
            FieldField<Field, scalar>& coarseInterfaceIntCoeffs
        ) const;
189

190 191 192 193
        //- Collect matrices from other processors
        void gatherMatrices
        (
            const labelList& procIDs,
194 195 196 197 198 199 200 201 202
            const lduMesh& dummyMesh,
            const label meshComm,

            const lduMatrix& mat,
            const FieldField<Field, scalar>& interfaceBouCoeffs,
            const FieldField<Field, scalar>& interfaceIntCoeffs,
            const lduInterfaceFieldPtrsList& interfaces,

            PtrList<lduMatrix>& otherMats,
203 204
            PtrList<FieldField<Field, scalar>>& otherBouCoeffs,
            PtrList<FieldField<Field, scalar>>& otherIntCoeffs,
205
            List<boolList>& otherTransforms,
206
            List<List<label>>& otherRanks
207 208
        ) const;

209 210 211 212 213
        //- Agglomerate processor matrices
        void procAgglomerateMatrix
        (
            // Agglomeration information
            const labelList& procAgglomMap,
214
            const List<label>& agglomProcIDs,
215

216 217
            const label levelI,

218 219 220 221 222 223
            // Resulting matrix
            autoPtr<lduMatrix>& allMatrixPtr,
            FieldField<Field, scalar>& allInterfaceBouCoeffs,
            FieldField<Field, scalar>& allInterfaceIntCoeffs,
            PtrList<lduInterfaceField>& allPrimitiveInterfaces,
            lduInterfaceFieldPtrsList& allInterfaces
224 225 226 227 228 229
        ) const;

        //- Agglomerate processor matrices
        void procAgglomerateMatrix
        (
            const labelList& procAgglomMap,
230
            const List<label>& agglomProcIDs,
231
            const label levelI
232 233
        );

234 235 236
        //- Interpolate the correction after injected prolongation
        void interpolate
        (
237 238
            solveScalarField& psi,
            solveScalarField& Apsi,
239 240 241 242 243 244 245
            const lduMatrix& m,
            const FieldField<Field, scalar>& interfaceBouCoeffs,
            const lduInterfaceFieldPtrsList& interfaces,
            const direction cmpt
        ) const;

        //- Interpolate the correction after injected prolongation and
246
        //  re-normalise
247 248
        void interpolate
        (
249 250
            solveScalarField& psi,
            solveScalarField& Apsi,
251 252 253
            const lduMatrix& m,
            const FieldField<Field, scalar>& interfaceBouCoeffs,
            const lduInterfaceFieldPtrsList& interfaces,
254
            const labelList& restrictAddressing,
255
            const solveScalarField& psiC,
256 257 258 259
            const direction cmpt
        ) const;

        //- Calculate and apply the scaling factor from Acf, coarseSource
260 261
        //  and coarseField.
        //  At the same time do a Jacobi iteration on the coarseField using
262
        //  the Acf provided after the coarseField values are used for the
263
        //  scaling factor.
264
        void scale
265
        (
266 267
            solveScalarField& field,
            solveScalarField& Acf,
268 269 270
            const lduMatrix& A,
            const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
            const lduInterfaceFieldPtrsList& interfaceLevel,
271
            const solveScalarField& source,
272 273 274 275 276 277
            const direction cmpt
        ) const;

        //- Initialise the data structures for the V-cycle
        void initVcycle
        (
278 279
            PtrList<solveScalarField>& coarseCorrFields,
            PtrList<solveScalarField>& coarseSources,
280
            PtrList<lduMatrix::smoother>& smoothers,
281 282
            solveScalarField& scratch1,
            solveScalarField& scratch2
283 284 285 286 287 288 289
        ) const;


        //- Perform a single GAMG V-cycle with pre, post and finest smoothing.
        void Vcycle
        (
            const PtrList<lduMatrix::smoother>& smoothers,
290
            solveScalarField& psi,
291
            const scalarField& source,
292 293 294
            solveScalarField& Apsi,
            solveScalarField& finestCorrection,
            solveScalarField& finestResidual,
295

296 297
            solveScalarField& scratch1,
            solveScalarField& scratch2,
298

299 300
            PtrList<solveScalarField>& coarseCorrFields,
            PtrList<solveScalarField>& coarseSources,
301 302 303
            const direction cmpt=0
        ) const;

304 305 306 307 308 309 310 311
        //- Create and return the dictionary to specify the PCG solver
        //  to solve the coarsest level
        dictionary PCGsolverDict
        (
            const scalar tol,
            const scalar relTol
        ) const;

312
        //- Create and return the dictionary to specify the PBiCGStab solver
313
        //  to solve the coarsest level
314
        dictionary PBiCGStabSolverDict
315 316 317 318
        (
            const scalar tol,
            const scalar relTol
        ) const;
319 320 321 322

        //- Solve the coarsest level with either an iterative or direct solver
        void solveCoarsestLevel
        (
323 324
            solveScalarField& coarsestCorrField,
            const solveScalarField& coarsestSource
325 326 327 328 329 330 331 332 333 334 335 336 337
        ) const;


public:

    friend class GAMGPreconditioner;

    //- Runtime type information
    TypeName("GAMG");


    // Constructors

338
        //- Construct from lduMatrix and solver controls
339 340 341 342 343 344 345
        GAMGSolver
        (
            const word& fieldName,
            const lduMatrix& matrix,
            const FieldField<Field, scalar>& interfaceBouCoeffs,
            const FieldField<Field, scalar>& interfaceIntCoeffs,
            const lduInterfaceFieldPtrsList& interfaces,
346
            const dictionary& solverControls
347 348 349
        );


350 351
    //- Destructor
    virtual ~GAMGSolver();
352 353 354 355 356


    // Member Functions

        //- Solve
357
        virtual solverPerformance solve
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
        (
            scalarField& psi,
            const scalarField& source,
            const direction cmpt=0
        ) const;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //