From f8d14b0b289d653f5be64f2b7b99487e1b0e869c Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Sun, 16 Dec 2012 22:29:33 +0000
Subject: [PATCH] symGaussSeidelSmoother: symmetric Gauss-Seidel smoother
---
.../symGaussSeidel/symGaussSeidelSmoother.C | 241 ++++++++++++++++++
.../symGaussSeidel/symGaussSeidelSmoother.H | 108 ++++++++
2 files changed, 349 insertions(+)
create mode 100644 src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
create mode 100644 src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.H
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C b/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
new file mode 100644
index 00000000000..c8794ba5d6b
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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 "symGaussSeidelSmoother.H"
+
+// * * * * * * * * * * * * * * 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_,
+ scalarField& psi,
+ const lduMatrix& matrix_,
+ const scalarField& source,
+ const FieldField<Field, scalar>& interfaceBouCoeffs_,
+ const lduInterfaceFieldPtrsList& interfaces_,
+ const direction cmpt,
+ const label nSweeps
+)
+{
+ register scalar* __restrict__ psiPtr = psi.begin();
+
+ register const label nCells = psi.size();
+
+ scalarField bPrime(nCells);
+ register scalar* __restrict__ bPrimePtr = bPrime.begin();
+
+ register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
+ register const scalar* const __restrict__ upperPtr =
+ matrix_.upper().begin();
+ register const scalar* const __restrict__ lowerPtr =
+ matrix_.lower().begin();
+
+ register const label* const __restrict__ uPtr =
+ matrix_.lduAddr().upperAddr().begin();
+
+ register const label* const __restrict__ ownStartPtr =
+ 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.
+
+ FieldField<Field, scalar>& mBouCoeffs =
+ const_cast<FieldField<Field, scalar>&>
+ (
+ interfaceBouCoeffs_
+ );
+
+ forAll(mBouCoeffs, patchi)
+ {
+ if (interfaces_.set(patchi))
+ {
+ mBouCoeffs[patchi].negate();
+ }
+ }
+
+
+ for (label sweep=0; sweep<nSweeps; sweep++)
+ {
+ bPrime = source;
+
+ matrix_.initMatrixInterfaces
+ (
+ mBouCoeffs,
+ interfaces_,
+ psi,
+ bPrime,
+ cmpt
+ );
+
+ matrix_.updateMatrixInterfaces
+ (
+ mBouCoeffs,
+ interfaces_,
+ psi,
+ bPrime,
+ cmpt
+ );
+
+ register scalar psii;
+ register label fStart;
+ register label fEnd = ownStartPtr[0];
+
+ for (register label celli=0; celli<nCells; celli++)
+ {
+ // 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
+ for (register label facei=fStart; facei<fEnd; facei++)
+ {
+ psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
+ }
+
+ // Finish current psi
+ psii /= diagPtr[celli];
+
+ // Distribute the neighbour side using current psi
+ for (register label facei=fStart; facei<fEnd; facei++)
+ {
+ bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
+ }
+
+ psiPtr[celli] = psii;
+ }
+
+ fStart = ownStartPtr[nCells];
+
+ for (register label celli=nCells-1; celli>=0; celli--)
+ {
+ // Start and end of this row
+ fEnd = fStart;
+ fStart = ownStartPtr[celli];
+
+ // Get the accumulated neighbour side
+ psii = bPrimePtr[celli];
+
+ // Accumulate the owner product side
+ for (register label facei=fStart; facei<fEnd; facei++)
+ {
+ psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
+ }
+
+ // Finish psi for this cell
+ psii /= diagPtr[celli];
+
+ // Distribute the neighbour side using psi for this cell
+ for (register label facei=fStart; facei<fEnd; facei++)
+ {
+ bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
+ }
+
+ psiPtr[celli] = psii;
+ }
+ }
+
+ // Restore interfaceBouCoeffs_
+ forAll(mBouCoeffs, patchi)
+ {
+ if (interfaces_.set(patchi))
+ {
+ mBouCoeffs[patchi].negate();
+ }
+ }
+}
+
+
+void Foam::symGaussSeidelSmoother::smooth
+(
+ scalarField& psi,
+ const scalarField& source,
+ const direction cmpt,
+ const label nSweeps
+) const
+{
+ smooth
+ (
+ fieldName_,
+ psi,
+ matrix_,
+ source,
+ interfaceBouCoeffs_,
+ interfaces_,
+ cmpt,
+ nSweeps
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.H b/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.H
new file mode 100644
index 00000000000..a9a25e651fb
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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/>.
+
+Class
+ Foam::symGaussSeidelSmoother
+
+Description
+ A lduMatrix::smoother for symmetric Gauss-Seidel
+
+SourceFiles
+ symGaussSeidelSmoother.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef symGaussSeidelSmoother_H
+#define symGaussSeidelSmoother_H
+
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class symGaussSeidelSmoother Declaration
+\*---------------------------------------------------------------------------*/
+
+class symGaussSeidelSmoother
+:
+ public lduMatrix::smoother
+{
+
+public:
+
+ //- Runtime type information
+ TypeName("symGaussSeidel");
+
+
+ // Constructors
+
+ //- Construct from components
+ symGaussSeidelSmoother
+ (
+ const word& fieldName,
+ const lduMatrix& matrix,
+ const FieldField<Field, scalar>& interfaceBouCoeffs,
+ const FieldField<Field, scalar>& interfaceIntCoeffs,
+ const lduInterfaceFieldPtrsList& interfaces
+ );
+
+
+ // Member Functions
+
+ //- Smooth for the given number of sweeps
+ static void smooth
+ (
+ const word& fieldName,
+ scalarField& psi,
+ const lduMatrix& matrix,
+ const scalarField& source,
+ const FieldField<Field, scalar>& interfaceBouCoeffs,
+ const lduInterfaceFieldPtrsList& interfaces,
+ const direction cmpt,
+ const label nSweeps
+ );
+
+
+ //- Smooth the solution for a given number of sweeps
+ virtual void smooth
+ (
+ scalarField& psi,
+ const scalarField& Source,
+ const direction cmpt,
+ const label nSweeps
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--
GitLab