symGaussSeidelSmoother.C 6.82 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
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8 9
    Copyright (C) 2012-2015 OpenFOAM Foundation
    Copyright (C) 2017-2019 OpenCFD Ltd.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    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.

    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
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "symGaussSeidelSmoother.H"
30
#include "PrecisionAdaptor.H"
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(symGaussSeidelSmoother, 0);

    lduMatrix::smoother::addsymMatrixConstructorToTable<symGaussSeidelSmoother>
        addsymGaussSeidelSmootherSymMatrixConstructorToTable_;

    lduMatrix::smoother::addasymMatrixConstructorToTable<symGaussSeidelSmoother>
        addsymGaussSeidelSmootherAsymMatrixConstructorToTable_;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::symGaussSeidelSmoother::symGaussSeidelSmoother
(
    const word& fieldName,
    const lduMatrix& matrix,
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const FieldField<Field, scalar>& interfaceIntCoeffs,
    const lduInterfaceFieldPtrsList& interfaces
)
:
    lduMatrix::smoother
    (
        fieldName,
        matrix,
        interfaceBouCoeffs,
        interfaceIntCoeffs,
        interfaces
    )
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::symGaussSeidelSmoother::smooth
(
    const word& fieldName_,
73
    solveScalarField& psi,
74
    const lduMatrix& matrix_,
75
    const solveScalarField& source,
76 77 78 79 80 81
    const FieldField<Field, scalar>& interfaceBouCoeffs_,
    const lduInterfaceFieldPtrsList& interfaces_,
    const direction cmpt,
    const label nSweeps
)
{
82
    solveScalar* __restrict__ psiPtr = psi.begin();
83

84
    const label nCells = psi.size();
85

86 87
    solveScalarField bPrime(nCells);
    solveScalar* __restrict__ bPrimePtr = bPrime.begin();
88

89 90
    const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
    const scalar* const __restrict__ upperPtr =
91
        matrix_.upper().begin();
92
    const scalar* const __restrict__ lowerPtr =
93 94
        matrix_.lower().begin();

95
    const label* const __restrict__ uPtr =
96 97
        matrix_.lduAddr().upperAddr().begin();

98
    const label* const __restrict__ ownStartPtr =
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
        matrix_.lduAddr().ownerStartAddr().begin();


    // Parallel boundary initialisation.  The parallel boundary is treated
    // as an effective jacobi interface in the boundary.
    // Note: there is a change of sign in the coupled
    // interface update.  The reason for this is that the
    // internal coefficients are all located at the l.h.s. of
    // the matrix whereas the "implicit" coefficients on the
    // coupled boundaries are all created as if the
    // coefficient contribution is of a source-kind (i.e. they
    // have a sign as if they are on the r.h.s. of the matrix.
    // To compensate for this, it is necessary to turn the
    // sign of the contribution.

    for (label sweep=0; sweep<nSweeps; sweep++)
    {
        bPrime = source;

118 119
        const label startRequest = Pstream::nRequests();

120 121
        matrix_.initMatrixInterfaces
        (
122 123
            false,
            interfaceBouCoeffs_,
124 125 126 127 128 129 130 131
            interfaces_,
            psi,
            bPrime,
            cmpt
        );

        matrix_.updateMatrixInterfaces
        (
132 133
            false,
            interfaceBouCoeffs_,
134 135 136
            interfaces_,
            psi,
            bPrime,
137 138
            cmpt,
            startRequest
139 140
        );

141
        solveScalar psii;
142 143
        label fStart;
        label fEnd = ownStartPtr[0];
144

145
        for (label celli=0; celli<nCells; celli++)
146 147 148 149 150 151 152 153 154
        {
            // Start and end of this row
            fStart = fEnd;
            fEnd = ownStartPtr[celli + 1];

            // Get the accumulated neighbour side
            psii = bPrimePtr[celli];

            // Accumulate the owner product side
155
            for (label facei=fStart; facei<fEnd; facei++)
156 157 158 159 160 161 162 163
            {
                psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
            }

            // Finish current psi
            psii /= diagPtr[celli];

            // Distribute the neighbour side using current psi
164
            for (label facei=fStart; facei<fEnd; facei++)
165 166 167 168 169 170 171 172 173
            {
                bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
            }

            psiPtr[celli] = psii;
        }

        fStart = ownStartPtr[nCells];

174
        for (label celli=nCells-1; celli>=0; celli--)
175 176 177 178 179 180 181 182 183
        {
            // Start and end of this row
            fEnd = fStart;
            fStart = ownStartPtr[celli];

            // Get the accumulated neighbour side
            psii = bPrimePtr[celli];

            // Accumulate the owner product side
184
            for (label facei=fStart; facei<fEnd; facei++)
185 186 187 188 189 190 191 192
            {
                psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
            }

            // Finish psi for this cell
            psii /= diagPtr[celli];

            // Distribute the neighbour side using psi for this cell
193
            for (label facei=fStart; facei<fEnd; facei++)
194 195 196 197 198 199 200 201 202 203
            {
                bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
            }

            psiPtr[celli] = psii;
        }
    }
}


204
void Foam::symGaussSeidelSmoother::scalarSmooth
205
(
206 207
    solveScalarField& psi,
    const solveScalarField& source,
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
    const direction cmpt,
    const label nSweeps
) const
{
    smooth
    (
        fieldName_,
        psi,
        matrix_,
        source,
        interfaceBouCoeffs_,
        interfaces_,
        cmpt,
        nSweeps
    );
}


226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
void Foam::symGaussSeidelSmoother::smooth
(
    solveScalarField& psi,
    const scalarField& source,
    const direction cmpt,
    const label nSweeps
) const
{
    scalarSmooth
    (
        psi,
        ConstPrecisionAdaptor<solveScalar, scalar>(source),
        cmpt,
        nSweeps
    );
}


244
// ************************************************************************* //