zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Removed templated solvers. Will be added to next version.

Browse Source
This commit is contained in:
henry
2008-06-09 22:57:49 +01:00
parent 28d5c784dc
commit 107ff36066
40 changed files with 0 additions and 5927 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/OpenFOAM/Make/files
View File

@ -256,11 +256,6 @@ $(pairGAMGAgglomeration)/pairGAMGAgglomerationCombineLevels.C
algebraicPairGAMGAgglomeration = $(GAMGAgglomerations)/algebraicPairGAMGAgglomeration algebraicPairGAMGAgglomeration = $(GAMGAgglomerations)/algebraicPairGAMGAgglomeration
$(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C $(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C
matrices/LduMatrix/LduMatrix/lduMatrices.C
matrices/LduMatrix/Preconditioners/lduPreconditioners.C
matrices/LduMatrix/Smoothers/lduSmoothers.C
matrices/LduMatrix/Solvers/lduSolvers.C
meshes/lduMesh/lduMesh.C meshes/lduMesh/lduMesh.C
primitiveShapes = meshes/primitiveShapes primitiveShapes = meshes/primitiveShapes

66
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.C
View File

@ -1,66 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "CyclicLduInterfaceField.H"
#include "diagTensorField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//defineTypeNameAndDebug(CyclicLduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::CyclicLduInterfaceField<Type>::~CyclicLduInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::CyclicLduInterfaceField<Type>::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
label sizeby2 = f.size()/2;
for (label facei=0; facei<sizeby2; facei++)
{
f[facei] = transform(f[facei], forwardT()[0]);
f[facei + sizeby2] = transform(f[facei + sizeby2], forwardT()[0]);
}
}
}
// ************************************************************************* //

103
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.H
View File

@ -1,103 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::CyclicLduInterfaceField
Description
Abstract base class for cyclic coupled interfaces.
SourceFiles
CyclicLduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef CyclicLduInterfaceField_H
#define CyclicLduInterfaceField_H
#include "primitiveFieldsFwd.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class CyclicLduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class CyclicLduInterfaceField
{
public:
//- Runtime type information
TypeName("CyclicLduInterfaceField");
// Constructors
//- Construct given coupled patch
CyclicLduInterfaceField()
{}
// Destructor
virtual ~CyclicLduInterfaceField();
// Member Functions
// Access
//- Is the transform required
virtual bool doTransform() const = 0;
//- Return face transformation tensor
virtual const tensorField& forwardT() const = 0;
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const = 0;
//- Return rank of component for transform
virtual int rank() const = 0;
//- Transform given patch internal field
void transformCoupleField(Field<Type>& f) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

48
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.C
View File

@ -1,48 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "lduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(lduInterfaceField, 0);
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
lduInterfaceField::~lduInterfaceField()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

147
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H
View File

@ -1,147 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::LduInterfaceField
Description
An abstract base class for implicitly-coupled interface fields
e.g. processor and cyclic patch fields.
SourceFiles
LduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef LduInterfaceField_H
#define LduInterfaceField_H
#include "lduInterface.H"
#include "Field.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class LduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class LduInterfaceField
{
// Private data
//- Reference to the coupled patch this field is defined for
const lduInterface& interface_;
// Private Member Functions
//- Disallow default bitwise copy construct
LduInterfaceField(const LduInterfaceField&);
//- Disallow default bitwise assignment
void operator=(const LduInterfaceField&);
public:
class Amultiplier
{
public:
Amultiplier()
{}
virtual void addAmul
(
Field<Type>& Apsi,
const Field<Type>& psi
) const = 0;
};
//- Runtime type information
TypeName("LduInterfaceField");
// Constructors
//- Construct given interface
LduInterfaceField(const lduInterface& interface)
:
interface_(interface)
{}
// Destructor
virtual ~LduInterfaceField();
// Member Functions
// Access
//- Return the interface
const lduInterface& interface() const
{
return interface_;
}
// Coupled interface matrix update
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
Field<Type>& Apsi,
const Field<Type>& psi,
const Amultiplier&,
const Pstream::commsTypes commsType
) const
{}
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>& Apsi,
const Field<Type>& psi,
const Amultiplier&,
const Pstream::commsTypes commsType
) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

71
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceFieldPtrsList.H
View File

@ -1,71 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Type
lduInterfaceFieldPtrsList
Description
List of coupled interface fields to be used in coupling.
\*---------------------------------------------------------------------------*/
#ifndef LduInterfaceFieldPtrsList_H
#define LduInterfaceFieldPtrsList_H
#include "LduInterfaceField.H"
#include "UPtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class LduInterfaceFieldPtrsList Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class LduInterfaceFieldPtrsList
:
public UPtrList<const LduInterfaceField<Type> >
{
public:
LduInterfaceFieldPtrsList(label size)
:
UPtrList<const LduInterfaceField<Type> >(size)
{}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

67
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.C
View File

@ -1,67 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "processorLduInterfaceField.H"
#include "diagTensorField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//defineTypeNameAndDebug(processorLduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::processorLduInterfaceField<Type>::~processorLduInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::processorLduInterfaceField<Type>::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
if (forwardT().size() == 1)
{
transform(f, forwardT()[0], f);
}
else
{
transform(f, forwardT(), f);
}
}
}
// ************************************************************************* //

106
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.H
View File

