sensitivityVolBSplinesIncompressible.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) 2007-2020 PCOpt/NTUA
9 Copyright (C) 2013-2020 FOSS GP
10 Copyright (C) 2019 OpenCFD Ltd.
11-------------------------------------------------------------------------------
12License
13 This file is part of OpenFOAM.
14
15 OpenFOAM is free software: you can redistribute it and/or modify it
16 under the terms of the GNU General Public License as published by
17 the Free Software Foundation, either version 3 of the License, or
18 (at your option) any later version.
19
20 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27
28\*---------------------------------------------------------------------------*/
29
32#include "IOmanip.H"
33
34// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35
36namespace Foam
37{
38
39namespace incompressible
40{
41
42// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43
46(
50);
51
52// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
53
55{
57 {
58 label passedCPs = 0;
60 forAll(boxes, iNURB)
61 {
62 label nb = boxes[iNURB].getControlPoints().size();
63 for (label cpI = 0; cpI < nb; cpI++)
64 {
65 vector dSdbSensCP(Zero);
66 vector dndbSensCP(Zero);
67 for (const label patchI : sensitivityPatchIDs_)
68 {
69 tensorField dSdb
70 (
71 boxes[iNURB].dndbBasedSensitivities(patchI, cpI)
72 );
73 dSdbSensCP += gSum(dSfdbMult_()[patchI] & dSdb);
74
75 tensorField dndb
76 (
77 boxes[iNURB].dndbBasedSensitivities
78 (
79 patchI,
80 cpI,
81 false
82 )
83 );
84 dndbSensCP += gSum((dnfdbMult_()[patchI] & dndb));
85 }
86 dSdbSens_[passedCPs + cpI] = dSdbSensCP;
87 dndbSens_[passedCPs + cpI] = dndbSensCP;
88 }
89 passedCPs += nb;
90 }
93
94 passedCPs = 0;
95 forAll(boxes, iNURB)
96 {
97 vectorField sensDxDbDirect =
98 boxes[iNURB].computeControlPointSensitivities
99 (
101 sensitivityPatchIDs_.toc()
102 );
103
104 // Transfer to global list
105 forAll(sensDxDbDirect, cpI)
106 {
107 dxdbDirectSens_[passedCPs + cpI] = sensDxDbDirect[cpI];
108 }
109 passedCPs += sensDxDbDirect.size();
110 }
112 }
113}
114
115
117{
118 label passedCPs = 0;
120 forAll(boxes, iNURB)
121 {
122 vectorField sensBcsDxDb =
123 boxes[iNURB].computeControlPointSensitivities
124 (
125 bcDxDbMult_(),
126 sensitivityPatchIDs_.toc()
127 );
128
129 // Transfer to global list
130 forAll(sensBcsDxDb, cpI)
131 {
132 bcSens_[passedCPs + cpI] = sensBcsDxDb[cpI];
133 }
134 passedCPs += sensBcsDxDb.size();
135 }
137}
138
139
140// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
141
143(
144 const fvMesh& mesh,
145 const dictionary& dict,
147)
148:
150 volBSplinesBase_
151 (
153 ),
154
155 flowSens_(0),
156 dSdbSens_(0),
157 dndbSens_(0),
158 dxdbDirectSens_(0),
159 bcSens_(0),
160
161 derivativesFolder_("optimisation"/type() + "Derivatives")
162{
163 // No boundary field pointers need to be allocated
165 derivatives_ = scalarField(3*nCPs, Zero);
166 flowSens_ = vectorField(nCPs, Zero);
167 dSdbSens_ = vectorField(nCPs, Zero);
168 dndbSens_ = vectorField(nCPs, Zero);
170 bcSens_ = vectorField(nCPs, Zero);
171
172 // Create folder to store sensitivities
174}
175
176
177// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
178
180{
181 // Assemble the sensitivity map
182 // Solves for the post-processing equations and adds their contribution to
183 // the sensitivity map
185
186 // Finalise sensitivities including dxFace/db
187 const boundaryVectorField& faceSens =
188 surfaceSensitivity_.getWallFaceSensVecBoundary();
189
190 label passedCPs(0);
192 forAll(boxes, iNURB)
193 {
194 vectorField sens =
195 boxes[iNURB].computeControlPointSensitivities
196 (
197 faceSens,
198 sensitivityPatchIDs_.toc()
199 );
200 // Transfer to global list
201 forAll(sens, cpI)
202 {
203 flowSens_[passedCPs + cpI] = sens[cpI];
204 }
205 passedCPs += sens.size();
206 }
208
209 // Contribution from objective function
210 // Note:
211 // includeObjectiveContribution has to be set to false (false by default)
212 // in surfaceSensitivity, in order to avoid computing this term twice.
213 // Optionally avoided altogether if includeObjectiveContribution is set to
214 // false for sensitivityVolBSplines
216
218
219 // Transform sensitivites to scalarField in order to cooperate with
220 // updateMethod
221 forAll(flowSens_, cpI)
222 {
223 derivatives_[3*cpI] =
224 flowSens_[cpI].x()
225 + dSdbSens_[cpI].x()
226 + dndbSens_[cpI].x()
227 + dxdbDirectSens_[cpI].x()
228 + bcSens_[cpI].x();
229 derivatives_[3*cpI + 1] =
230 flowSens_[cpI].y()
231 + dSdbSens_[cpI].y()
232 + dndbSens_[cpI].y()
233 + dxdbDirectSens_[cpI].y()
234 + bcSens_[cpI].y();
235 derivatives_[3*cpI + 2] =
236 flowSens_[cpI].z()
237 + dSdbSens_[cpI].z()
238 + dndbSens_[cpI].z()
239 + dxdbDirectSens_[cpI].z()
240 + bcSens_[cpI].z();
241 }
242}
243
244
246{
252
254}
255
256
258{
259 Info<< "Writing control point sensitivities to file" << endl;
260 // Write sensitivity map
261 SIBase::write(baseName);
262 // Write control point sensitivities
263 if (Pstream::master())
264 {
265 OFstream derivFile
266 (
268 baseName + adjointVars_.solverName() + mesh_.time().timeName()
269 );
270 unsigned int widthDV =
271 max(int(Foam::name(derivatives_.size()).size()), int(3));
272 unsigned int width = IOstream::defaultPrecision() + 7;
273 derivFile
274 << setw(widthDV) << "#cp" << " "
275 << setw(width) << "total::x"<< " "
276 << setw(width) << "total::y"<< " "
277 << setw(width) << "total::z"<< " "
278 << setw(width) << "flow::x" << " "
279 << setw(width) << "flow::y" << " "
280 << setw(width) << "flow::z" << " "
281 << setw(width) << "dSdb::x" << " "
282 << setw(width) << "dSdb::y" << " "
283 << setw(width) << "dSdb::z" << " "
284 << setw(width) << "dndb::x" << " "
285 << setw(width) << "dndb::y" << " "
286 << setw(width) << "dndb::z" << " "
287 << setw(width) << "dxdbDirect::x" << " "
288 << setw(width) << "dxdbDirect::y" << " "
289 << setw(width) << "dxdbDirect::z" << " "
290 << setw(width) << "dvdb::x" << " "
291 << setw(width) << "dvdb::y" << " "
292 << setw(width) << "dvdb::z" << endl;
293
294 label passedCPs(0);
295 label lastActive(-1);
297 forAll(boxes, iNURB)
298 {
299 label nb = boxes[iNURB].getControlPoints().size();
300 const boolList& activeCPs = boxes[iNURB].getActiveCPs();
301 for (label iCP = 0; iCP < nb; iCP++)
302 {
303 if (activeCPs[iCP])
304 {
305 label globalCP = passedCPs + iCP;
306 if (globalCP!=lastActive + 1)
307 {
308 derivFile << "\n";
309 }
310 lastActive = globalCP;
311
312 derivFile
313 << setw(widthDV) << globalCP << " "
314 << setw(width) << derivatives_[3*globalCP] << " "
315 << setw(width) << derivatives_[3*globalCP + 1] << " "
316 << setw(width) << derivatives_[3*globalCP + 2] << " "
317 << setw(width) << flowSens_[globalCP].x() << " "
318 << setw(width) << flowSens_[globalCP].y() << " "
319 << setw(width) << flowSens_[globalCP].z() << " "
320 << setw(width) << dSdbSens_[globalCP].x() << " "
321 << setw(width) << dSdbSens_[globalCP].y() << " "
322 << setw(width) << dSdbSens_[globalCP].z() << " "
323 << setw(width) << dndbSens_[globalCP].x() << " "
324 << setw(width) << dndbSens_[globalCP].y() << " "
325 << setw(width) << dndbSens_[globalCP].z() << " "
326 << setw(width) << dxdbDirectSens_[globalCP].x() << " "
327 << setw(width) << dxdbDirectSens_[globalCP].y() << " "
328 << setw(width) << dxdbDirectSens_[globalCP].z() << " "
329 << setw(width) << bcSens_[globalCP].x() << " "
330 << setw(width) << bcSens_[globalCP].y() << " "
331 << setw(width) << bcSens_[globalCP].z()
332 << endl;
333 }
334 }
335 passedCPs += nb;
336 }
337 }
338}
339
340
341// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
342
343} // End namespace incompressible
344} // End namespace Foam
345
346// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Generic GeometricBoundaryField class.
static unsigned int defaultPrecision() noexcept
Return the default precision.
Definition: IOstream.H:342
Output to file stream, using an OSstream.
Definition: OFstream.H:57
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Base class for adjoint solvers.
Definition: adjointSolver.H:60
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
virtual bool write()
Write the output fields.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
Base class for incompressibleAdjoint solvers.
Base class for Surface Integral-based sensitivity derivatives.
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
sensitivitySurface surfaceSensitivity_
Surface sensitivities.
bool includeObjective_
Whether to include direct sensitivities or not.
Abstract base class for adjoint-based sensitivities in incompressible flows.
virtual void assembleSensitivities()
Assemble sensitivities.
Calculation of adjoint based sensitivities at vol B-Splines control points using the SI or e-SI appro...
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
vectorField dSdbSens_
Term depending on delta(n dS)/delta b.
vectorField bcSens_
Term dependng on the differentiation of boundary conditions.
vectorField dndbSens_
Term depending on delta (n)/delta b.
volBSplinesBase & volBSplinesBase_
Reference to underlaying volumetric B-Splines morpher.
autoPtr< boundaryVectorField > dSfdbMult_
Fields related to direct sensitivities.
const fvMesh & mesh_
Definition: sensitivity.H:69
splitCell * master() const
Definition: splitCell.H:113
const word & solverName() const
Return solver name.
Definition: variablesSet.C:84
Class constructing a number of volumetric B-Splines boxes, read from dynamicMeshDict....
PtrList< NURBS3DVolume > & boxesRef()
Get non-const reference to the vol. B-splines boxes.
label getTotalControlPointsNumber() const
Get cumulative number of control points from all boxes.
void boundControlPointMovement(vectorField &controlPointsMovement) const
Bound control points movement.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Type gSum(const FieldField< Field, Type > &f)
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
Definition: MSwindows.C:515
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Field< vector > vectorField
Specialisation of Field<T> for vector.
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
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
A non-counting (dummy) refCount.
Definition: refCount.H:59