fvMatrixSolve.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) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2016-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
27\*---------------------------------------------------------------------------*/
28
29#include "LduMatrix.H"
30#include "diagTensorField.H"
31#include "profiling.H"
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 if (psi_.needReference())
46 {
47 if (Pstream::master())
48 {
49 internalCoeffs_[patchi][facei].component(cmpt) +=
50 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
51
52 boundaryCoeffs_[patchi][facei].component(cmpt) +=
53 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
54 *value;
55 }
56 }
57}
58
59
60template<class Type>
62(
63 const dictionary& solverControls
64)
65{
67 if (psi_.mesh().name() != polyMesh::defaultRegion)
68 {
69 regionName = psi_.mesh().name() + "::";
70 }
71 addProfiling(solve, "fvMatrix::solve." + regionName + psi_.name());
72
73 if (debug)
74 {
75 Info.masterStream(this->mesh().comm())
76 << "fvMatrix<Type>::solveSegregatedOrCoupled"
77 "(const dictionary& solverControls) : "
78 "solving fvMatrix<Type>"
79 << endl;
80 }
81
82 // Do not solve if maxIter == 0
83 if (solverControls.getOrDefault<label>("maxIter", -1) == 0)
84 {
86 }
87
88 word type(solverControls.getOrDefault<word>("type", "segregated"));
89
90 if (type == "segregated")
91 {
92 return solveSegregated(solverControls);
93 }
94 else if (type == "coupled")
95 {
96 return solveCoupled(solverControls);
97 }
98 else
99 {
100 FatalIOErrorInFunction(solverControls)
101 << "Unknown type " << type
102 << "; currently supported solver types are segregated and coupled"
103 << exit(FatalIOError);
104
106 }
107}
108
109
110template<class Type>
112(
113 const dictionary& solverControls
114)
115{
116 if (useImplicit_)
117 {
119 << "Implicit option is not allowed for type: " << Type::typeName
120 << exit(FatalError);
121 }
122
123 if (debug)
124 {
125 Info.masterStream(this->mesh().comm())
126 << "fvMatrix<Type>::solveSegregated"
127 "(const dictionary& solverControls) : "
128 "solving fvMatrix<Type>"
129 << endl;
130 }
131
132 const int logLevel =
133 solverControls.getOrDefault<int>
134 (
135 "log",
137 );
138
139 auto& psi =
141
142 SolverPerformance<Type> solverPerfVec
143 (
144 "fvMatrix<Type>::solveSegregated",
145 psi.name()
146 );
147
148 scalarField saveDiag(diag());
149
150 Field<Type> source(source_);
151
152 // At this point include the boundary source from the coupled boundaries.
153 // This is corrected for the implicit part by updateMatrixInterfaces within
154 // the component loop.
155 addBoundarySource(source);
156
157 typename Type::labelType validComponents
158 (
159 psi.mesh().template validComponents<Type>()
160 );
161
162 for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
163 {
164 if (validComponents[cmpt] == -1) continue;
165
166 // copy field and source
167
168 scalarField psiCmpt(psi.primitiveField().component(cmpt));
169 addBoundaryDiag(diag(), cmpt);
170
171 scalarField sourceCmpt(source.component(cmpt));
172
173 FieldField<Field, scalar> bouCoeffsCmpt
174 (
175 boundaryCoeffs_.component(cmpt)
176 );
177
178 FieldField<Field, scalar> intCoeffsCmpt
179 (
180 internalCoeffs_.component(cmpt)
181 );
182
183 lduInterfaceFieldPtrsList interfaces =
184 psi.boundaryField().scalarInterfaces();
185
186 // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
187 // bouCoeffsCmpt for the explicit part of the coupled boundary
188 // conditions
189 {
190 PrecisionAdaptor<solveScalar, scalar> sourceCmpt_ss(sourceCmpt);
192
193 const label startRequest = Pstream::nRequests();
194
195 initMatrixInterfaces
196 (
197 true,
198 bouCoeffsCmpt,
199 interfaces,
200 psiCmpt_ss(),
201 sourceCmpt_ss.ref(),
202 cmpt
203 );
204
205 updateMatrixInterfaces
206 (
207 true,
208 bouCoeffsCmpt,
209 interfaces,
210 psiCmpt_ss(),
211 sourceCmpt_ss.ref(),
212 cmpt,
213 startRequest
214 );
215 }
216
217 solverPerformance solverPerf;
218
219 // Solver call
220 solverPerf = lduMatrix::solver::New
221 (
222 psi.name() + pTraits<Type>::componentNames[cmpt],
223 *this,
224 bouCoeffsCmpt,
225 intCoeffsCmpt,
226 interfaces,
227 solverControls
228 )->solve(psiCmpt, sourceCmpt, cmpt);
229
230 if (logLevel)
231 {
232 solverPerf.print(Info.masterStream(this->mesh().comm()));
233 }
234
235 solverPerfVec.replace(cmpt, solverPerf);
236 solverPerfVec.solverName() = solverPerf.solverName();
237
238 psi.primitiveFieldRef().replace(cmpt, psiCmpt);
239 diag() = saveDiag;
240 }
241
242 psi.correctBoundaryConditions();
243
244 psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);
245
246 return solverPerfVec;
247}
248
249
250template<class Type>
252(
253 const dictionary& solverControls
254)
255{
256 if (debug)
257 {
258 Info.masterStream(this->mesh().comm())
259 << "fvMatrix<Type>::solveCoupled"
260 "(const dictionary& solverControls) : "
261 "solving fvMatrix<Type>"
262 << endl;
263 }
264
265 const int logLevel =
266 solverControls.getOrDefault<int>
267 (
268 "log",
270 );
271
272 auto& psi =
274
275 LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
276 coupledMatrix.diag() = diag();
277 coupledMatrix.upper() = upper();
278 coupledMatrix.lower() = lower();
279 coupledMatrix.source() = source();
280
281 addBoundaryDiag(coupledMatrix.diag(), 0);
282 addBoundarySource(coupledMatrix.source(), false);
283
284 coupledMatrix.interfaces() = psi.boundaryFieldRef().interfaces();
285 coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
286 coupledMatrix.interfacesLower() = internalCoeffs().component(0);
287
289 coupledMatrixSolver
290 (
292 (
293 psi.name(),
294 coupledMatrix,
295 solverControls
296 )
297 );
298
299 SolverPerformance<Type> solverPerf
300 (
301 coupledMatrixSolver->solve(psi)
302 );
303
304 if (logLevel)
305 {
306 solverPerf.print(Info.masterStream(this->mesh().comm()));
307 }
308
309 psi.correctBoundaryConditions();
310
311 psi.mesh().setSolverPerformance(psi.name(), solverPerf);
312
313 return solverPerf;
314}
315
316
317template<class Type>
319(
320 const dictionary& solverControls
321)
322{
323 return psi_.mesh().solve(*this, solverControls);
324}
325
326
327template<class Type>
330{
331 return solver(solverDict());
332}
333
334
335template<class Type>
337{
338 return solve(fvMat_.solverDict());
339}
340
341
342template<class Type>
344{
345 return this->solve(solverDict());
346}
347
348
349template<class Type>
351{
352 tmp<Field<Type>> tres(new Field<Type>(source_));
353 Field<Type>& res = tres.ref();
354
356
357 // Loop over field components
358 for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
359 {
360 scalarField psiCmpt(psi_.primitiveField().component(cmpt));
361
362 scalarField boundaryDiagCmpt(psi_.size(), Zero);
363 addBoundaryDiag(boundaryDiagCmpt, cmpt);
364
365 FieldField<Field, scalar> bouCoeffsCmpt
366 (
367 boundaryCoeffs_.component(cmpt)
368 );
369
370 res.replace
371 (
372 cmpt,
374 (
375 psiCmpt,
376 res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
377 bouCoeffsCmpt,
379 cmpt
380 )
381 );
382 }
383
384 return tres;
385}
386
387
388// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
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.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
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.
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:88
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
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
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
tmp< Field< Type > > residual() const
Return the matrix residual.
SolverPerformance< Type > solve()
Solve returning the solution statistics.
autoPtr< fvSolver > solver()
Construct and return the solver.
SolverPerformance< Type > solve()
Solve returning the solution statistics.
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: fvMatrix.C:177
void setComponentReference(const label patchi, const label facei, const direction cmpt, const scalar value)
Definition: fvMatrixSolve.C:38
SolverPerformance< Type > solveSegregatedOrCoupled(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:62
SolverPerformance< Type > solveCoupled(const dictionary &)
Solve coupled returning the solution statistics.
tmp< Field< Type > > residual() const
Return the matrix residual.
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: fvMatrix.C:124
SolverPerformance< Type > solveSegregated(const dictionary &)
Solve segregated returning the solution statistics.
const dictionary & solverDict() const
Return the solver dictionary taking into account finalIteration.
Definition: fvMatrix.C:1558
Abstract base-class for lduMatrix solvers.
Definition: lduMatrix.H:99
OSstream & masterStream(const label communicator)
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
A class for handling words, derived from Foam::string.
Definition: word.H:68
const volScalarField & psi
dynamicFvMesh & mesh
Foam::word regionName(Foam::polyMesh::defaultRegion)
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
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
IOerror FatalIOError
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
CEqn solve()
static const char *const typeName
The type name used in ensight case files.