@ -1,106 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::ProcessorLduInterfaceField
Description
Abstract base class for processor coupled interfaces.
SourceFiles
ProcessorLduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef ProcessorLduInterfaceField_H
#define ProcessorLduInterfaceField_H
#include "primitiveFieldsFwd.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ProcessorLduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class ProcessorLduInterfaceField
{
public:
//- Runtime type information
TypeName("ProcessorLduInterfaceField");
// Constructors
//- Construct given coupled patch
ProcessorLduInterfaceField()
{}
// Destructor
virtual ~ProcessorLduInterfaceField();
// Member Functions
// Access
//- Return processor number
virtual int myProcNo() const = 0;
//- Return neigbour processor number
virtual int neighbProcNo() const = 0;
//- Is the transform required
virtual bool doTransform() const = 0;
//- Return face transformation tensor
virtual const tensorField& forwardT() const = 0;
//- Return rank of component for transform
virtual int rank() const = 0;
//- Transform given patch component field
void transformCoupleField(Field<Type>& f) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

391
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
View File

@ -1,391 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "lduMatrix.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const lduMesh& mesh)
:
lduMesh_(mesh),
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
sourcePtr_(NULL),
interfaces_(0),
interfacesUpper_(0),
interfacesLower_(0)
{}
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const LduMatrix& A)
:
lduMesh_(A.lduMesh_),
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
sourcePtr_(NULL),
interfaces_(0),
interfacesUpper_(0),
interfacesLower_(0)
{
if (A.diagPtr_)
{
diagPtr_ = new Field<DType>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
upperPtr_ = new Field<LUType>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
}
if (A.sourcePtr_)
{
sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
}
}
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reUse)
:
lduMesh_(A.lduMesh_),
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
sourcePtr_(NULL),
interfaces_(0),
interfacesUpper_(0),
interfacesLower_(0)
{
if (reUse)
{
if (A.diagPtr_)
{
diagPtr_ = A.diagPtr_;
A.diagPtr_ = NULL;
}
if (A.upperPtr_)
{
upperPtr_ = A.upperPtr_;
A.upperPtr_ = NULL;
}
if (A.lowerPtr_)
{
lowerPtr_ = A.lowerPtr_;
A.lowerPtr_ = NULL;
}
if (A.sourcePtr_)
{
sourcePtr_ = A.sourcePtr_;
A.sourcePtr_ = NULL;
}
}
else
{
if (A.diagPtr_)
{
diagPtr_ = new Field<DType>(*(A.diagPtr_));
}
if (A.upperPtr_)
{
upperPtr_ = new Field<LUType>(*(A.upperPtr_));
}
if (A.lowerPtr_)
{
lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
}
if (A.sourcePtr_)
{
sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
}
}
}
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::LduMatrix
(
const lduMesh& mesh,
Istream& is
)
:
lduMesh_(mesh),
diagPtr_(new Field<DType>(is)),
upperPtr_(new Field<LUType>(is)),
lowerPtr_(new Field<LUType>(is)),
sourcePtr_(new Field<Type>(is)),
interfaces_(0),
interfacesUpper_(0),
interfacesLower_(0)
{}
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
{
if (diagPtr_)
{
delete diagPtr_;
}
if (upperPtr_)
{
delete upperPtr_;
}
if (lowerPtr_)
{
delete lowerPtr_;
}
if (sourcePtr_)
{
delete sourcePtr_;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
{
if (!diagPtr_)
{
diagPtr_ = new Field<DType>(lduAddr().size(), pTraits<DType>::zero);
}
return *diagPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
{
if (!upperPtr_)
{
if (lowerPtr_)
{
upperPtr_ = new Field<LUType>(*lowerPtr_);
}
else
{
upperPtr_ = new Field<LUType>
(
lduAddr().lowerAddr().size(),
pTraits<LUType>::zero
);
}
}
return *upperPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
{
if (!lowerPtr_)
{
if (upperPtr_)
{
lowerPtr_ = new Field<LUType>(*upperPtr_);
}
else
{
lowerPtr_ = new Field<LUType>
(
lduAddr().lowerAddr().size(),
pTraits<LUType>::zero
);
}
}
return *lowerPtr_;
}
template<class Type, class DType, class LUType>
Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
{
if (!sourcePtr_)
{
sourcePtr_ = new Field<Type>(lduAddr().size(), pTraits<Type>::zero);
}
return *sourcePtr_;
}
template<class Type, class DType, class LUType>
const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
{
if (!diagPtr_)
{
FatalErrorIn
(
"const Field<DType>& LduMatrix<Type, DType, LUType>::diag() const"
) << "diagPtr_ unallocated"
<< abort(FatalError);
}
return *diagPtr_;
}
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
{
if (!lowerPtr_ && !upperPtr_)
{
FatalErrorIn
(
"const Field<LUType>& LduMatrix<Type, DType, LUType>::upper() const"
) << "lowerPtr_ or upperPtr_ unallocated"
<< abort(FatalError);
}
if (upperPtr_)
{
return *upperPtr_;
}
else
{
return *lowerPtr_;
}
}
template<class Type, class DType, class LUType>
const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
{
if (!lowerPtr_ && !upperPtr_)
{
FatalErrorIn
(
"const Field<LUType>& LduMatrix<Type, DType, LUType>::lower() const"
) << "lowerPtr_ or upperPtr_ unallocated"
<< abort(FatalError);
}
if (lowerPtr_)
{
return *lowerPtr_;
}
else
{
return *upperPtr_;
}
}
template<class Type, class DType, class LUType>
const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
{
if (!sourcePtr_)
{
FatalErrorIn
(
"const Field<Type>& LduMatrix<Type, DType, LUType>::source() const"
) << "sourcePtr_ unallocated"
<< abort(FatalError);
}
return *sourcePtr_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const LduMatrix<Type, DType, LUType>& ldum
)
{
if (ldum.diagPtr_)
{
os << "Diagonal = "
<< *ldum.diagPtr_
<< endl << endl;
}
if (ldum.upperPtr_)
{
os << "Upper triangle = "
<< *ldum.upperPtr_
<< endl << endl;
}
if (ldum.lowerPtr_)
{
os << "Lower triangle = "
<< *ldum.lowerPtr_
<< endl << endl;
}
if (ldum.sourcePtr_)
{
os << "Source = "
<< *ldum.sourcePtr_
<< endl << endl;
}
os.check("Ostream& operator<<(Ostream&, const LduMatrix&");
return os;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "LduMatrixOperations.C"
#include "LduMatrixATmul.C"
#include "LduMatrixUpdateMatrixInterfaces.C"
#include "LduMatrixTests.C"
#include "LduMatrixPreconditioner.C"
#include "LduMatrixSmoother.C"
#include "LduMatrixSolver.C"
// ************************************************************************* //

935
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
View File

@ -1,935 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::LduMatrix
Description
LduMatrix is a general matrix class in which the coefficients are
stored as three arrays, one for the upper triangle, one for the
lower triangle and a third for the diagonal.
Addressing arrays must be supplied for the upper and lower triangles.
Note
It might be better if this class were organised as a hierachy starting
from an empty matrix, then deriving diagonal, symmetric and asymmetric
matrices.
SourceFiles
LduMatrixATmul.C
LduMatrix.C
LduMatrixOperations.C
LduMatrixSolver.C
LduMatrixPreconditioner.C
LduMatrixTests.C
LduMatrixUpdateMatrixInterfaces.C
\*---------------------------------------------------------------------------*/
#ifndef LduMatrix_H
#define LduMatrix_H
#include "lduMesh.H"
#include "Field.H"
#include "FieldField.H"
#include "LduInterfaceFieldPtrsList.H"
#include "typeInfo.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
template<class Type, class DType, class LUType>
class LduMatrix;
template<class Type, class DType, class LUType>
Ostream& operator<<
(
Ostream&,
const LduMatrix<Type, DType, LUType>&
);
/*---------------------------------------------------------------------------*\
Class LduMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class LduMatrix
{
// private data
//- LDU mesh reference
const lduMesh& lduMesh_;
//- Diagonal coefficients
Field<DType> *diagPtr_;
//- Off-diagonal coefficients
Field<LUType> *upperPtr_, *lowerPtr_;
//- Source
Field<Type> *sourcePtr_;
//- Field interfaces (processor patches etc.)
LduInterfaceFieldPtrsList<Type> interfaces_;
//- Off-diagonal coefficients for interfaces
FieldField<Field, LUType> interfacesUpper_, interfacesLower_;
public:
//- Class returned by the solver
// containing performance statistics
class solverPerformance
{
word solverName_;
word fieldName_;
Type initialResidual_;
Type finalResidual_;
label noIterations_;
bool converged_;
FixedList<bool, pTraits<Type>::nComponents> singular_;
public:
// Constructors
solverPerformance()
:
initialResidual_(pTraits<Type>::zero),
finalResidual_(pTraits<Type>::zero),
noIterations_(0),
converged_(false),
singular_(false)
{}
solverPerformance
(
const word& solverName,
const word& fieldName,
const Type& iRes = pTraits<Type>::zero,
const Type& fRes = pTraits<Type>::zero,
const label nIter = 0,
const bool converged = false,
const bool singular = false
)
:
solverName_(solverName),
fieldName_(fieldName),
initialResidual_(iRes),
finalResidual_(fRes),
noIterations_(nIter),
converged_(converged),
singular_(singular)
{}
// Member functions
//- Return solver name
const word& solverName() const
{
return solverName_;
}
//- Return initial residual
const Type& initialResidual() const
{
return initialResidual_;
}
//- Return initial residual
Type& initialResidual()
{
return initialResidual_;
}
//- Return final residual
const Type& finalResidual() const
{
return finalResidual_;
}
//- Return final residual
Type& finalResidual()
{
return finalResidual_;
}
//- Return number of iterations
label nIterations() const
{
return noIterations_;
}
//- Return number of iterations
label& nIterations()
{
return noIterations_;
}
//- Check, store and return singularity
bool singular(const Type& wApA);
//- Is the matrix singular?
bool singular() const;
//- Check, store and return convergence
bool converged
(
const Type& tolerance,
const Type& relTolerance
);
//- Has the solver converged?
bool converged() const
{
return converged_;
}
//- Print summary of solver performance to the given stream
void print(Ostream& os) const;
};
//- Abstract base-class for LduMatrix solvers
class solver
{
protected:
// Protected data
word fieldName_;
const LduMatrix<Type, DType, LUType>& matrix_;
//- dictionary of controls
dictionary controlDict_;
//- Maximum number of iterations in the solver
label maxIter_;
//- Final convergence tolerance
Type tolerance_;
//- Convergence tolerance relative to the initial
Type relTol_;
// Protected Member Functions
//- Read a control parameter from controlDict
template<class T>
inline void readControl
(
const dictionary& controlDict,
T& control,
const word& controlName
);
//- Read the control parameters from the controlDict_
virtual void readControls();
public:
//- Runtime type information
virtual const word& type() const = 0;
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
solver,
symMatrix,
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
),
(
fieldName,
matrix,
solverDict
)
);
declareRunTimeSelectionTable
(
autoPtr,
solver,
asymMatrix,
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
),
(
fieldName,
matrix,
solverDict
)
);
// Constructors
solver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Selectors
//- Return a new solver
static autoPtr<solver> New
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Destructor
virtual ~solver()
{}
// Member functions
// Access
const word& fieldName() const
{
return fieldName_;
}
const LduMatrix<Type, DType, LUType>& matrix() const
{
return matrix_;
}
//- Read and reset the solver parameters from the given dictionary
virtual void read(const dictionary& solverDict);
virtual solverPerformance solve
(
Field<Type>& psi
) const = 0;
//- Return the matrix norm used to normalise the residual for the
// stopping criterion
Type normFactor
(
const Field<Type>& psi,
const Field<Type>& Apsi,
Field<Type>& tmpField
) const;
};
//- Abstract base-class for LduMatrix smoothers
class smoother
{
protected:
// Protected data
word fieldName_;
const LduMatrix<Type, DType, LUType>& matrix_;
public:
//- Runtime type information
virtual const word& type() const = 0;
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
smoother,
symMatrix,
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
),
(
fieldName,
matrix
)
);
declareRunTimeSelectionTable
(
autoPtr,
smoother,
asymMatrix,
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
),
(
fieldName,
matrix
)
);
// Constructors
smoother
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
);
// Selectors
//- Return a new smoother
static autoPtr<smoother> New
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& smootherDict
);
// Destructor
virtual ~smoother()
{}
// Member functions
// Access
const word& fieldName() const
{
return fieldName_;
}
const LduMatrix<Type, DType, LUType>& matrix() const
{
return matrix_;
}
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
Field<Type>& psi,
const label nSweeps
) const = 0;
};
//- Abstract base-class for LduMatrix preconditioners
class preconditioner
{
protected:
// Protected data
//- Reference to the base-solver this preconditioner is used with
const solver& solver_;
public:
//- Runtime type information
virtual const word& type() const = 0;
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
preconditioner,
symMatrix,
(
const solver& sol,
const dictionary& preconditionerDict
),
(sol, preconditionerDict)
);
declareRunTimeSelectionTable
(
autoPtr,
preconditioner,
asymMatrix,
(
const solver& sol,
const dictionary& preconditionerDict
),
(sol, preconditionerDict)
);
// Constructors
preconditioner
(
const solver& sol
)
:
solver_(sol)
{}
// Selectors
//- Return a new preconditioner
static autoPtr<preconditioner> New
(
const solver& sol,
const dictionary& preconditionerDict
);
// Destructor
virtual ~preconditioner()
{}
// Member functions
//- Read and reset the preconditioner parameters
// from the given dictionary
virtual void read(const dictionary& preconditionerDict)
{}
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const = 0;
//- Return wT the transpose-matrix preconditioned form of
// residual rT.
// This is only required for preconditioning asymmetric matrices.
virtual void preconditionT
(
Field<Type>& wT,
const Field<Type>& rT
) const
{
notImplemented
(
type() +"::preconditionT"
"(Field<Type>& wT, const Field<Type>& rT)"
);
}
};
// Static data
// Declare name of the class and its debug switch
ClassName("LduMatrix");
//- Large Type for the use in solvers
static const scalar great_;
//- Small Type for the use in solvers
static const scalar small_;
//- Very small Type for the use in solvers
static const scalar vsmall_;
// Constructors
//- Construct given an LDU addressed mesh.
// The coefficients are initially empty for subsequent setting.
LduMatrix(const lduMesh&);
//- Construct as copy
LduMatrix(const LduMatrix<Type, DType, LUType>&);
//- Construct as copy or re-use as specified.
LduMatrix(LduMatrix<Type, DType, LUType>&, bool reUse);
//- Construct given an LDU addressed mesh and an Istream
// from which the coefficients are read
LduMatrix(const lduMesh&, Istream&);
// Destructor
~LduMatrix();
// Member functions
// Access to addressing
//- Return the LDU mesh from which the addressing is obtained
const lduMesh& mesh() const
{
return lduMesh_;
}
//- Return the LDU addressing
const lduAddressing& lduAddr() const
{
return lduMesh_.lduAddr();
}
//- Return the patch evaluation schedule
const lduSchedule& patchSchedule() const
{
return lduAddr().patchSchedule();
}
//- Return interfaces
const LduInterfaceFieldPtrsList<Type>& interfaces() const
{
return interfaces_;
}
// Access to coefficients
Field<DType>& diag();
Field<LUType>& upper();
Field<LUType>& lower();
Field<Type>& source();
FieldField<Field, LUType>& interfacesUpper()
{
return interfacesUpper_;
}
FieldField<Field, LUType>& interfacesLower()
{
return interfacesLower_;
}
const Field<DType>& diag() const;
const Field<LUType>& upper() const;
const Field<LUType>& lower() const;
const Field<Type>& source() const;
const FieldField<Field, LUType>& interfacesUpper() const
{
return interfacesUpper_;
}
const FieldField<Field, LUType>& interfacesLower() const
{
return interfacesLower_;
}
bool hasDiag() const
{
return (diagPtr_);
}
bool hasUpper() const
{
return (upperPtr_);
}
bool hasLower() const
{
return (lowerPtr_);
}
bool hasSource() const
{
return (sourcePtr_);
}
bool diagonal() const
{
return (diagPtr_ && !lowerPtr_ && !upperPtr_);
}
bool symmetric() const
{
return (diagPtr_ && (!lowerPtr_ && upperPtr_));
}
bool asymmetric() const
{
return (diagPtr_ && lowerPtr_ && upperPtr_);
}
// operations
void sumDiag();
void negSumDiag();
void sumMagOffDiag(Field<LUType>& sumOff) const;
//- Matrix multiplication
void Amul(Field<Type>&, const tmp<Field<Type> >&) const;
//- Matrix transpose multiplication
void Tmul(Field<Type>&, const tmp<Field<Type> >&) const;
//- Sum the coefficients on each row of the matrix
void sumA(Field<Type>&) const;
void residual(Field<Type>& rA, const Field<Type>& psi) const;
tmp<Field<Type> > residual(const Field<Type>& psi) const;
//- Initialise the update of interfaced interfaces
// for matrix operations
void initMatrixInterfaces
(
const FieldField<Field, LUType>& interfaceCoeffs,
const Field<Type>& psiif,
Field<Type>& result
) const;
//- Update interfaced interfaces for matrix operations
void updateMatrixInterfaces
(
const FieldField<Field, LUType>& interfaceCoeffs,
const Field<Type>& psiif,
Field<Type>& result
) const;
tmp<Field<Type> > H(const Field<Type>&) const;
tmp<Field<Type> > H(const tmp<Field<Type> >&) const;
tmp<Field<Type> > faceH(const Field<Type>&) const;
tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
// Member operators
void operator=(const LduMatrix<Type, DType, LUType>&);
void negate();
void operator+=(const LduMatrix<Type, DType, LUType>&);
void operator-=(const LduMatrix<Type, DType, LUType>&);
void operator*=(const scalarField&);
void operator*=(scalar);
// Ostream operator
friend Ostream& operator<< <Type, DType, LUType>
(
Ostream&,
const LduMatrix<Type, DType, LUType>&
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeLduMatrix(Type, DType, LUType) \
\
typedef Foam::LduMatrix<Type, DType, LUType> \
ldu##Type##DType##LUType##Matrix; \
\
defineNamedTemplateTypeNameAndDebug(ldu##Type##DType##LUType##Matrix, 0); \
\
\
template<> \
const scalar ldu##Type##DType##LUType##Matrix::great_(1e15); \
\
template<> \
const scalar ldu##Type##DType##LUType##Matrix::small_(1e-15); \
\
template<> \
const scalar ldu##Type##DType##LUType##Matrix::vsmall_(VSMALL); \
\
\
typedef LduMatrix<Type, DType, LUType>::smoother \
ldu##Type##DType##LUType##Smoother; \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Smoother, \
symMatrix \
); \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Smoother, \
asymMatrix \
); \
\
\
typedef LduMatrix<Type, DType, LUType>::preconditioner \
ldu##Type##DType##LUType##Preconditioner; \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Preconditioner, \
symMatrix \
); \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Preconditioner, \
asymMatrix \
); \
\
\
typedef LduMatrix<Type, DType, LUType>::solver \
ldu##Type##DType##LUType##Solver; \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Solver, \
symMatrix \
); \
\
defineTemplateRunTimeSelectionTable \
( \
ldu##Type##DType##LUType##Solver, \
asymMatrix \
);
#define makeLduPreconditioner(Precon, Type, DType, LUType) \
\
typedef Precon<Type, DType, LUType> \
Precon##Type##DType##LUType##Preconditioner; \
defineNamedTemplateTypeNameAndDebug \
( \
Precon##Type##DType##LUType##Preconditioner, \
0 \
);
#define makeLduSymPreconditioner(Precon, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::preconditioner:: \
addsymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner> \
add##Precon##Type##DType##LUType##PreconditionerSymMatrixConstructorToTable_;
#define makeLduAsymPreconditioner(Precon, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::preconditioner:: \
addasymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner> \
add##Precon##Type##DType##LUType##PreconditionerAsymMatrixConstructorToTable_;
#define makeLduSmoother(Smoother, Type, DType, LUType) \
\
typedef Smoother<Type, DType, LUType> \
Smoother##Type##DType##LUType##Smoother; \
\
defineNamedTemplateTypeNameAndDebug \
( \
Smoother##Type##DType##LUType##Smoother, \
0 \
);
#define makeLduSymSmoother(Smoother, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::smoother:: \
addsymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother> \
add##Smoother##Type##DType##LUType##SymMatrixConstructorToTable_;
#define makeLduAsymSmoother(Smoother, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::smoother:: \
addasymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother> \
add##Smoother##Type##DType##LUType##AsymMatrixConstructorToTable_;
#define makeLduSolver(Solver, Type, DType, LUType) \
\
typedef Solver<Type, DType, LUType> \
Solver##Type##DType##LUType##Solver; \
\
defineNamedTemplateTypeNameAndDebug \
( \
Solver##Type##DType##LUType##Solver, \
0 \
);
#define makeLduSymSolver(Solver, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::solver:: \
addsymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver> \
add##Solver##Type##DType##LUType##SymMatrixConstructorToTable_;
#define makeLduAsymSolver(Solver, Type, DType, LUType) \
\
LduMatrix<Type, DType, LUType>::solver:: \
addasymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver> \
add##Solver##Type##DType##LUType##AsymMatrixConstructorToTable_;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "LduMatrixI.H"
# include "LduMatrix.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

291
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixATmul.C
View File

@ -1,291 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type, class LUType>
class Amultiplier
:
public LduInterfaceField<Type>::Amultiplier
{
const Field<LUType>& A_;
public:
Amultiplier(const Field<LUType>& A)
:
A_(A)
{}
virtual void addAmul(Field<Type>& Apsi, const Field<Type>& psi) const
{
Apsi += A_*psi;
}
};
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::Amul
(
Field<Type>& Apsi,
const tmp<Field<Type> >& tpsi
) const
{
Type* __restrict__ ApsiPtr = Apsi.begin();
const Field<Type>& psi = tpsi();
const Type* const __restrict__ psiPtr = psi.begin();
const DType* const __restrict__ diagPtr = diag().begin();
const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
const LUType* const __restrict__ upperPtr = upper().begin();
const LUType* const __restrict__ lowerPtr = lower().begin();
// Initialise the update of interfaced interfaces
initMatrixInterfaces
(
interfacesUpper_,
psi,
Apsi
);
register const label nCells = diag().size();
for (register label cell=0; cell<nCells; cell++)
{
ApsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
}
register const label nFaces = upper().size();
for (register label face=0; face<nFaces; face++)
{
ApsiPtr[uPtr[face]] += dot(lowerPtr[face], psiPtr[lPtr[face]]);
ApsiPtr[lPtr[face]] += dot(upperPtr[face], psiPtr[uPtr[face]]);
}
// Update interface interfaces
updateMatrixInterfaces
(
interfacesUpper_,
psi,
Apsi
);
tpsi.clear();
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::Tmul
(
Field<Type>& Tpsi,
const tmp<Field<Type> >& tpsi
) const
{
Type* __restrict__ TpsiPtr = Tpsi.begin();
const Field<Type>& psi = tpsi();
const Type* const __restrict__ psiPtr = psi.begin();
const DType* const __restrict__ diagPtr = diag().begin();
const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
const LUType* const __restrict__ lowerPtr = lower().begin();
const LUType* const __restrict__ upperPtr = upper().begin();
// Initialise the update of interfaced interfaces
initMatrixInterfaces
(
interfacesLower_,
psi,
Tpsi
);
register const label nCells = diag().size();
for (register label cell=0; cell<nCells; cell++)
{
TpsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
}
register const label nFaces = upper().size();
for (register label face=0; face<nFaces; face++)
{
TpsiPtr[uPtr[face]] += dot(upperPtr[face], psiPtr[lPtr[face]]);
TpsiPtr[lPtr[face]] += dot(lowerPtr[face], psiPtr[uPtr[face]]);
}
// Update interface interfaces
updateMatrixInterfaces
(
interfacesLower_,
psi,
Tpsi
);
tpsi.clear();
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::sumA
(
Field<Type>& sumA
) const
{
Type* __restrict__ sumAPtr = sumA.begin();
const DType* __restrict__ diagPtr = diag().begin();
const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
const LUType* __restrict__ lowerPtr = lower().begin();
const LUType* __restrict__ upperPtr = upper().begin();
register const label nCells = diag().size();
register const label nFaces = upper().size();
for (register label cell=0; cell<nCells; cell++)
{
sumAPtr[cell] = dot(diagPtr[cell], pTraits<Type>::one);
}
for (register label face=0; face<nFaces; face++)
{
sumAPtr[uPtr[face]] += dot(lowerPtr[face], pTraits<Type>::one);
sumAPtr[lPtr[face]] += dot(upperPtr[face], pTraits<Type>::one);
}
// Add the interface internal coefficients to diagonal
// and the interface boundary coefficients to the sum-off-diagonal
forAll(interfaces_, patchI)
{
if (interfaces_.set(patchI))
{
const unallocLabelList& pa = lduAddr().patchAddr(patchI);
const Field<LUType>& pCoeffs = interfacesUpper_[patchI];
forAll(pa, face)
{
sumAPtr[pa[face]] -= dot(pCoeffs[face], pTraits<Type>::one);
}
}
}
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::residual
(
Field<Type>& rA,
const Field<Type>& psi
) const
{
Type* __restrict__ rAPtr = rA.begin();
const Type* const __restrict__ psiPtr = psi.begin();
const DType* const __restrict__ diagPtr = diag().begin();
const Type* const __restrict__ sourcePtr = source().begin();
const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
const LUType* const __restrict__ upperPtr = upper().begin();
const LUType* const __restrict__ lowerPtr = lower().begin();
// Parallel boundary initialisation.
// Note: there is a change of sign in the coupled
// interface update to add the contibution to the r.h.s.
FieldField<Field, LUType> mBouCoeffs(interfacesUpper_.size());
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs.set(patchi, -interfacesUpper_[patchi]);
}
}
// Initialise the update of interfaced interfaces
initMatrixInterfaces
(
mBouCoeffs,
psi,
rA
);
register const label nCells = diag().size();
for (register label cell=0; cell<nCells; cell++)
{
rAPtr[cell] = sourcePtr[cell] - dot(diagPtr[cell], psiPtr[cell]);
}
register const label nFaces = upper().size();
for (register label face=0; face<nFaces; face++)
{
rAPtr[uPtr[face]] -= dot(lowerPtr[face], psiPtr[lPtr[face]]);
rAPtr[lPtr[face]] -= dot(upperPtr[face], psiPtr[uPtr[face]]);
}
// Update interface interfaces
updateMatrixInterfaces
(
mBouCoeffs,
psi,
rA
);
}
template<class Type, class DType, class LUType>
Foam::tmp<Foam::Field<Type> > Foam::LduMatrix<Type, DType, LUType>::residual
(
const Field<Type>& psi
) const
{
tmp<Field<Type> > trA(new Field<Type>(psi.size()));
residual(trA(), psi);
return trA;
}
// ************************************************************************* //

45
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixI.H
View File

@ -1,45 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
template<class T>
inline void Foam::LduMatrix<Type, DType, LUType>::solver::readControl
(
const dictionary& controlDict,
T& control,
const word& controlName
)
{
if (controlDict.found(controlName))
{
controlDict.lookup(controlName) >> control;
}
}
// ************************************************************************* //

481
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
View File

@ -1,481 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "lduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::sumDiag()
{
const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
Field<DType>& Diag = diag();
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();
for (register label face=0; face<l.size(); face++)
{
Diag[l[face]] += Lower[face];
Diag[u[face]] += Upper[face];
}
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::negSumDiag()
{
const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
Field<DType>& Diag = diag();
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();
for (register label face=0; face<l.size(); face++)
{
Diag[l[face]] -= Lower[face];
Diag[u[face]] -= Upper[face];
}
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::sumMagOffDiag
(
Field<LUType>& sumOff
) const
{
const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();
for (register label face = 0; face < l.size(); face++)
{
sumOff[u[face]] += cmptMag(Lower[face]);
sumOff[l[face]] += cmptMag(Upper[face]);
}
}
template<class Type, class DType, class LUType>
Foam::tmp<Foam::Field<Type> >
Foam::LduMatrix<Type, DType, LUType>::H(const Field<Type>& psi) const
{
tmp<Field<Type> > tHpsi
(
new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
);
if (lowerPtr_ || upperPtr_)
{
Field<Type> & Hpsi = tHpsi();
Type* __restrict__ HpsiPtr = Hpsi.begin();
const Type* __restrict__ psiPtr = psi.begin();
const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
const LUType* __restrict__ lowerPtr = lower().begin();
const LUType* __restrict__ upperPtr = upper().begin();
register const label nFaces = upper().size();
for (register label face=0; face<nFaces; face++)
{
HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
}
}
return tHpsi;
}
template<class Type, class DType, class LUType>
Foam::tmp<Foam::Field<Type> >
Foam::LduMatrix<Type, DType, LUType>::H(const tmp<Field<Type> >& tpsi) const
{
tmp<Field<Type> > tHpsi(H(tpsi()));
tpsi.clear();
return tHpsi;
}
template<class Type, class DType, class LUType>
Foam::tmp<Foam::Field<Type> >
Foam::LduMatrix<Type, DType, LUType>::faceH(const Field<Type>& psi) const
{
const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
// Take refereces to addressing
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();
tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
Field<Type> & faceHpsi = tfaceHpsi();
for (register label face=0; face<l.size(); face++)
{
faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
}
return tfaceHpsi;
}
template<class Type, class DType, class LUType>
Foam::tmp<Foam::Field<Type> >
Foam::LduMatrix<Type, DType, LUType>::faceH(const tmp<Field<Type> >& tpsi) const
{
tmp<Field<Type> > tfaceHpsi(faceH(tpsi()));
tpsi.clear();
return tfaceHpsi;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::operator=(const LduMatrix& A)
{
if (this == &A)
{
FatalErrorIn
(
"LduMatrix<Type, DType, LUType>::operator=(const LduMatrix&)"
) << "attempted assignment to self"
<< abort(FatalError);
}
if (A.diagPtr_)
{
diag() = A.diag();
}
if (A.upperPtr_)
{
upper() = A.upper();
}
else if (upperPtr_)
{
delete upperPtr_;
upperPtr_ = NULL;
}
if (A.lowerPtr_)
{
lower() = A.lower();
}
else if (lowerPtr_)
{
delete lowerPtr_;
lowerPtr_ = NULL;
}
if (A.sourcePtr_)
{
source() = A.source();
}
interfacesUpper_ = A.interfacesUpper_;
interfacesLower_ = A.interfacesLower_;
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::negate()
{
if (diagPtr_)
{
diagPtr_->negate();
}
if (upperPtr_)
{
upperPtr_->negate();
}
if (lowerPtr_)
{
lowerPtr_->negate();
}
if (sourcePtr_)
{
sourcePtr_->negate();
}
negate(interfacesUpper_);
negate(interfacesLower_);
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
{
if (A.diagPtr_)
{
diag() += A.diag();
}
if (A.sourcePtr_)
{
source() += A.source();
}
if (symmetric() && A.symmetric())
{
upper() += A.upper();
}
else if (symmetric() && A.asymmetric())
{
if (upperPtr_)
{
lower();
}
else
{
upper();
}
upper() += A.upper();
lower() += A.lower();
}
else if (asymmetric() && A.symmetric())
{
if (A.upperPtr_)
{
lower() += A.upper();
upper() += A.upper();
}
else
{
lower() += A.lower();
upper() += A.lower();
}
}
else if (asymmetric() && A.asymmetric())
{
lower() += A.lower();
upper() += A.upper();
}
else if (diagonal())
{
if (A.upperPtr_)
{
upper() = A.upper();
}
if (A.lowerPtr_)
{
lower() = A.lower();
}
}
else if (A.diagonal())
{
}
else
{
FatalErrorIn
(
"LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)"
) << "Unknown matrix type combination"
<< abort(FatalError);
}
interfacesUpper_ += A.interfacesUpper_;
interfacesLower_ += A.interfacesLower_;
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
{
if (A.diagPtr_)
{
diag() -= A.diag();
}
if (A.sourcePtr_)
{
source() -= A.source();
}
if (symmetric() && A.symmetric())
{
upper() -= A.upper();
}
else if (symmetric() && A.asymmetric())
{
if (upperPtr_)
{
lower();
}
else
{
upper();
}
upper() -= A.upper();
lower() -= A.lower();
}
else if (asymmetric() && A.symmetric())
{
if (A.upperPtr_)
{
lower() -= A.upper();
upper() -= A.upper();
}
else
{
lower() -= A.lower();
upper() -= A.lower();
}
}
else if (asymmetric() && A.asymmetric())
{
lower() -= A.lower();
upper() -= A.upper();
}
else if (diagonal())
{
if (A.upperPtr_)
{
upper() = -A.upper();
}
if (A.lowerPtr_)
{
lower() = -A.lower();
}
}
else if (A.diagonal())
{
}
else
{
FatalErrorIn
(
"LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)"
) << "Unknown matrix type combination"
<< abort(FatalError);
}
interfacesUpper_ -= A.interfacesUpper_;
interfacesLower_ -= A.interfacesLower_;
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::operator*=
(
const scalarField& sf
)
{
if (diagPtr_)
{
*diagPtr_ *= sf;
}
if (sourcePtr_)
{
*sourcePtr_ *= sf;
}
if (upperPtr_)
{
Field<LUType>& upper = *upperPtr_;
const unallocLabelList& l = lduAddr().lowerAddr();
for (register label face=0; face<upper.size(); face++)
{
upper[face] *= sf[l[face]];
}
}
if (lowerPtr_)
{
Field<LUType>& lower = *lowerPtr_;
const unallocLabelList& u = lduAddr().upperAddr();
for (register label face=0; face<lower.size(); face++)
{
lower[face] *= sf[u[face]];
}
}
FatalErrorIn
(
"LduMatrix<Type, DType, LUType>::operator*=(const scalarField& sf)"
) << "Scaling a matrix by scalarField is not currently supported\n"
"because scaling interfacesUpper_ and interfacesLower_ "
"require special transfers"
<< abort(FatalError);
//interfacesUpper_ *= ;
//interfacesLower_ *= sf;
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::operator*=(scalar s)
{
if (diagPtr_)
{
*diagPtr_ *= s;
}
if (sourcePtr_)
{
*sourcePtr_ *= s;
}
if (upperPtr_)
{
*upperPtr_ *= s;
}
if (lowerPtr_)
{
*lowerPtr_ *= s;
}
interfacesUpper_ *= s;
interfacesLower_ *= s;
}
// ************************************************************************* //

116
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixPreconditioner.C
View File

@ -1,116 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::preconditioner>
Foam::LduMatrix<Type, DType, LUType>::preconditioner::New
(
const solver& sol,
const dictionary& preconditionerDict
)
{
word preconditionerName = preconditionerDict.lookup("preconditioner");
if (sol.matrix().symmetric())
{
typename symMatrixConstructorTable::iterator constructorIter =
symMatrixConstructorTablePtr_->find(preconditionerName);
if (constructorIter == symMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::preconditioner::New"
"(const solver&, Istream&)",
preconditionerDict
) << "Unknown symmetric matrix preconditioner "
<< preconditionerName << endl << endl
<< "Valid symmetric matrix preconditioners are :" << endl
<< symMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
(
constructorIter()
(
sol,
preconditionerDict
)
);
}
else if (sol.matrix().asymmetric())
{
typename asymMatrixConstructorTable::iterator constructorIter =
asymMatrixConstructorTablePtr_->find(preconditionerName);
if (constructorIter == asymMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::preconditioner::New"
"(const solver&, Istream&)",
preconditionerDict
) << "Unknown asymmetric matrix preconditioner "
<< preconditionerName << endl << endl
<< "Valid asymmetric matrix preconditioners are :" << endl
<< asymMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
(
constructorIter()
(
sol,
preconditionerDict
)
);
}
else
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::preconditioner::New"
"(const solver&, Istream&)",
preconditionerDict
) << "cannot preconditione incomplete matrix, "
"no diagonal or off-diagonal coefficient"
<< exit(FatalIOError);
return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
(
NULL
);
}
}
// ************************************************************************* //

121
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSmoother.C
View File

@ -1,121 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::smoother>
Foam::LduMatrix<Type, DType, LUType>::smoother::New
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& smootherDict
)
{
word smootherName = smootherDict.lookup("smoother");
if (matrix.symmetric())
{
typename symMatrixConstructorTable::iterator constructorIter =
symMatrixConstructorTablePtr_->find(smootherName);
if (constructorIter == symMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
) << "Unknown symmetric matrix smoother " << smootherName
<< endl << endl
<< "Valid symmetric matrix smoothers are :" << endl
<< symMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
(
constructorIter()
(
fieldName,
matrix
)
);
}
else if (matrix.asymmetric())
{
typename asymMatrixConstructorTable::iterator constructorIter =
asymMatrixConstructorTablePtr_->find(smootherName);
if (constructorIter == asymMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
) << "Unknown asymmetric matrix smoother " << smootherName
<< endl << endl
<< "Valid asymmetric matrix smoothers are :" << endl
<< asymMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
(
constructorIter()
(
fieldName,
matrix
)
);
}
else
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
) << "cannot solve incomplete matrix, no off-diagonal coefficients"
<< exit(FatalIOError);
return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>(NULL);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::smoother::smoother
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
)
:
fieldName_(fieldName),
matrix_(matrix)
{}
// ************************************************************************* //

191
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
View File

@ -1,191 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
#include "DiagonalSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::solver>
Foam::LduMatrix<Type, DType, LUType>::solver::New
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
{
word solverName = solverDict.lookup("solver");
if (matrix.diagonal())
{
return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
(
new DiagonalSolver<Type, DType, LUType>
(
fieldName,
matrix,
solverDict
)
);
}
else if (matrix.symmetric())
{
typename symMatrixConstructorTable::iterator constructorIter =
symMatrixConstructorTablePtr_->find(solverName);
if (constructorIter == symMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::solver::New", solverDict
) << "Unknown symmetric matrix solver " << solverName
<< endl << endl
<< "Valid symmetric matrix solvers are :" << endl
<< symMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
(
constructorIter()
(
fieldName,
matrix,
solverDict
)
);
}
else if (matrix.asymmetric())
{
typename asymMatrixConstructorTable::iterator constructorIter =
asymMatrixConstructorTablePtr_->find(solverName);
if (constructorIter == asymMatrixConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::solver::New", solverDict
) << "Unknown asymmetric matrix solver " << solverName
<< endl << endl
<< "Valid asymmetric matrix solvers are :" << endl
<< asymMatrixConstructorTablePtr_->toc()
<< exit(FatalIOError);
}
return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
(
constructorIter()
(
fieldName,
matrix,
solverDict
)
);
}
else
{
FatalIOErrorIn
(
"LduMatrix<Type, DType, LUType>::solver::New", solverDict
) << "cannot solve incomplete matrix, "
"no diagonal or off-diagonal coefficient"
<< exit(FatalIOError);
return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>(NULL);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::LduMatrix<Type, DType, LUType>::solver::solver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
fieldName_(fieldName),
matrix_(matrix),
controlDict_(solverDict),
maxIter_(1000),
tolerance_(1e-6*pTraits<Type>::one),
relTol_(pTraits<Type>::zero)
{
readControls();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::solver::readControls()
{
readControl(controlDict_, maxIter_, "maxIter");
readControl(controlDict_, tolerance_, "tolerance");
readControl(controlDict_, relTol_, "relTol");
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::solver::read
(
const dictionary& solverDict
)
{
controlDict_ = solverDict;
readControls();
}
template<class Type, class DType, class LUType>
Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
(
const Field<Type>& psi,
const Field<Type>& Apsi,
Field<Type>& tmpField
) const
{
// --- Calculate A dot reference value of psi
matrix_.sumA(tmpField);
cmptMultiply(tmpField, tmpField, gAverage(psi));
return stabilise
(
gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
matrix_.small_
);
// At convergence this simpler method is equivalent to the above
// return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
}
// ************************************************************************* //

119
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C
View File

@ -1,119 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular
(
const Type& wApA
)
{
for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
{
singular_[cmpt] = component(wApA, cmpt) < vsmall_;
}
return singular();
}
template<class Type, class DType, class LUType>
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const
{
for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
{
if (!singular_[cmpt]) return false;
}
return true;
}
template<class Type, class DType, class LUType>
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
(
const Type& Tolerance,
const Type& RelTolerance
)
{
if (debug >= 2)
{
Info<< solverName_
<< ": Iteration " << noIterations_
<< " residual = " << finalResidual_
<< endl;
}
if
(
finalResidual_ < Tolerance
|| (
RelTolerance > small_*pTraits<Type>::one
&& finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
)
)
{
converged_ = true;
}
else
{
converged_ = false;
}
return converged_;
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
(
Ostream& os
) const
{
for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
{
os << solverName_ << ": Solving for "
<< word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
if (singular_[cmpt])
{
os << ": solution singularity" << endl;
}
else
{
os << ", Initial residual = " << component(initialResidual_, cmpt)
<< ", Final residual = " << component(finalResidual_, cmpt)
<< ", No Iterations " << noIterations_
<< endl;
}
}
}
// ************************************************************************* //

195
src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C
View File

@ -1,195 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
#include "LduInterfaceField.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
(
const FieldField<Field, LUType>& interfaceCoeffs,
const Field<Type>& psiif,
Field<Type>& result
) const
{
if
(
Pstream::defaultCommsType == Pstream::blocking
|| Pstream::defaultCommsType == Pstream::nonBlocking
)
{
forAll (interfaces_, interfaceI)
{
if (interfaces_.set(interfaceI))
{
interfaces_[interfaceI].initInterfaceMatrixUpdate
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::defaultCommsType
);
}
}
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();
// Loop over the "global" patches are on the list of interfaces but
// beyond the end of the schedule which only handles "normal" patches
for
(
label interfaceI=patchSchedule.size()/2;
interfaceI<interfaces_.size();
interfaceI++
)
{
if (interfaces_.set(interfaceI))
{
interfaces_[interfaceI].initInterfaceMatrixUpdate
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::blocking
);
}
}
}
else
{
FatalErrorIn("LduMatrix<Type, DType, LUType>::initMatrixInterfaces")
<< "Unsuported communications type "
<< Pstream::commsTypeNames[Pstream::defaultCommsType]
<< exit(FatalError);
}
}
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
(
const FieldField<Field, LUType>& interfaceCoeffs,
const Field<Type>& psiif,
Field<Type>& result
) const
{
if
(
Pstream::defaultCommsType == Pstream::blocking
|| Pstream::defaultCommsType == Pstream::nonBlocking
)
{
// Block until all sends/receives have been finished
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
}
forAll (interfaces_, interfaceI)
{
if (interfaces_.set(interfaceI))
{
interfaces_[interfaceI].updateInterfaceMatrix
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::defaultCommsType
);
}
}
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();
// Loop over all the "normal" interfaces relating to standard patches
forAll (patchSchedule, i)
{
label interfaceI = patchSchedule[i].patch;
if (interfaces_.set(interfaceI))
{
if (patchSchedule[i].init)
{
interfaces_[interfaceI].initInterfaceMatrixUpdate
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::scheduled
);
}
else
{
interfaces_[interfaceI].updateInterfaceMatrix
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::scheduled
);
}
}
}
// Loop over the "global" patches are on the list of interfaces but
// beyond the end of the schedule which only handles "normal" patches
for
(
label interfaceI=patchSchedule.size()/2;
interfaceI<interfaces_.size();
interfaceI++
)
{
if (interfaces_.set(interfaceI))
{
interfaces_[interfaceI].updateInterfaceMatrix
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::blocking
);
}
}
}
else
{
FatalErrorIn("LduMatrix<Type, DType, LUType>::updateMatrixInterfaces")
<< "Unsuported communications type "
<< Pstream::commsTypeNames[Pstream::defaultCommsType]
<< exit(FatalError);
}
}
// ************************************************************************* //

11
src/OpenFOAM/matrices/LduMatrix/LduMatrix/lduMatrices.C
View File

@ -1,11 +0,0 @@
#include "LduMatrix.H"
#include "fieldTypes.H"
namespace Foam
{
makeLduMatrix(scalar, scalar, scalar);
makeLduMatrix(vector, scalar, scalar);
makeLduMatrix(sphericalTensor, scalar, scalar);
makeLduMatrix(symmTensor, scalar, scalar);
makeLduMatrix(tensor, scalar, scalar);
};

180
src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.C
View File

@ -1,180 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "TDILUPreconditioner.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::TDILUPreconditioner<Type, DType, LUType>::TDILUPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary&
)
:
LduMatrix<Type, DType, LUType>::preconditioner(sol),
rD_(sol.matrix().diag())
{
calcInvD(rD_, sol.matrix());
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::TDILUPreconditioner<Type, DType, LUType>::calcInvD
(
Field<DType>& rD,
const LduMatrix<Type, DType, LUType>& matrix
)
{
DType* __restrict__ rDPtr = rD.begin();
const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
const LUType* const __restrict__ upperPtr = matrix.upper().begin();
const LUType* const __restrict__ lowerPtr = matrix.lower().begin();
register label nFaces = matrix.upper().size();
for (register label face=0; face<nFaces; face++)
{
rDPtr[uPtr[face]] -=
dot(dot(upperPtr[face], lowerPtr[face]), inv(rDPtr[lPtr[face]]));
}
// Calculate the reciprocal of the preconditioned diagonal
register label nCells = rD.size();
for (register label cell=0; cell<nCells; cell++)
{
rDPtr[cell] = inv(rDPtr[cell]);
}
}
template<class Type, class DType, class LUType>
void Foam::TDILUPreconditioner<Type, DType, LUType>::precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const
{
Type* __restrict__ wAPtr = wA.begin();
const Type* __restrict__ rAPtr = rA.begin();
const DType* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
this->solver_.matrix().lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr =
this->solver_.matrix().lduAddr().lowerAddr().begin();
const label* const __restrict__ losortPtr =
this->solver_.matrix().lduAddr().losortAddr().begin();
const LUType* const __restrict__ upperPtr =
this->solver_.matrix().upper().begin();
const LUType* const __restrict__ lowerPtr =
this->solver_.matrix().lower().begin();
register label nCells = wA.size();
register label nFaces = this->solver_.matrix().upper().size();
register label nFacesM1 = nFaces - 1;
for (register label cell=0; cell<nCells; cell++)
{
wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
}
register label sface;
for (register label face=0; face<nFaces; face++)
{
sface = losortPtr[face];
wAPtr[uPtr[sface]] -=
dot(rDPtr[uPtr[sface]], dot(lowerPtr[sface], wAPtr[lPtr[sface]]));
}
for (register label face=nFacesM1; face>=0; face--)
{
wAPtr[lPtr[face]] -=
dot(rDPtr[lPtr[face]], dot(upperPtr[face], wAPtr[uPtr[face]]));
}
}
template<class Type, class DType, class LUType>
void Foam::TDILUPreconditioner<Type, DType, LUType>::preconditionT
(
Field<Type>& wT,
const Field<Type>& rT
) const
{
Type* __restrict__ wTPtr = wT.begin();
const Type* __restrict__ rTPtr = rT.begin();
const DType* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
this->solver_.matrix().lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr =
this->solver_.matrix().lduAddr().lowerAddr().begin();
const label* const __restrict__ losortPtr =
this->solver_.matrix().lduAddr().losortAddr().begin();
const LUType* const __restrict__ upperPtr =
this->solver_.matrix().upper().begin();
const LUType* const __restrict__ lowerPtr =
this->solver_.matrix().lower().begin();
register label nCells = wT.size();
register label nFaces = this->solver_.matrix().upper().size();
register label nFacesM1 = nFaces - 1;
for (register label cell=0; cell<nCells; cell++)
{
wTPtr[cell] = dot(rDPtr[cell], rTPtr[cell]);
}
for (register label face=0; face<nFaces; face++)
{
wTPtr[uPtr[face]] -=
dot(rDPtr[uPtr[face]], dot(upperPtr[face], wTPtr[lPtr[face]]));
}
register label sface;
for (register label face=nFacesM1; face>=0; face--)
{
sface = losortPtr[face];
wTPtr[lPtr[sface]] -=
dot(rDPtr[lPtr[sface]], dot(lowerPtr[sface], wTPtr[uPtr[sface]]));
}
}
// ************************************************************************* //

127
src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.H
View File

@ -1,127 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::TDILUPreconditioner
Description
Simplified diagonal-based incomplete LU preconditioner for asymmetric
matrices.
The inverse (reciprocal for scalar) of the preconditioned diagonal is
calculated and stored.
SourceFiles
TDILUPreconditioner.C
\*---------------------------------------------------------------------------*/
#ifndef TDILUPreconditioner_H
#define TDILUPreconditioner_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TDILUPreconditioner Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class TDILUPreconditioner
:
public LduMatrix<Type, DType, LUType>::preconditioner
{
// Private data
//- The inverse (reciprocal for scalar) preconditioned diagonal
Field<DType> rD_;
public:
//- Runtime type information
TypeName("DILU");
// Constructors
//- Construct from matrix components and preconditioner data dictionary
TDILUPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary& preconditionerDict
);
// Destructor
virtual ~TDILUPreconditioner()
{}
// Member Functions
//- Calculate the reciprocal of the preconditioned diagonal
static void calcInvD
(
Field<DType>& rD,
const LduMatrix<Type, DType, LUType>& matrix
);
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const;
//- Return wT the transpose-matrix preconditioned form of
// residual rT.
virtual void preconditionT
(
Field<Type>& wT,
const Field<Type>& rT
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "TDILUPreconditioner.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

81
src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.C
View File

@ -1,81 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "DiagonalPreconditioner.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::DiagonalPreconditioner<Type, DType, LUType>::DiagonalPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary&
)
:
LduMatrix<Type, DType, LUType>::preconditioner(sol),
rD(sol.matrix().diag().size())
{
DType* __restrict__ rDPtr = rD.begin();
const DType* __restrict__ DPtr = this->solver_.matrix().diag().begin();
register label nCells = rD.size();
// Generate inverse (reciprocal for scalar) diagonal
for (register label cell=0; cell<nCells; cell++)
{
rDPtr[cell] = inv(DPtr[cell]);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::DiagonalPreconditioner<Type, DType, LUType>::read(const dictionary&)
{}
template<class Type, class DType, class LUType>
void Foam::DiagonalPreconditioner<Type, DType, LUType>::precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const
{
Type* __restrict__ wAPtr = wA.begin();
const Type* __restrict__ rAPtr = rA.begin();
const DType* __restrict__ rDPtr = rD.begin();
register label nCells = wA.size();
for (register label cell=0; cell<nCells; cell++)
{
wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
}
}
// ************************************************************************* //

135
src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.H
View File

@ -1,135 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::DiagonalPreconditioner
Description
Diagonal preconditioner for both symmetric and asymmetric matrices.
The inverse (reciprocal for scalar) of the diagonal is calculated and
stored.
SourceFiles
DiagonalPreconditioner.C
\*---------------------------------------------------------------------------*/
#ifndef DiagonalPreconditioner_H
#define DiagonalPreconditioner_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class DiagonalPreconditioner Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class DiagonalPreconditioner
:
public LduMatrix<Type, DType, LUType>::preconditioner
{
// Private data
//- The inverse (reciprocal for scalar) diagonal
Field<DType> rD;
// Private Member Functions
//- Disallow default bitwise copy construct
DiagonalPreconditioner(const DiagonalPreconditioner&);
//- Disallow default bitwise assignment
void operator=(const DiagonalPreconditioner&);
public:
//- Runtime type information
TypeName("diagonal");
// Constructors
//- Construct from matrix components and preconditioner data dictionary
DiagonalPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary& preconditionerDict
);
// Destructor
virtual ~DiagonalPreconditioner()
{}
// Member Functions
//- Read and reset the preconditioner parameters from the given
// dictionary
virtual void read(const dictionary& preconditionerDict);
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const;
//- Return wT the transpose-matrix preconditioned form of
// residual rT.
virtual void preconditionT
(
Field<Type>& wT,
const Field<Type>& rT
) const
{
return(precondition(wT, rT));
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "DiagonalPreconditioner.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

60
src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.C
View File

@ -1,60 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "NoPreconditioner.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::NoPreconditioner<Type, DType, LUType>::NoPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary&
)
:
LduMatrix<Type, DType, LUType>::preconditioner(sol)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::NoPreconditioner<Type, DType, LUType>::read(const dictionary&)
{}
template<class Type, class DType, class LUType>
void Foam::NoPreconditioner<Type, DType, LUType>::precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const
{
wA = rA;
}
// ************************************************************************* //

126
src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.H
View File

@ -1,126 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::NoPreconditioner
Description
Null preconditioner for both symmetric and asymmetric matrices.
SourceFiles
NoPreconditioner.C
\*---------------------------------------------------------------------------*/
#ifndef NoPreconditioner_H
#define NoPreconditioner_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class NoPreconditioner Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class NoPreconditioner
:
public LduMatrix<Type, DType, LUType>::preconditioner
{
// Private Member Functions
//- Disallow default bitwise copy construct
NoPreconditioner(const NoPreconditioner&);
//- Disallow default bitwise assignment
void operator=(const NoPreconditioner&);
public:
//- Runtime type information
TypeName("none");
// Constructors
//- Construct from matrix components and preconditioner data dictionary
NoPreconditioner
(
const typename LduMatrix<Type, DType, LUType>::solver& sol,
const dictionary& preconditionerDict
);
// Destructor
virtual ~NoPreconditioner()
{}
// Member Functions
//- Read and reset the preconditioner parameters from the given
// dictionary
virtual void read(const dictionary& preconditionerDict);
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
Field<Type>& wA,
const Field<Type>& rA
) const;
//- Return wT the transpose-matrix preconditioned form of
// residual rT.
virtual void preconditionT
(
Field<Type>& wT,
const Field<Type>& rT
) const
{
return(precondition(wT, rT));
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "NoPreconditioner.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

26
src/OpenFOAM/matrices/LduMatrix/Preconditioners/lduPreconditioners.C
View File

@ -1,26 +0,0 @@
#include "NoPreconditioner.H"
#include "DiagonalPreconditioner.H"
#include "TDILUPreconditioner.H"
#include "fieldTypes.H"
#define makeLduPreconditioners(Type, DType, LUType) \
\
makeLduPreconditioner(NoPreconditioner, Type, DType, LUType); \
makeLduSymPreconditioner(NoPreconditioner, Type, DType, LUType); \
makeLduAsymPreconditioner(NoPreconditioner, Type, DType, LUType); \
\
makeLduPreconditioner(DiagonalPreconditioner, Type, DType, LUType); \
makeLduSymPreconditioner(DiagonalPreconditioner, Type, DType, LUType); \
makeLduAsymPreconditioner(DiagonalPreconditioner, Type, DType, LUType); \
\
makeLduPreconditioner(TDILUPreconditioner, Type, DType, LUType); \
makeLduAsymPreconditioner(TDILUPreconditioner, Type, DType, LUType);
namespace Foam
{
makeLduPreconditioners(scalar, scalar, scalar);
makeLduPreconditioners(vector, scalar, scalar);
makeLduPreconditioners(sphericalTensor, scalar, scalar);
makeLduPreconditioners(symmTensor, scalar, scalar);
makeLduPreconditioners(tensor, scalar, scalar);
};

168
src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.C
View File

@ -1,168 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "TGaussSeidelSmoother.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::TGaussSeidelSmoother<Type, DType, LUType>::TGaussSeidelSmoother
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
)
:
LduMatrix<Type, DType, LUType>::smoother
(
fieldName,
matrix
),
rD_(matrix.diag().size())
{
register const label nCells = matrix.diag().size();
register const DType* const __restrict__ diagPtr = matrix.diag().begin();
register DType* __restrict__ rDPtr = rD_.begin();
for (register label cellI=0; cellI<nCells; cellI++)
{
rDPtr[cellI] = inv(diagPtr[cellI]);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
(
const word& fieldName_,
Field<Type>& psi,
const LduMatrix<Type, DType, LUType>& matrix_,
const Field<DType>& rD_,
const label nSweeps
)
{
register Type* __restrict__ psiPtr = psi.begin();
register const label nCells = psi.size();
Field<Type> bPrime(nCells);
register Type* __restrict__ bPrimePtr = bPrime.begin();
register const DType* const __restrict__ rDPtr = rD_.begin();
register const LUType* const __restrict__ upperPtr =
matrix_.upper().begin();
register const LUType* 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 to add the contibution to the r.h.s.
FieldField<Field, LUType> mBouCoeffs(matrix_.interfacesUpper().size());
forAll(mBouCoeffs, patchi)
{
if (matrix_.interfaces().set(patchi))
{
mBouCoeffs.set(patchi, -matrix_.interfacesUpper()[patchi]);
}
}
for (label sweep=0; sweep<nSweeps; sweep++)
{
bPrime = matrix_.source();
matrix_.initMatrixInterfaces
(
mBouCoeffs,
psi,
bPrime
);
matrix_.updateMatrixInterfaces
(
mBouCoeffs,
psi,
bPrime
);
Type curPsi;
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
curPsi = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
curPsi -= dot(upperPtr[curFace], psiPtr[uPtr[curFace]]);
}
// Finish current psi
curPsi = dot(rDPtr[cellI], curPsi);
// Distribute the neighbour side using current psi
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= dot(lowerPtr[curFace], curPsi);
}
psiPtr[cellI] = curPsi;
}
}
}
template<class Type, class DType, class LUType>
void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
(
Field<Type>& psi,
const label nSweeps
) const
{
smooth(this->fieldName_, psi, this->matrix_, rD_, nSweeps);
}
// ************************************************************************* //

113
src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.H
View File

@ -1,113 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::TGaussSeidelSmoother
Description
Foam::TGaussSeidelSmoother
SourceFiles
TGaussSeidelSmoother.C
\*---------------------------------------------------------------------------*/
#ifndef TGaussSeidelSmoother_H
#define TGaussSeidelSmoother_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TGaussSeidelSmoother Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class TGaussSeidelSmoother
:
public LduMatrix<Type, DType, LUType>::smoother
{
// Private data
//- The inverse (reciprocal for scalars) diagonal
Field<DType> rD_;
public:
//- Runtime type information
TypeName("GaussSeidel");
// Constructors
//- Construct from components
TGaussSeidelSmoother
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix
);
// Member Functions
//- Smooth for the given number of sweeps
static void smooth
(
const word& fieldName,
Field<Type>& psi,
const LduMatrix<Type, DType, LUType>& matrix,
const Field<DType>& rD,
const label nSweeps
);
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
Field<Type>& psi,
const label nSweeps
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "TGaussSeidelSmoother.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

17
src/OpenFOAM/matrices/LduMatrix/Smoothers/lduSmoothers.C
View File

@ -1,17 +0,0 @@
#include "TGaussSeidelSmoother.H"
#include "fieldTypes.H"
#define makeLduSmoothers(Type, DType, LUType) \
\
makeLduSmoother(TGaussSeidelSmoother, Type, DType, LUType); \
makeLduSymSmoother(TGaussSeidelSmoother, Type, DType, LUType); \
makeLduAsymSmoother(TGaussSeidelSmoother, Type, DType, LUType);
namespace Foam
{
makeLduSmoothers(scalar, scalar, scalar);
makeLduSmoothers(vector, scalar, scalar);
makeLduSmoothers(sphericalTensor, scalar, scalar);
makeLduSmoothers(symmTensor, scalar, scalar);
makeLduSmoothers(tensor, scalar, scalar);
};

80
src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
View File

@ -1,80 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "DiagonalSolver.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::DiagonalSolver<Type, DType, LUType>::DiagonalSolver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
)
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::DiagonalSolver<Type, DType, LUType>::read
(
const dictionary&
)
{}
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
Foam::DiagonalSolver<Type, DType, LUType>::solve
(
Field<Type>& psi
) const
{
psi = this->matrix_.source()/this->matrix_.diag();
return typename LduMatrix<Type, DType, LUType>::solverPerformance
(
typeName,
this->fieldName_,
pTraits<Type>::zero,
pTraits<Type>::zero,
0,
true,
false
);
}
// ************************************************************************* //

108
src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
View File

@ -1,108 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::DiagonalSolver
Description
Foam::DiagonalSolver
SourceFiles
DiagonalSolver.C
\*---------------------------------------------------------------------------*/
#ifndef DiagonalSolver_H
#define DiagonalSolver_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class DiagonalSolver Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class DiagonalSolver
:
public LduMatrix<Type, DType, LUType>::solver
{
// Private Member Functions
//- Disallow default bitwise copy construct
DiagonalSolver(const DiagonalSolver&);
//- Disallow default bitwise assignment
void operator=(const DiagonalSolver&);
public:
//- Runtime type information
TypeName("diagonal");
// Constructors
//- Construct from matrix
DiagonalSolver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Member Functions
//- Read and reset the solver parameters from the given dictionary
void read(const dictionary& solverDict);
//- Solve the matrix with this solver
typename LduMatrix<Type, DType, LUType>::solverPerformance solve
(
Field<Type>& psi
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "DiagonalSolver.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

195
src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.C
View File

@ -1,195 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "TPBiCG.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::TPBiCG<Type, DType, LUType>::TPBiCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
Foam::TPBiCG<Type, DType, LUType>::solve(Field<Type>& psi) const
{
word preconditionerName(this->controlDict_.lookup("preconditioner"));
// --- Setup class containing solver performance data
typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
(
preconditionerName + typeName,
this->fieldName_
);
register label nCells = psi.size();
Type* __restrict__ psiPtr = psi.begin();
Field<Type> pA(nCells);
Type* __restrict__ pAPtr = pA.begin();
Field<Type> pT(nCells, pTraits<Type>::zero);
Type* __restrict__ pTPtr = pT.begin();
Field<Type> wA(nCells);
Type* __restrict__ wAPtr = wA.begin();
Field<Type> wT(nCells);
Type* __restrict__ wTPtr = wT.begin();
Type wArT = this->matrix_.great_*pTraits<Type>::one;
Type wArTold = wArT;
// --- Calculate A.psi and T.psi
this->matrix_.Amul(wA, psi);
this->matrix_.Tmul(wT, psi);
// --- Calculate initial residual and transpose residual fields
Field<Type> rA(this->matrix_.source() - wA);
Field<Type> rT(this->matrix_.source() - wT);
Type* __restrict__ rAPtr = rA.begin();
Type* __restrict__ rTPtr = rT.begin();
// --- Calculate normalisation factor
Type normFactor = this->normFactor(psi, wA, pA);
if (LduMatrix<Type, DType, LUType>::debug >= 2)
{
Info<< " Normalisation factor = " << normFactor << endl;
}
// --- Calculate normalised residual norm
solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
solverPerf.finalResidual() = solverPerf.initialResidual();
// --- Check convergence, solve if not converged
if (!solverPerf.converged(this->tolerance_, this->relTol_))
{
// --- Select and construct the preconditioner
autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
(
*this,
this->controlDict_
);
// --- Solver iteration
do
{
// --- Store previous wArT
wArTold = wArT;
// --- Precondition residuals
preconPtr->precondition(wA, rA);
preconPtr->preconditionT(wT, rT);
// --- Update search directions:
wArT = gSumCmptProd(wA, rT);
if (solverPerf.nIterations() == 0)
{
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell];
pTPtr[cell] = wTPtr[cell];
}
}
else
{
Type beta = cmptDivide
(
wArT,
stabilise(wArTold, this->matrix_.vsmall_)
);
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
pTPtr[cell] = wTPtr[cell] + cmptMultiply(beta, pTPtr[cell]);
}
}
// --- Update preconditioned residuals
this->matrix_.Amul(wA, pA);
this->matrix_.Tmul(wT, pT);
Type wApT = gSumCmptProd(wA, pT);
// --- Test for singularity
if (solverPerf.singular(cmptDivide(cmptMag(wApT), normFactor)))
{
break;
}
// --- Update solution and residual:
Type alpha = cmptDivide
(
wArT,
stabilise(wApT, this->matrix_.vsmall_)
);
for (register label cell=0; cell<nCells; cell++)
{
psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
rTPtr[cell] -= cmptMultiply(alpha, wTPtr[cell]);
}
solverPerf.finalResidual() =
cmptDivide(gSumCmptMag(rA), normFactor);
} while
(
solverPerf.nIterations()++ < this->maxIter_
&& !(solverPerf.converged(this->tolerance_, this->relTol_))
);
}
return solverPerf;
}
// ************************************************************************* //

112
src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.H
View File

@ -1,112 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::TPBiCG
Description
Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
using a run-time selectable preconditiioner.
SourceFiles
TPBiCG.C
\*---------------------------------------------------------------------------*/
#ifndef TPBiCG_H
#define TPBiCG_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TPBiCG Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class TPBiCG
:
public LduMatrix<Type, DType, LUType>::solver
{
// Private Member Functions
//- Disallow default bitwise copy construct
TPBiCG(const TPBiCG&);
//- Disallow default bitwise assignment
void operator=(const TPBiCG&);
public:
//- Runtime type information
TypeName("PBiCG");
// Constructors
//- Construct from matrix components and solver data dictionary
TPBiCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Destructor
virtual ~TPBiCG()
{}
// Member Functions
//- Solve the matrix with this solver
virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
(
Field<Type>& psi
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "TPBiCG.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

191
src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.C
View File

@ -1,191 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "TPBiCG.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::TPBiCG<Type, DType, LUType>::TPBiCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
Foam::TPBiCG<Type, DType, LUType>::solve
(
Field<Type>& psi
) const
{
word preconditionerName(this->controlDict_.lookup("preconditioner"));
// --- Setup class containing solver performance data
typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
(
preconditionerName + typeName,
this->fieldName_
);
register label nCells = psi.size();
Type* __restrict__ psiPtr = psi.begin();
Field<Type> pA(nCells);
Type* __restrict__ pAPtr = pA.begin();
Field<Type> pT(nCells, pTraits<Type>::zero);
Type* __restrict__ pTPtr = pT.begin();
Field<Type> wA(nCells);
Type* __restrict__ wAPtr = wA.begin();
Field<Type> wT(nCells);
Type* __restrict__ wTPtr = wT.begin();
scalar wArT = 1e15; //this->matrix_.great_;
scalar wArTold = wArT;
// --- Calculate A.psi and T.psi
this->matrix_.Amul(wA, psi);
this->matrix_.Tmul(wT, psi);
// --- Calculate initial residual and transpose residual fields
Field<Type> rA(this->matrix_.source() - wA);
Field<Type> rT(this->matrix_.source() - wT);
Type* __restrict__ rAPtr = rA.begin();
Type* __restrict__ rTPtr = rT.begin();
// --- Calculate normalisation factor
Type normFactor = this->normFactor(psi, wA, pA);
if (LduMatrix<Type, DType, LUType>::debug >= 2)
{
Info<< " Normalisation factor = " << normFactor << endl;
}
// --- Calculate normalised residual norm
solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
solverPerf.finalResidual() = solverPerf.initialResidual();
// --- Check convergence, solve if not converged
if (!solverPerf.converged(this->tolerance_, this->relTol_))
{
// --- Select and construct the preconditioner
autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
(
*this,
this->controlDict_
);
// --- Solver iteration
do
{
// --- Store previous wArT
wArTold = wArT;
// --- Precondition residuals
preconPtr->precondition(wA, rA);
preconPtr->preconditionT(wT, rT);
// --- Update search directions:
//wArT = gSumProd(wA, rT);
wArT = sum(wA & rT);
if (solverPerf.nIterations() == 0)
{
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell];
pTPtr[cell] = wTPtr[cell];
}
}
else
{
scalar beta = cmptDivide(wArT, wArTold);
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell] + (beta* pAPtr[cell]);
pTPtr[cell] = wTPtr[cell] + (beta* pTPtr[cell]);
}
}
// --- Update preconditioned residuals
this->matrix_.Amul(wA, pA);
this->matrix_.Tmul(wT, pT);
scalar wApT = sum(wA & pT);
// --- Test for singularity
//if
//(
// solverPerf.checkSingularity(ratio(mag(wApT)normFactor))
//) break;
// --- Update solution and residual:
scalar alpha = cmptDivide(wArT, wApT);
for (register label cell=0; cell<nCells; cell++)
{
psiPtr[cell] += (alpha* pAPtr[cell]);
rAPtr[cell] -= (alpha* wAPtr[cell]);
rTPtr[cell] -= (alpha* wTPtr[cell]);
}
solverPerf.finalResidual() =
cmptDivide(gSumCmptMag(rA), normFactor);
} while
(
solverPerf.nIterations()++ < this->maxIter_
&& !(solverPerf.converged(this->tolerance_, this->relTol_))
);
}
return solverPerf;
}
// ************************************************************************* //

112
src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.H
View File

@ -1,112 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::TPBiCG
Description
Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
using a run-time selectable preconditiioner.
SourceFiles
TPBiCG.C
\*---------------------------------------------------------------------------*/
#ifndef TPBiCG_H
#define TPBiCG_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TPBiCG Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class TPBiCG
:
public LduMatrix<Type, DType, LUType>::solver
{
// Private Member Functions
//- Disallow default bitwise copy construct
TPBiCG(const TPBiCG&);
//- Disallow default bitwise assignment
void operator=(const TPBiCG&);
public:
//- Runtime type information
TypeName("PBiCG");
// Constructors
//- Construct from matrix components and solver data dictionary
TPBiCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Destructor
virtual ~TPBiCG()
{}
// Member Functions
//- Solve the matrix with this solver
virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
(
Field<Type>& psi
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "PBiCGScalarAlpha.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

181
src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.C
View File

@ -1,181 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "TPCG.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::TPCG<Type, DType, LUType>::TPCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
Foam::TPCG<Type, DType, LUType>::solve(Field<Type>& psi) const
{
word preconditionerName(this->controlDict_.lookup("preconditioner"));
// --- Setup class containing solver performance data
typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
(
preconditionerName + typeName,
this->fieldName_
);
register label nCells = psi.size();
Type* __restrict__ psiPtr = psi.begin();
Field<Type> pA(nCells);
Type* __restrict__ pAPtr = pA.begin();
Field<Type> wA(nCells);
Type* __restrict__ wAPtr = wA.begin();
Type wArA = this->matrix_.great_*pTraits<Type>::one;
Type wArAold = wArA;
// --- Calculate A.psi
this->matrix_.Amul(wA, psi);
// --- Calculate initial residual field
Field<Type> rA(this->matrix_.source() - wA);
Type* __restrict__ rAPtr = rA.begin();
// --- Calculate normalisation factor
Type normFactor = this->normFactor(psi, wA, pA);
if (LduMatrix<Type, DType, LUType>::debug >= 2)
{
Info<< " Normalisation factor = " << normFactor << endl;
}
// --- Calculate normalised residual norm
solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
solverPerf.finalResidual() = solverPerf.initialResidual();
// --- Check convergence, solve if not converged
if (!solverPerf.converged(this->tolerance_, this->relTol_))
{
// --- Select and construct the preconditioner
autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
(
*this,
this->controlDict_
);
// --- Solver iteration
do
{
// --- Store previous wArA
wArAold = wArA;
// --- Precondition residual
preconPtr->precondition(wA, rA);
// --- Update search directions:
wArA = gSumCmptProd(wA, rA);
if (solverPerf.nIterations() == 0)
{
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell];
}
}
else
{
Type beta = cmptDivide
(
wArA,
stabilise(wArAold, this->matrix_.vsmall_)
);
for (register label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
}
}
// --- Update preconditioned residual
this->matrix_.Amul(wA, pA);
Type wApA = gSumCmptProd(wA, pA);
// --- Test for singularity
if (solverPerf.singular(cmptDivide(cmptMag(wApA), normFactor)))
{
break;
}
// --- Update solution and residual:
Type alpha = cmptDivide
(
wArA,
stabilise(wApA, this->matrix_.vsmall_)
);
for (register label cell=0; cell<nCells; cell++)
{
psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
}
solverPerf.finalResidual() =
cmptDivide(gSumCmptMag(rA), normFactor);
} while
(
solverPerf.nIterations()++ < this->maxIter_
&& !(solverPerf.converged(this->tolerance_, this->relTol_))
);
}
return solverPerf;
}
// ************************************************************************* //

112
src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.H
View File

@ -1,112 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::TPCG
Description
Preconditioned conjugate gradient solver for symmetric lduMatrices
using a run-time selectable preconditiioner.
SourceFiles
TPCG.C
\*---------------------------------------------------------------------------*/
#ifndef TPCG_H
#define TPCG_H
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class TPCG Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class TPCG
:
public LduMatrix<Type, DType, LUType>::solver
{
// Private Member Functions
//- Disallow default bitwise copy construct
TPCG(const TPCG&);
//- Disallow default bitwise assignment
void operator=(const TPCG&);
public:
//- Runtime type information
TypeName("PCG");
// Constructors
//- Construct from matrix components and solver data dictionary
TPCG
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Destructor
virtual ~TPCG()
{}
// Member Functions
//- Solve the matrix with this solver
virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
(
Field<Type>& psi
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "TPCG.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

154
src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
View File

@ -1,154 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "SmoothSolver.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::SmoothSolver<Type, DType, LUType>::SmoothSolver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
),
nSweeps_(1)
{
readControls();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::SmoothSolver<Type, DType, LUType>::readControls()
{
LduMatrix<Type, DType, LUType>::solver::readControls();
readControl(this->controlDict_, nSweeps_, "nSweeps");
}
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const
{
// --- Setup class containing solver performance data
typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
(
typeName,
this->fieldName_
);
// If the nSweeps_ is negative do a fixed number of sweeps
if (nSweeps_ < 0)
{
autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
(
this->fieldName_,
this->matrix_,
this->controlDict_
);
smootherPtr->smooth(psi, -nSweeps_);
solverPerf.nIterations() -= nSweeps_;
}
else
{
Type normFactor = pTraits<Type>::zero;
{
Field<Type> Apsi(psi.size());
Field<Type> temp(psi.size());
// Calculate A.psi
this->matrix_.Amul(Apsi, psi);
// Calculate normalisation factor
normFactor = this->normFactor(psi, Apsi, temp);
// Calculate residual magnitude
solverPerf.initialResidual() = cmptDivide
(
gSumCmptMag(this->matrix_.source() - Apsi),
normFactor
);
solverPerf.finalResidual() = solverPerf.initialResidual();
}
if (LduMatrix<Type, DType, LUType>::debug >= 2)
{
Info<< " Normalisation factor = " << normFactor << endl;
}
// Check convergence, solve if not converged
if (!solverPerf.converged(this->tolerance_, this->relTol_))
{
autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
(
this->fieldName_,
this->matrix_,
this->controlDict_
);
// Smoothing loop
do
{
smootherPtr->smooth
(
psi,
nSweeps_
);
// Calculate the residual to check convergence
solverPerf.finalResidual() = cmptDivide
(
gSumCmptMag(this->matrix_.residual(psi)),
normFactor
);
} while
(
(solverPerf.nIterations() += nSweeps_) < this->maxIter_
&& !(solverPerf.converged(this->tolerance_, this->relTol_))
);
}
}
return solverPerf;
}
// ************************************************************************* //

111
src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
View File

@ -1,111 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::SmoothSolver
Description
Iterative solver for symmetric and assymetric matrices which uses a
run-time selected smoother e.g. GaussSeidel to converge the solution to
the required tolerance. To improve efficiency, the residual is evaluated
after every nSweeps smoothing iterations.
SourceFiles
SmoothSolver.C
\*---------------------------------------------------------------------------*/
#ifndef SmoothSolver_H
#define SmoothSolver_H
#include "lduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class SmoothSolver Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class DType, class LUType>
class SmoothSolver
:
public LduMatrix<Type, DType, LUType>::solver
{
protected:
// Protected data
//- Number of sweeps before the evaluation of residual
label nSweeps_;
//- Read the control parameters from the controlDict_
virtual void readControls();
public:
//- Runtime type information
TypeName("SmoothSolver");
// Constructors
//- Construct from matrix components and solver data dictionary
SmoothSolver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
);
// Member Functions
//- Solve the matrix with this solver
virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
(
Field<Type>& psi
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "SmoothSolver.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

29
src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C
View File

@ -1,29 +0,0 @@
#include "TPCG.H"
#include "TPBiCG.H"
#include "SmoothSolver.H"
#include "fieldTypes.H"
#define makeLduSolvers(Type, DType, LUType) \
\
makeLduSolver(DiagonalSolver, Type, DType, LUType); \
makeLduSymSolver(DiagonalSolver, Type, DType, LUType); \
makeLduAsymSolver(DiagonalSolver, Type, DType, LUType); \
\
makeLduSolver(TPCG, Type, DType, LUType); \
makeLduSymSolver(TPCG, Type, DType, LUType); \
\
makeLduSolver(TPBiCG, Type, DType, LUType); \
makeLduAsymSolver(TPBiCG, Type, DType, LUType); \
\
makeLduSolver(SmoothSolver, Type, DType, LUType); \
makeLduSymSolver(SmoothSolver, Type, DType, LUType); \
makeLduAsymSolver(SmoothSolver, Type, DType, LUType);
namespace Foam
{
makeLduSolvers(scalar, scalar, scalar);
makeLduSolvers(vector, scalar, scalar);
makeLduSolvers(sphericalTensor, scalar, scalar);
makeLduSolvers(symmTensor, scalar, scalar);
makeLduSolvers(tensor, scalar, scalar);
};