fvScalarMatrix.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-2017 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 "fvScalarMatrix.H"
31#include "profiling.H"
32#include "PrecisionAdaptor.H"
34#include "cyclicPolyPatch.H"
35#include "cyclicAMIPolyPatch.H"
36
37// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
38
39template<>
41(
42 const label patchi,
43 const label facei,
44 const direction,
45 const scalar value
46)
47{
48 if (psi_.needReference())
49 {
50 if (Pstream::master())
51 {
52 internalCoeffs_[patchi][facei] +=
53 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
54
55 boundaryCoeffs_[patchi][facei] +=
56 diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
57 *value;
58 }
59 }
60}
61
62
63template<>
66(
67 const dictionary& solverControls
68)
69{
71 if (psi_.mesh().name() != polyMesh::defaultRegion)
72 {
73 regionName = psi_.mesh().name() + "::";
74 }
75 addProfiling(solve, "fvMatrix::solve." + regionName + psi_.name());
76
77 if (debug)
78 {
79 Info.masterStream(this->mesh().comm())
80 << "fvMatrix<scalar>::solver(const dictionary& solverControls) : "
81 "solver for fvMatrix<scalar>"
82 << endl;
83 }
84
85 scalarField saveDiag(diag());
86 addBoundaryDiag(diag(), 0);
87
88 lduInterfaceFieldPtrsList interfaces =
89 psi_.boundaryField().scalarInterfaces();
90
91
93 (
95 (
96 *this,
97 lduMatrix::solver::New
98 (
99 psi_.name(),
100 *this,
101 boundaryCoeffs_,
102 internalCoeffs_,
103 interfaces,
104 solverControls
105 )
106 )
107 );
108
109 diag() = saveDiag;
110
111 return solverPtr;
112}
113
114
115template<>
117(
118 const dictionary& solverControls
119)
120{
121 const int logLevel =
122 solverControls.getOrDefault<int>
123 (
124 "log",
125 solverPerformance::debug
126 );
127
128 auto& psi =
130 (
131 fvMat_.psi()
132 );
133
134 scalarField saveDiag(fvMat_.diag());
135 fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
136
137 scalarField totalSource(fvMat_.source());
138 fvMat_.addBoundarySource(totalSource, false);
139
140 // Assign new solver controls
141 solver_->read(solverControls);
142
143 solverPerformance solverPerf = solver_->solve
144 (
146 totalSource
147 );
148
149 if (logLevel)
150 {
151 solverPerf.print(Info.masterStream(fvMat_.mesh().comm()));
152 }
153
154 fvMat_.diag() = saveDiag;
155
157
158 psi.mesh().setSolverPerformance(psi.name(), solverPerf);
159
160 return solverPerf;
161}
162
163
164template<>
166(
167 const dictionary& solverControls
168)
169{
170 if (debug)
171 {
172 Info.masterStream(this->mesh().comm())
173 << "fvMatrix<scalar>::solveSegregated"
174 "(const dictionary& solverControls) : "
175 "solving fvMatrix<scalar>"
176 << endl;
177 }
178
179 const int logLevel =
180 solverControls.getOrDefault<int>
181 (
182 "log",
183 solverPerformance::debug
184 );
185
186 scalarField saveLower;
187 scalarField saveUpper;
188
189 if (useImplicit_)
190 {
192
193 if (psi_.mesh().fluxRequired(psi_.name()))
194 {
195 // Save lower/upper for flux calculation
196 if (asymmetric())
197 {
198 saveLower = lower();
199 }
200 saveUpper = upper();
201 }
202
206 direction cmpt = 0;
207 manipulateMatrix(cmpt);
208 }
209
210 scalarField saveDiag(diag());
211 addBoundaryDiag(diag(), 0);
212
213 scalarField totalSource(source_);
214 addBoundarySource(totalSource, false);
215
216 lduInterfaceFieldPtrsList interfaces;
217 PtrDynList<lduInterfaceField> newInterfaces;
218 if (!useImplicit_)
219 {
220 interfaces = this->psi(0).boundaryField().scalarInterfaces();
221 }
222 else
223 {
224 setInterfaces(interfaces, newInterfaces);
225 }
226
227 tmp<scalarField> tpsi;
228 if (!useImplicit_)
229 {
230 tpsi.ref
231 (
233 (
234 psi_
236 );
237 }
238 else
239 {
240 tpsi = tmp<scalarField>::New(lduAddr().size(), Zero);
241 scalarField& psi = tpsi.ref();
242
243 for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
244 {
245 const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
246 const auto& psiInternal = this->psi(fieldi).primitiveField();
247
248 forAll(psiInternal, localCellI)
249 {
250 psi[cellOffset + localCellI] = psiInternal[localCellI];
251 }
252 }
253 }
254 scalarField& psi = tpsi.ref();
255
256 // Solver call
258 (
259 this->psi(0).name(),
260 *this,
261 boundaryCoeffs_,
262 internalCoeffs_,
263 interfaces,
264 solverControls
265 )->solve(psi, totalSource);
266
267 if (useImplicit_)
268 {
269 for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
270 {
271 auto& psiInternal =
273 (
274 this->psi(fieldi)
276
277 const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
278
279 forAll(psiInternal, localCellI)
280 {
281 psiInternal[localCellI] = psi[localCellI + cellOffset];
282 }
283 }
284 }
285
286 if (logLevel)
287 {
288 solverPerf.print(Info.masterStream(mesh().comm()));
289 }
290
291 diag() = saveDiag;
292
293 if (useImplicit_)
294 {
295 if (psi_.mesh().fluxRequired(psi_.name()))
296 {
297 // Restore lower/upper
298 if (asymmetric())
299 {
300 lower().setSize(saveLower.size());
301 lower() = saveLower;
302 }
303
304 upper().setSize(saveUpper.size());
305 upper() = saveUpper;
306 }
307 // Set the original lduMesh
308 setLduMesh(psi_.mesh());
309 }
310
311 for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
312 {
313 auto& localPsi =
315 (
316 this->psi(fieldi)
317 );
318
319 localPsi.correctBoundaryConditions();
320 localPsi.mesh().setSolverPerformance(localPsi.name(), solverPerf);
321 }
322
323 return solverPerf;
324}
325
326
327template<>
329{
330 scalarField boundaryDiag(psi_.size(), Zero);
331 addBoundaryDiag(boundaryDiag, 0);
332
333 const scalarField& psif = psi_.primitiveField();
335 const solveScalarField& psi = tpsi();
336
338 (
340 (
341 psi,
342 source_ - boundaryDiag*psif,
343 boundaryCoeffs_,
345 0
346 )
347 );
348
350 addBoundarySource(tres_s.ref());
351
352 return tres_s;
353}
354
355
356template<>
358{
360 (
362 (
364 (
365 "H("+psi_.name()+')',
366 psi_.instance(),
367 psi_.mesh(),
370 ),
371 psi_.mesh(),
372 dimensions_/dimVol,
373 extrapolatedCalculatedFvPatchScalarField::typeName
374 )
375 );
376 volScalarField& Hphi = tHphi.ref();
377
378 Hphi.primitiveFieldRef() = (lduMatrix::H(psi_.primitiveField()) + source_);
380
381 Hphi.primitiveFieldRef() /= psi_.mesh().V();
383
384 return tHphi;
385}
386
387
388template<>
390{
392 (
394 (
396 (
397 "H(1)",
398 psi_.instance(),
399 psi_.mesh(),
402 ),
403 psi_.mesh(),
404 dimensions_/(dimVol*psi_.dimensions()),
405 extrapolatedCalculatedFvPatchScalarField::typeName
406 )
407 );
408 volScalarField& H1_ = tH1.ref();
409
411 //addBoundarySource(Hphi.primitiveField());
412
413 H1_.primitiveFieldRef() /= psi_.mesh().V();
415
416 return tH1;
417}
418
419
420// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
const dimensionSet & dimensions() const
Return dimensions.
const Mesh & mesh() const
Return mesh.
lduInterfaceFieldPtrsList scalarInterfaces() const
Generic GeometricField class.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
const fileName & instance() const noexcept
Read access to instance path component.
Definition: IOobjectI.H:196
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A dynamically resizable PtrList with allocation management.
Definition: PtrDynList.H:64
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.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
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
void setSolverPerformance(const word &name, const SolverPerformance< Type > &sp) const
Add/set the solverPerformance entry for the named field.
Definition: dataTemplates.C:36
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< GeometricField< Type, faPatchField, areaMesh > > H() const
Return the H operation source.
Definition: faMatrix.C:633
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.
void manipulateMatrix(direction cmp)
Manipulate matrix.
Definition: fvMatrix.C:790
void transferFvMatrixCoeffs()
Transfer lower, upper, diag and source to this fvMatrix.
Definition: fvMatrix.C:813
tmp< volScalarField > H1() const
Return H(1)
Definition: fvMatrix.C:1394
void createOrUpdateLduPrimitiveAssembly()
Create or update ldu assembly.
Definition: fvMatrix.C:906
label nMatrices() const
Definition: fvMatrix.H:339
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
void setInterfaces(lduInterfaceFieldPtrsList &, PtrDynList< lduInterfaceField > &newInterfaces)
Set interfaces.
Definition: fvMatrix.C:470
tmp< GeometricField< Type, fvPatchField, volMesh > > H() const
Return the H operation source.
Definition: fvMatrix.C:1332
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.
lduPrimitiveMeshAssembly * lduMeshPtr()
Access to lduPrimitiveMeshAssembly.
Definition: fvMatrix.C:880
void setBounAndInterCoeffs()
Manipulate boundary/internal coeffs for coupling.
Definition: fvMatrix.C:660
tmp< scalarField > H1() const
scalarField & upper()
Definition: lduMatrix.C:203
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:578
void setLduMesh(const lduMesh &m)
Set the LDU mesh containing the addressing is obtained.
Definition: lduMatrix.H:572
const lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:566
scalarField & lower()
Definition: lduMatrix.C:174
scalarField & diag()
Definition: lduMatrix.C:192
bool asymmetric() const
Definition: lduMatrix.H:633
const labelList & cellOffsets() const
Return cellOffsets.
OSstream & masterStream(const label communicator)
T & ref() const
Definition: refPtrI.H:203
bool fluxRequired(const word &name) const
Set flux-required for given name (mutable)
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)
A scalar instance of fvMatrix.
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 > &)
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:64
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
conserve primitiveFieldRef()+