Commit 973e2d4e authored by Mark Olesen's avatar Mark Olesen

COMP: remnant foam-extend code in faMatrix::setValues (fixes #1834)

BUG: faMatrix::residual changes source vector (fixes #1835)

ENH: improve code alignment between faMatrix and fvMatrix

- support setValues() with a single value
parent bf3b4fab
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -30,6 +30,8 @@ License
#include "edgeFields.H"
#include "calculatedFaPatchFields.H"
#include "zeroGradientFaPatchFields.H"
#include "UIndirectList.H"
#include "UniformList.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
......@@ -220,10 +222,10 @@ Foam::faMatrix<Type>::faMatrix
}
// Update the boundary coefficients of psi without changing its event No.
GeometricField<Type, faPatchField, areaMesh>& psiRef =
auto& psiRef =
const_cast<GeometricField<Type, faPatchField, areaMesh>&>(psi_);
label currentStatePsi = psiRef.eventNo();
const label currentStatePsi = psiRef.eventNo();
psiRef.boundaryFieldRef().updateCoeffs();
psiRef.eventNo() = currentStatePsi;
}
......@@ -298,7 +300,6 @@ Foam::faMatrix<Type>::faMatrix
)
);
}
}
......@@ -319,8 +320,7 @@ template<class Type>
Foam::faMatrix<Type>::~faMatrix()
{
DebugInFunction
<< "destroying faMatrix<Type> for field " << psi_.name()
<< endl;
<< "Destroying faMatrix<Type> for field " << psi_.name() << endl;
deleteDemandDrivenData(faceFluxCorrectionPtr_);
}
......@@ -329,19 +329,22 @@ Foam::faMatrix<Type>::~faMatrix()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::faMatrix<Type>::setValues
template<template<class> class ListType>
void Foam::faMatrix<Type>::setValuesFromList
(
const labelUList& faceLabels,
const UList<Type>& values
const ListType<Type>& values
)
{
const faMesh& mesh = psi_.mesh();
#if 0 /* Specific to foam-extend */
// Record face labels of eliminated equations
for (const label i : faceLabels)
{
this->eliminatedEqns().insert(i);
}
#endif
const faMesh& mesh = psi_.mesh();
const labelListList& edges = mesh.patch().faceEdges();
const labelUList& own = mesh.owner();
......@@ -352,35 +355,31 @@ void Foam::faMatrix<Type>::setValues
const_cast
<
GeometricField<Type, faPatchField, areaMesh>&
>(psi_).internalField();
>(psi_).primitiveFieldRef();
forAll(faceLabels, i)
{
label facei = faceLabels[i];
const label facei = faceLabels[i];
const Type& value = values[i];
psi[facei] = values[i];
source_[facei] = values[i]*Diag[facei];
psi[facei] = value;
source_[facei] = value*Diag[facei];
if (symmetric() || asymmetric())
{
const labelList& c= edges[facei];
forAll(c, j)
for (const label edgei : edges[facei])
{
label edgei = c[j];
if (mesh.isInternalEdge(edgei))
{
if (symmetric())
{
if (facei == own[edgei])
{
source_[nei[edgei]] -= upper()[edgei]*values[i];
source_[nei[edgei]] -= upper()[edgei]*value;
}
else
{
source_[own[edgei]] -= upper()[edgei]*values[i];
source_[own[edgei]] -= upper()[edgei]*value;
}
upper()[edgei] = 0.0;
......@@ -389,11 +388,11 @@ void Foam::faMatrix<Type>::setValues
{
if (facei == own[edgei])
{
source_[nei[edgei]] -= lower()[edgei]*values[i];
source_[nei[edgei]] -= lower()[edgei]*value;
}
else
{
source_[own[edgei]] -= upper()[edgei]*values[i];
source_[own[edgei]] -= upper()[edgei]*value;
}
upper()[edgei] = 0.0;
......@@ -402,18 +401,15 @@ void Foam::faMatrix<Type>::setValues
}
else
{
label patchi = mesh.boundary().whichPatch(edgei);
const label patchi = mesh.boundary().whichPatch(edgei);
if (internalCoeffs_[patchi].size())
{
label patchEdgei =
const label patchEdgei =
mesh.boundary()[patchi].whichEdge(edgei);
internalCoeffs_[patchi][patchEdgei] =
pTraits<Type>::zero;
boundaryCoeffs_[patchi][patchEdgei] =
pTraits<Type>::zero;
internalCoeffs_[patchi][patchEdgei] = Zero;
boundaryCoeffs_[patchi][patchEdgei] = Zero;
}
}
}
......@@ -422,14 +418,48 @@ void Foam::faMatrix<Type>::setValues
}
template<class Type>
void Foam::faMatrix<Type>::setValues
(
const labelUList& faceLabels,
const Type& value
)
{
this->setValuesFromList(faceLabels, UniformList<Type>(value));
}
template<class Type>
void Foam::faMatrix<Type>::setValues
(
const labelUList& faceLabels,
const UList<Type>& values
)
{
this->setValuesFromList(faceLabels, values);
}
template<class Type>
void Foam::faMatrix<Type>::setValues
(
const labelUList& faceLabels,
const UIndirectList<Type>& values
)
{
this->setValuesFromList(faceLabels, values);
}
template<class Type>
void Foam::faMatrix<Type>::setReference
(
const label facei,
const Type& value
const Type& value,
const bool forceReference
)
{
if (psi_.needReference())
if ((forceReference || psi_.needReference()) && facei >= 0)
{
if (Pstream::master())
{
......@@ -440,6 +470,52 @@ void Foam::faMatrix<Type>::setReference
}
template<class Type>
void Foam::faMatrix<Type>::setReferences
(
const labelUList& faceLabels,
const Type& value,
const bool forceReference
)
{
if (forceReference || psi_.needReference())
{
forAll(faceLabels, facei)
{
const label faceId = faceLabels[facei];
if (faceId >= 0)
{
source()[faceId] += diag()[faceId]*value;
diag()[faceId] += diag()[faceId];
}
}
}
}
template<class Type>
void Foam::faMatrix<Type>::setReferences
(
const labelUList& faceLabels,
const UList<Type>& values,
const bool forceReference
)
{
if (forceReference || psi_.needReference())
{
forAll(faceLabels, facei)
{
const label faceId = faceLabels[facei];
if (faceId >= 0)
{
source()[faceId] += diag()[faceId]*values[facei];
diag()[faceId] += diag()[faceId];
}
}
}
}
template<class Type>
void Foam::faMatrix<Type>::relax(const scalar alpha)
{
......
......@@ -56,10 +56,9 @@ Author
namespace Foam
{
// Forward declaration of friend functions and operators
template<class Type>
class faMatrix;
// Forward Declarations
template<class Type> class faMatrix;
template<class T> class UIndirectList;
template<class Type>
Ostream& operator<<(Ostream&, const faMatrix<Type>&);
......@@ -75,9 +74,10 @@ class faMatrix
public refCount,
public lduMatrix
{
// Private data
// Private Data
// Reference to GeometricField<Type, faPatchField, areaMesh>
//- Const reference to field
// Converted into a non-const reference at the point of solution.
const GeometricField<Type, faPatchField, areaMesh>& psi_;
//- Dimension set
......@@ -87,20 +87,25 @@ class faMatrix
Field<Type> source_;
//- Boundary scalar field containing pseudo-matrix coeffs
// for internal faces
//- for internal faces
FieldField<Field, Type> internalCoeffs_;
//- Boundary scalar field containing pseudo-matrix coeffs
// for boundary faces
//- for boundary faces
FieldField<Field, Type> boundaryCoeffs_;
//- Face flux field for non-orthogonal correction
mutable GeometricField<Type, faePatchField, edgeMesh>
*faceFluxCorrectionPtr_;
// Private member functions
protected:
//- Declare friendship with the faSolver class
friend class faSolver;
// Protected Member Functions
//- Add patch contribution to internal field
template<class Type2>
......@@ -154,9 +159,21 @@ class faMatrix
) const;
// Matrix manipulation functionality
//- Set solution in given faces to the specified values
template<template<class> class ListType>
void setValuesFromList
(
const labelUList& faceLabels,
const ListType<Type>& values
);
public:
//- Solver class returned by the solver function
//- used for systems in which it is useful to cache the solver for reuse.
class faSolver
{
faMatrix<Type>& faMat_;
......@@ -174,11 +191,11 @@ public:
{}
// Member functions
// Member Functions
//- Solve returning the solution statistics.
// Solver controls read from dictionary
SolverPerformance<Type> solve(const dictionary&);
SolverPerformance<Type> solve(const dictionary& solverControls);
//- Solve returning the solution statistics.
// Solver controls read from faSolution
......@@ -194,18 +211,18 @@ public:
//- Construct given a field to solve for
faMatrix
(
const GeometricField<Type, faPatchField, areaMesh>&,
const dimensionSet&
const GeometricField<Type, faPatchField, areaMesh>& psi,
const dimensionSet& ds
);
//- Construct as copy
//- Copy construct
faMatrix(const faMatrix<Type>&);
//- Construct from Istream given field to solve for
faMatrix
(
const GeometricField<Type, faPatchField, areaMesh>&,
Istream&
const GeometricField<Type, faPatchField, areaMesh>& psi,
Istream& is
);
//- Clone
......@@ -216,7 +233,7 @@ public:
virtual ~faMatrix();
// Member functions
// Member Functions
// Access
......@@ -241,14 +258,28 @@ public:
}
//- faBoundary scalar field containing pseudo-matrix coeffs
// for internal cells
//- for internal cells
const FieldField<Field, Type>& internalCoeffs() const
{
return internalCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for internal cells
FieldField<Field, Type>& internalCoeffs()
{
return internalCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
// for boundary cells
//- for boundary cells
const FieldField<Field, Type>& boundaryCoeffs() const
{
return boundaryCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for boundary cells
FieldField<Field, Type>& boundaryCoeffs()
{
return boundaryCoeffs_;
......@@ -268,23 +299,56 @@ public:
// Operations
//- Set solution in given cells and eliminate corresponding
// equations from the matrix
//- Set solution in given faces to the specified value
//- and eliminate the corresponding equations from the matrix.
void setValues
(
const labelUList& faceLabels,
const Type& value
);
//- Set solution in given faces to the specified values
//- and eliminate the corresponding equations from the matrix.
void setValues
(
const labelUList& faces,
const labelUList& faceLabels,
const UList<Type>& values
);
//- Set solution in given faces to the specified values
//- and eliminate the corresponding equations from the matrix.
void setValues
(
const labelUList& faceLabels,
const UIndirectList<Type>& values
);
//- Set reference level for solution
void setReference
(
const label facei,
const Type& value
const Type& value,
const bool forceReference = false
);
//- Set reference level for solution
void setReferences
(
const labelUList& faceLabels,
const Type& value,
const bool forceReference = false
);
//- Set reference level for solution
void setReferences
(
const labelUList& faceLabels,
const UList<Type>& values,
const bool forceReference = false
);
//- Set reference level for a component of the solution
// on a given patch face
//- on a given patch face
void setComponentReference
(
const label patchi,
......@@ -297,6 +361,7 @@ public:
// alpha = 1 : diagonally equal
// alpha < 1 : ,, dominant
// alpha = 0 : do nothing
// Note: Requires positive diagonal.
void relax(const scalar alpha);
//- Relax matrix (for steady-state solution).
......@@ -327,7 +392,7 @@ public:
tmp<GeometricField<Type, faePatchField, edgeMesh>> flux() const;
// Member operators
// Member Operators
void operator=(const faMatrix<Type>&);
void operator=(const tmp<faMatrix<Type>>&);
......@@ -355,7 +420,7 @@ public:
void operator*=(const dimensioned<scalar>&);
// Ostream operator
// Ostream Operator
friend Ostream& operator<< <Type>
(
......@@ -365,7 +430,7 @@ public:
};
// * * * * * * * * * * * * * * * Global functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
template<class Type>
void checkMethod
......@@ -399,7 +464,7 @@ SolverPerformance<Type> solve(faMatrix<Type>&, Istream&);
//- Solve returning the solution statistics given convergence tolerance,
// deleting temporary matrix after solution.
//- deleting temporary matrix after solution.
// Solver controls read Istream
template<class Type>
SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&, Istream&);
......@@ -412,13 +477,13 @@ SolverPerformance<Type> solve(faMatrix<Type>&);
//- Solve returning the solution statistics given convergence tolerance,
// deleting temporary matrix after solution.
//- deleting temporary matrix after solution.
// Solver controls read faSolution
template<class Type>
SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&);
// * * * * * * * * * * * * * * * Global operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class Type>
tmp<faMatrix<Type>> operator-
......@@ -780,6 +845,8 @@ tmp<faMatrix<Type>> operator*
#include "faMatrix.C"
#endif
// Specialisation for scalars
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
......
......@@ -60,10 +60,13 @@ Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve
<< "solving faMatrix<Type>"
<< endl;
auto& psi =
const_cast<GeometricField<Type, faPatchField, areaMesh>&>(psi_);
SolverPerformance<Type> solverPerfVec
(
"faMatrix<Type>::solve",
psi_.name()
psi.name()
);
scalarField saveDiag(diag());
......@@ -75,10 +78,6 @@ Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve
lduInterfaceFieldPtrsList interfaces =
psi_.boundaryField().scalarInterfaces();
// Cast into a non-const to solve
GeometricField<Type, faPatchField, areaMesh>& psi =
const_cast<GeometricField<Type, faPatchField, areaMesh>&>(psi_);
for (direction cmpt = 0; cmpt < Type::nComponents; ++cmpt)
{
// copy field and source
......@@ -179,14 +178,11 @@ Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve()
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::faMatrix<Type>::residual() const
{
tmp<Field<Type>> tres(source_);
tmp<Field<Type>> tres(new Field<Type>(source_));
Field<Type>& res = tres().ref();
addBoundarySource(res);
lduInterfaceFieldPtrsList interfaces =
psi_.boundaryField().scalarInterfaces();
// Loop over field components
for (direction cmpt = 0; cmpt < Type::nComponents; ++cmpt)
{
......@@ -208,7 +204,7 @@ Foam::tmp<Foam::Field<Type>> Foam::faMatrix<Type>::residual() const
psiCmpt,
res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
bouCoeffsCmpt,
interfaces,
psi_.boundaryField().scalarInterfaces(),
cmpt
)
);
......
......@@ -32,6 +32,7 @@ License
#include "extrapolatedCalculatedFvPatchFields.H"
#include "coupledFvPatchFields.H"
#include "UIndirectList.H"
#include "UniformList.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
......@@ -207,12 +208,8 @@ void Foam::fvMatrix<Type>::setValuesFromList
if (symmetric() || asymmetric())
{
const cell& c = cells[celli];
forAll(c, j)
for (const label facei : cells[celli])
{
const label facei = c[j];
if (mesh.isInternalFace(facei))
{
if (symmetric())
......@@ -245,18 +242,15 @@ void Foam::fvMatrix<Type>::setValuesFromList
}
else
{
label patchi = mesh.boundaryMesh().whichPatch(facei);
const label patchi = mesh.boundaryMesh().whichPatch(facei);
if (internalCoeffs_[patchi].size())
{
label patchFacei =
const label patchFacei =
mesh.boundaryMesh()[patchi].whichFace(facei);
internalCoeffs_[patchi][patchFacei] =
Zero;
boundaryCoeffs_[patchi][patchFacei] =
Zero;
internalCoeffs_[patchi][patchFacei] = Zero;
boundaryCoeffs_[patchi][patchFacei] = Zero;
}
}
}
......@@ -310,10 +304,10 @@ Foam::fvMatrix<Type>::fvMatrix
}
// Update the boundary coefficients of psi without changing its event No.
GeometricField<Type, fvPatchField, volMesh>& psiRef =
const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
auto& psiRef =
const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
label currentStatePsi = psiRef.eventNo();
const label currentStatePsi = psiRef.eventNo();
psiRef.boundaryFieldRef().updateCoeffs();
psiRef.eventNo() = currentStatePsi;
}
......@@ -336,11 +330,11 @@ Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
if (fvm.faceFluxCorrectionPtr_)