lduMatrixATmul.C 9.28 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) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2017-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 30 31 32 33 34 35 36 37 38

Description
    Multiply a given vector (second argument) by the matrix or its transpose
    and return the result in the first argument.

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

#include "lduMatrix.H"

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

void Foam::lduMatrix::Amul
(
39 40
    solveScalarField& Apsi,
    const tmp<solveScalarField>& tpsi,
41 42 43 44 45
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const lduInterfaceFieldPtrsList& interfaces,
    const direction cmpt
) const
{
46
    solveScalar* __restrict__ ApsiPtr = Apsi.begin();
47

48 49
    const solveScalarField& psi = tpsi();
    const solveScalar* const __restrict__ psiPtr = psi.begin();
50 51 52 53 54 55 56 57 58

    const scalar* const __restrict__ diagPtr = diag().begin();

    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();

    const scalar* const __restrict__ upperPtr = upper().begin();
    const scalar* const __restrict__ lowerPtr = lower().begin();

59 60
    const label startRequest = Pstream::nRequests();

61 62 63
    // Initialise the update of interfaced interfaces
    initMatrixInterfaces
    (
64
        true,
65 66 67
        interfaceBouCoeffs,
        interfaces,
        psi,
henry's avatar
henry committed
68
        Apsi,
69 70 71
        cmpt
    );

72 73
    const label nCells = diag().size();
    for (label cell=0; cell<nCells; cell++)
74 75 76 77 78
    {
        ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
    }


79
    const label nFaces = upper().size();
henry's avatar
henry committed
80

81
    for (label face=0; face<nFaces; face++)
82 83 84 85 86 87 88 89
    {
        ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
        ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
    }

    // Update interface interfaces
    updateMatrixInterfaces
    (
90
        true,
91 92 93 94
        interfaceBouCoeffs,
        interfaces,
        psi,
        Apsi,
95 96
        cmpt,
        startRequest
97 98 99 100 101 102 103 104
    );

    tpsi.clear();
}


void Foam::lduMatrix::Tmul
(
105 106
    solveScalarField& Tpsi,
    const tmp<solveScalarField>& tpsi,
107 108 109 110 111
    const FieldField<Field, scalar>& interfaceIntCoeffs,
    const lduInterfaceFieldPtrsList& interfaces,
    const direction cmpt
) const
{
112
    solveScalar* __restrict__ TpsiPtr = Tpsi.begin();
113

114 115
    const solveScalarField& psi = tpsi();
    const solveScalar* const __restrict__ psiPtr = psi.begin();
116 117 118 119 120 121 122 123 124

    const scalar* const __restrict__ diagPtr = diag().begin();

    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();

    const scalar* const __restrict__ lowerPtr = lower().begin();
    const scalar* const __restrict__ upperPtr = upper().begin();

125 126
    const label startRequest = Pstream::nRequests();

127 128 129
    // Initialise the update of interfaced interfaces
    initMatrixInterfaces
    (
130
        true,
131 132 133 134 135 136 137
        interfaceIntCoeffs,
        interfaces,
        psi,
        Tpsi,
        cmpt
    );

138 139
    const label nCells = diag().size();
    for (label cell=0; cell<nCells; cell++)
140 141 142 143
    {
        TpsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
    }

144 145
    const label nFaces = upper().size();
    for (label face=0; face<nFaces; face++)
146 147 148 149 150 151 152 153
    {
        TpsiPtr[uPtr[face]] += upperPtr[face]*psiPtr[lPtr[face]];
        TpsiPtr[lPtr[face]] += lowerPtr[face]*psiPtr[uPtr[face]];
    }

    // Update interface interfaces
    updateMatrixInterfaces
    (
154
        true,
155 156 157 158
        interfaceIntCoeffs,
        interfaces,
        psi,
        Tpsi,
159 160
        cmpt,
        startRequest
161 162 163 164 165 166 167 168
    );

    tpsi.clear();
}


