From 3e39990c9fa77db3acfd79cc0b73380b12efd670 Mon Sep 17 00:00:00 2001 From: Henry <Henry> Date: Thu, 17 Jan 2013 09:50:55 +0000 Subject: [PATCH] GAMG: Add correction interpolation after prolongation --- .../solvers/GAMG/GAMGSolverInterpolate.C | 86 +++++++++++++++++++ .../lduMatrix/solvers/GAMG/GAMGSolverScale.C | 78 +++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C new file mode 100644 index 00000000000..d7fe071922d --- /dev/null +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "GAMGSolver.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::GAMGSolver::interpolate +( + scalarField& psi, + scalarField& Apsi, + const lduMatrix& m, + const FieldField<Field, scalar>& interfaceBouCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + const scalarField& source, + const direction cmpt +) const +{ + scalar* __restrict__ psiPtr = psi.begin(); + + const label* const __restrict__ uPtr = m.lduAddr().upperAddr().begin(); + const label* const __restrict__ lPtr = m.lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ diagPtr = m.diag().begin(); + const scalar* const __restrict__ upperPtr = m.upper().begin(); + const scalar* const __restrict__ lowerPtr = m.lower().begin(); + + Apsi = 0; + scalar* __restrict__ ApsiPtr = Apsi.begin(); + + m.initMatrixInterfaces + ( + interfaceBouCoeffs, + interfaces, + psi, + Apsi, + cmpt + ); + + register const label nFaces = m.upper().size(); + for (register label face=0; face<nFaces; face++) + { + ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]]; + ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]]; + } + + m.updateMatrixInterfaces + ( + interfaceBouCoeffs, + interfaces, + psi, + Apsi, + cmpt + ); + + register const label nCells = m.diag().size(); + for (register label celli=0; celli<nCells; celli++) + { + psiPtr[celli] = -ApsiPtr[celli]/(diagPtr[celli]); + } +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C new file mode 100644 index 00000000000..7629a28d665 --- /dev/null +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C @@ -0,0 +1,78 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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 "GAMGSolver.H" +#include "vector2D.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::GAMGSolver::scale +( + scalarField& field, + scalarField& Acf, + const lduMatrix& A, + const FieldField<Field, scalar>& interfaceLevelBouCoeffs, + const lduInterfaceFieldPtrsList& interfaceLevel, + const scalarField& source, + const direction cmpt +) const +{ + A.Amul + ( + Acf, + field, + interfaceLevelBouCoeffs, + interfaceLevel, + cmpt + ); + + scalar scalingFactorNum = 0.0; + scalar scalingFactorDenom = 0.0; + + forAll(field, i) + { + scalingFactorNum += source[i]*field[i]; + scalingFactorDenom += Acf[i]*field[i]; + } + + vector2D scalingVector(scalingFactorNum, scalingFactorDenom); + reduce(scalingVector, sumOp<vector2D>()); + scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL); + + if (debug >= 2) + { + Pout<< sf << " "; + } + + const scalarField& D = A.diag(); + + forAll(field, i) + { + field[i] = sf*field[i] + (source[i] - sf*Acf[i])/D[i]; + } +} + + +// ************************************************************************* // -- GitLab