faMatrixSolve.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016-2017 Wikki Ltd
9 Copyright (C) 2019-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Description
28 Finite-Area matrix basic solvers.
29
30\*---------------------------------------------------------------------------*/
31
32#include "PrecisionAdaptor.H"
33
34// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
35
36template<class Type>
38(
39 const label patchi,
40 const label facei,
41 const direction cmpt,
42 const scalar value
43)
44{
45 internalCoeffs_[patchi][facei].component(cmpt) +=
46 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
47
48 boundaryCoeffs_[patchi][facei].component(cmpt) +=
49 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]*value;
50}
51
52
53template<class Type>
55(
56 const dictionary& solverControls
57)
58{
60 << "solving faMatrix<Type>"
61 << endl;
62
63 const int logLevel =
64 solverControls.getOrDefault<int>
65 (
66 "log",
68 );
69
70 auto& psi =
72
73 SolverPerformance<Type> solverPerfVec
74 (
75 "faMatrix<Type>::solve",
76 psi.name()
77 );
78
79 scalarField saveDiag(diag());
80
81 Field<Type> source(source_);
82 addBoundarySource(source);
83
84 // Note: make a copy of interfaces: no longer a reference
85 lduInterfaceFieldPtrsList interfaces =
86 psi_.boundaryField().scalarInterfaces();
87
88 for (direction cmpt = 0; cmpt < Type::nComponents; ++cmpt)
89 {
90 // copy field and source
91
92 scalarField psiCmpt(psi_.primitiveField().component(cmpt));
93 addBoundaryDiag(diag(), cmpt);
94
95 scalarField sourceCmpt(source.component(cmpt));
96
97 FieldField<Field, scalar> bouCoeffsCmpt
98 (
99 boundaryCoeffs_.component(cmpt)
100 );
101
102 FieldField<Field, scalar> intCoeffsCmpt
103 (
104 internalCoeffs_.component(cmpt)
105 );
106
107 // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
108 // bouCoeffsCmpt for the explicit part of the coupled boundary
109 // conditions
110 {
111 PrecisionAdaptor<solveScalar, scalar> sourceCmpt_ss(sourceCmpt);
113
114 const label startRequest = Pstream::nRequests();
115
116 initMatrixInterfaces
117 (
118 true,
119 bouCoeffsCmpt,
120 interfaces,
121 psiCmpt_ss(),
122 sourceCmpt_ss.ref(),
123 cmpt
124 );
125
126 updateMatrixInterfaces
127 (
128 true,
129 bouCoeffsCmpt,
130 interfaces,
131 psiCmpt_ss(),
132 sourceCmpt_ss.ref(),
133 cmpt,
134 startRequest
135 );
136 }
137
138 solverPerformance solverPerf;
139
140 // Solver call
141 solverPerf = lduMatrix::solver::New
142 (
143 psi_.name() + pTraits<Type>::componentNames[cmpt],
144 *this,
145 bouCoeffsCmpt,
146 intCoeffsCmpt,
147 interfaces,
148 solverControls
149 )->solve(psiCmpt, sourceCmpt, cmpt);
150
151 if (logLevel)
152 {
153 solverPerf.print(Info);
154 }
155
156 solverPerfVec.replace(cmpt, solverPerf);
157 solverPerfVec.solverName() = solverPerf.solverName();
158
159 psi.primitiveFieldRef().replace(cmpt, psiCmpt);
160 diag() = saveDiag;
161 }
162
163 psi.correctBoundaryConditions();
164
165 psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);
166
167 return solverPerfVec;
168}
169
170
171template<class Type>
173{
174 return solve(faMat_.psi().mesh().solverDict(faMat_.psi().name()));
175}
176
177
178template<class Type>
180{
181 return solve(this->psi().mesh().solverDict(this->psi().name()));
182}
183
184
185template<class Type>
187{
188 tmp<Field<Type>> tres(new Field<Type>(source_));
189 Field<Type>& res = tres().ref();
190
192
193 // Loop over field components
194 for (direction cmpt = 0; cmpt < Type::nComponents; ++cmpt)
195 {
196 scalarField psiCmpt(psi_.internalField().component(cmpt));
197
198 scalarField boundaryDiagCmpt(psi_.size(), Zero);
199 addBoundaryDiag(boundaryDiagCmpt, cmpt);
200
201 FieldField<Field, scalar> bouCoeffsCmpt
202 (
203 boundaryCoeffs_.component(cmpt)
204 );
205
206 res.replace
207 (
208 cmpt,
210 (
211 psiCmpt,
212 res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
213 bouCoeffsCmpt,
215 cmpt
216 )
217 );
218 }
219
220 return tres;
221}
222
223
224// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
tmp< DimensionedField< cmptType, GeoMesh > > component(const direction d) const
Return a component field of the field.
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:80
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:557
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:545
lduInterfaceFieldPtrsList scalarInterfaces() const
Generic GeometricField class.
void replace(const direction d, const GeometricField< cmptType, PatchField, GeoMesh > &gcf)
Replace specified field component with content from another field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const Internal & internalField() const
Return a const-reference to the dimensioned internal field.
A non-const Field/List wrapper with possible data conversion.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
void print(Ostream &os) const
Print summary of solver performance to the given stream.
const word & solverName() const noexcept
Return solver name.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
SolverPerformance< Type > solve()
Solve returning the solution statistics.
SolverPerformance< Type > solve()
Solve returning the solution statistics.
const GeometricField< Type, faPatchField, areaMesh > & psi() const
Definition: faMatrix.H:270
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: faMatrix.C:151
void setComponentReference(const label patchi, const label facei, const direction cmpt, const scalar value)
Definition: faMatrixSolve.C:38
tmp< Field< Type > > residual() const
Return the matrix residual.
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: faMatrix.C:117
const lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:566
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
A traits class, which is primarily used for primitives.
Definition: pTraits.H:59
T & ref() const
Definition: refPtrI.H:203
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
const volScalarField & psi
#define DebugInFunction
Report an information message using Foam::Info.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)