void Foam::lduMatrix::sumA
(
169
    solveScalarField& sumA,
170 171 172 173
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const lduInterfaceFieldPtrsList& interfaces
) const
{
174
    solveScalar* __restrict__ sumAPtr = sumA.begin();
175 176 177 178 179 180 181 182 183

    const scalar* __restrict__ diagPtr = diag().begin();

    const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
    const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();

    const scalar* __restrict__ lowerPtr = lower().begin();
    const scalar* __restrict__ upperPtr = upper().begin();

184 185
    const label nCells = diag().size();
    const label nFaces = upper().size();
186

187
    for (label cell=0; cell<nCells; cell++)
188 189 190 191
    {
        sumAPtr[cell] = diagPtr[cell];
    }

192
    for (label face=0; face<nFaces; face++)
193 194 195 196 197 198 199
    {
        sumAPtr[uPtr[face]] += lowerPtr[face];
        sumAPtr[lPtr[face]] += upperPtr[face];
    }

    // Add the interface internal coefficients to diagonal
    // and the interface boundary coefficients to the sum-off-diagonal
200
    forAll(interfaces, patchi)
201
    {
202
        if (interfaces.set(patchi))
203
        {
204 205
            const labelUList& pa = lduAddr().patchAddr(patchi);
            const scalarField& pCoeffs = interfaceBouCoeffs[patchi];
206 207 208 209 210 211 212 213 214 215 216 217

            forAll(pa, face)
            {
                sumAPtr[pa[face]] -= pCoeffs[face];
            }
        }
    }
}


void Foam::lduMatrix::residual
(
218 219
    solveScalarField& rA,
    const solveScalarField& psi,
220 221 222 223 224 225
    const scalarField& source,
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const lduInterfaceFieldPtrsList& interfaces,
    const direction cmpt
) const
{
226
    solveScalar* __restrict__ rAPtr = rA.begin();
227

228
    const solveScalar* const __restrict__ psiPtr = psi.begin();
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
    const scalar* const __restrict__ diagPtr = diag().begin();
    const scalar* const __restrict__ sourcePtr = source.begin();

    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();

    const scalar* const __restrict__ upperPtr = upper().begin();
    const scalar* const __restrict__ lowerPtr = lower().begin();

    // Parallel boundary initialisation.
    // 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.

249 250
    const label startRequest = Pstream::nRequests();

251 252 253
    // Initialise the update of interfaced interfaces
    initMatrixInterfaces
    (
254 255
        false,
        interfaceBouCoeffs,
256 257
        interfaces,
        psi,
henry's avatar
henry committed
258
        rA,
259 260 261
        cmpt
    );

262 263
    const label nCells = diag().size();
    for (label cell=0; cell<nCells; cell++)
264 265 266 267 268
    {
        rAPtr[cell] = sourcePtr[cell] - diagPtr[cell]*psiPtr[cell];
    }


269
    const label nFaces = upper().size();
henry's avatar
henry committed
270

271
    for (label face=0; face<nFaces; face++)
272 273 274 275 276 277 278 279
    {
        rAPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
        rAPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
    }

    // Update interface interfaces
    updateMatrixInterfaces
    (
280 281
        false,
        interfaceBouCoeffs,
282 283 284
        interfaces,
        psi,
        rA,
285 286
        cmpt,
        startRequest
287 288 289 290
    );
}


291
Foam::tmp<Foam::Field<Foam::solveScalar>> Foam::lduMatrix::residual
292
(
293
    const solveScalarField& psi,
294 295 296 297 298 299
    const scalarField& source,
    const FieldField<Field, scalar>& interfaceBouCoeffs,
    const lduInterfaceFieldPtrsList& interfaces,
    const direction cmpt
) const
{
300
    tmp<solveScalarField> trA(new solveScalarField(psi.size()));
301
    residual(trA.ref(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
302 303 304 305
    return trA;
}


306
Foam::tmp<Foam::scalarField> Foam::lduMatrix::H1() const
307
{
308
    auto tH1 = tmp<scalarField>::New(lduAddr().size(), Zero);
309 310 311

    if (lowerPtr_ || upperPtr_)
    {
312
        scalar* __restrict__ H1Ptr = tH1.ref().begin();
313 314 315 316 317 318 319

        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();

        const scalar* __restrict__ lowerPtr = lower().begin();
        const scalar* __restrict__ upperPtr = upper().begin();

320
        const label nFaces = upper().size();
321

322
        for (label face=0; face<nFaces; face++)
323 324 325 326 327 328 329 330 331 332
        {
            H1Ptr[uPtr[face]] -= lowerPtr[face];
            H1Ptr[lPtr[face]] -= upperPtr[face];
        }
    }

    return tH1;
}


333
// ************************************************************************* //