sensitivityBezierFIIncompressible.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-2020 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(
48);
49
50
51// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
52
54{
55 // Laplace solution controls
56 const dictionary dxdbDict = dict_.subOrEmptyDict("dxdbSolver");
57 meshMovementIters_ = dxdbDict.getOrDefault<label>("iters", 1000);
59 dxdbDict.getOrDefault<scalar>("tolerance", 1.e-07);
60
61 // Read variables related to the adjoint eikonal solver
63}
64
65
67(
68 const label iCP,
69 const label idir
70)
71{
72 read();
74 volVectorField& m = tm.ref();
75
76 // SOLVE FOR DXDB
77 //~~~~~~~~~~~~~~~~
78 // set boundary conditions
79 for (const label patchI : sensitivityPatchIDs_)
80 {
81 // interpolate parameterization info to faces
82 tmp<vectorField> tdxidXjFace = Bezier_.dxdbFace(patchI, iCP, idir);
83 const vectorField& dxidXjFace = tdxidXjFace();
84
85 m.boundaryFieldRef()[patchI] == dxidXjFace;
86 }
87
88 // iterate the adjoint to the eikonal equation
89 for (label iter = 0; iter < meshMovementIters_; iter++)
90 {
91 Info<< "Mesh Movement Propagation(direction, CP), ("
92 << idir << ", " << iCP << "), Iteration : "<< iter << endl;
93
95 (
97 );
98
99 // Scalar residual = max(mEqn.solve().initialResidual());
100 scalar residual = mag(mEqn.solve().initialResidual());
101
102 Info<< "Max dxdb " << gMax(mag(m)()) << endl;
103
105
106 // Check convergence
107 if (residual < meshMovementResidualLimit_)
108 {
109 Info<< "\n***Reached dxdb convergence limit, iteration " << iter
110 << "***\n\n";
111 break;
112 }
113 }
114
115 return tm;
116}
117
118
119// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
120
122(
123 const fvMesh& mesh,
124 const dictionary& dict,
126)
127:
129 //Bezier_(mesh, dict), // AJH Read locally?
130 Bezier_(mesh, mesh.lookupObject<IOdictionary>("optimisationDict")),
131 flowSens_(3*Bezier_.nBezier(), Zero),
132 dSdbSens_(3*Bezier_.nBezier(), Zero),
133 dndbSens_(3*Bezier_.nBezier(), Zero),
134 dxdbDirectSens_(3*Bezier_.nBezier(), Zero),
135 dVdbSens_(3*Bezier_.nBezier(), Zero),
136 distanceSens_(3*Bezier_.nBezier(), Zero),
137 optionsSens_(3*Bezier_.nBezier(), Zero),
138 bcSens_(3*Bezier_.nBezier(), Zero),
139
140 derivativesFolder_("optimisation"/type() + "Derivatives"),
141
142 meshMovementIters_(-1),
143 meshMovementResidualLimit_(1.e-7),
144 dxdb_
145 (
146 variablesSet::autoCreateMeshMovementField
147 (
148 mesh,
149 "mTilda",
151 )
152 )
153{
154 read();
155
157 // Create folder to store sensitivities
159}
160
161
162// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
163
165{
166 // Adjoint to the eikonal equation
167 autoPtr<volTensorField> distanceSensPtr(nullptr);
169 {
170 // Solver equation
171 eikonalSolver_->solve();
172
173 // Allocate memory and compute grad(dxdb) multiplier
174 distanceSensPtr.reset
175 (
176 createZeroFieldPtr<tensor>
177 (
178 mesh_,
179 "distanceSensPtr",
180 dimensionSet(0, 2, -3, 0, 0, 0, 0)
181 )
182 );
183 distanceSensPtr() = eikonalSolver_->getFISensitivityTerm()().T();
184 }
185
186 const label nBezier = Bezier_.nBezier();
187 const label nDVs = 3*nBezier;
188 for (label iDV = 0; iDV < nDVs; iDV++)
189 {
190 label iCP = iDV%nBezier;
191 label idir = iDV/nBezier;
192 if
193 (
194 (idir == 0 && Bezier_.confineXmovement()[iCP])
195 || (idir == 1 && Bezier_.confineYmovement()[iCP])
196 || (idir == 2 && Bezier_.confineZmovement()[iCP])
197 )
198 {
199 continue;
200 }
201 else
202 {
203 // Flow term
204 // ~~~~~~~~~~~
205 // compute dxdb and its grad
207 const volVectorField& m = tm();
208 volTensorField gradDxDb(fvc::grad(m, "grad(dxdb)"));
209
210 flowSens_[iDV] =
211 gSum
212 (
214 * mesh_.V()
215 );
216
217 for (const label patchI : sensitivityPatchIDs_)
218 {
219 // Contribution from objective function
220 // term from delta(n dS)/delta b and
221 // term from delta(n)/delta b
222 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
223 tmp<vectorField> tdSdb =
224 Bezier_.dndbBasedSensitivities(patchI, iCP, idir);
225 const vectorField& dSdb = tdSdb();
226 tmp<vectorField> tdndb =
227 Bezier_.dndbBasedSensitivities(patchI, iCP, idir, false);
228 const vectorField& dndb = tdndb();
229 dSdbSens_[iDV] += gSum(dSfdbMult_()[patchI] & dSdb);
230 dndbSens_[iDV] += gSum(dnfdbMult_()[patchI] & dndb);
231
232 // Contribution from objective function
233 // term from delta( x ) / delta b
234 // Only for objectives directly including
235 // x, like moments
236 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
237 tmp<vectorField> tdxdbFace =
238 Bezier_.dxdbFace(patchI, iCP, idir);
239 const vectorField& dxdbFace = tdxdbFace();
240 dxdbDirectSens_[iDV] +=
241 gSum(dxdbDirectMult_()[patchI] & dxdbFace);
242
243 // Contribution from boundary conditions
244 bcSens_[iDV] += gSum(bcDxDbMult_()[patchI] & dxdbFace);
245 }
246
247 // Contribution from delta (V) / delta b
248 // For objectives defined as volume integrals only
249 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250 dVdbSens_[iDV] =
251 gSum
252 (
254 * fvc::div(m)().primitiveField()
255 * mesh_.V()
256 );
257
258 // Distance dependent term
259 //~~~~~~~~~~~~~~~~~~~~~~~~~
261 {
262 distanceSens_[iDV] =
263 gSum
264 (
265 (
266 distanceSensPtr().primitiveField()
267 && gradDxDb.primitiveField()
268 )
269 *mesh_.V()
270 );
271 }
272
273 // Terms from fvOptions
274 optionsSens_[iDV] +=
276 }
277
278 // Sum contributions
281 + dSdbSens_
282 + dndbSens_
284 + dVdbSens_
287 + bcSens_;
288 }
289}
290
291
293{
294 flowSens_ = Zero;
295 dSdbSens_ = Zero;
296 dndbSens_ = Zero;
298 dVdbSens_ = Zero;
301 bcSens_ = Zero;
302
304}
305
306
307void sensitivityBezierFI::write(const word& baseName)
308{
309 Info<< "Writing control point sensitivities to file" << endl;
310 if (Pstream::master())
311 {
312 OFstream derivFile
313 (
315 baseName + adjointVars_.solverName() + mesh_.time().timeName()
316 );
317 unsigned int widthDV = max(int(name(flowSens_.size()).size()), int(3));
318 unsigned int width = IOstream::defaultPrecision() + 7;
319 derivFile
320 << setw(widthDV) << "#dv" << " "
321 << setw(width) << "total" << " "
322 << setw(width) << "flow" << " "
323 << setw(width) << "dSdb" << " "
324 << setw(width) << "dndb" << " "
325 << setw(width) << "dxdbDirect" << " "
326 << setw(width) << "dVdb" << " "
327 << setw(width) << "distance" << " "
328 << setw(width) << "options" << " "
329 << setw(width) << "dvdb" << endl;
330 const label nDVs = derivatives_.size();
331 const label nBezier = Bezier_.nBezier();
332 const boolListList& confineMovement = Bezier_.confineMovement();
333 label lastActive(-1);
334
335 for (label iDV = 0; iDV < nDVs; iDV++)
336 {
337 const label iCP(iDV%nBezier);
338 const label idir(iDV/nBezier);
339 if (!confineMovement[idir][iCP])
340 {
341 if (iDV!=lastActive + 1)
342 {
343 derivFile << "\n";
344 }
345 lastActive = iDV;
346 derivFile
347 << setw(widthDV) << iDV << " "
348 << setw(width) << derivatives_[iDV] << " "
349 << setw(width) << flowSens_[iDV] << " "
350 << setw(width) << dSdbSens_[iDV] << " "
351 << setw(width) << dndbSens_[iDV] << " "
352 << setw(width) << dxdbDirectSens_[iDV] << " "
353 << setw(width) << dVdbSens_[iDV] << " "
354 << setw(width) << distanceSens_[iDV] << " "
355 << setw(width) << optionsSens_[iDV] << " "
356 << setw(width) << bcSens_[iDV] << endl;
357 }
358 }
359 }
360}
361
362
363// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
364
365} // End namespace incompressible
366} // End namespace Foam
367
368// ************************************************************************* //
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.
label nBezier() const
Number of Bezier control points.
Definition: Bezier.C:127
const boolList & confineZmovement() const
Confine z movement.
Definition: Bezier.C:151
const boolList & confineYmovement() const
Confine y movement.
Definition: Bezier.C:145
const boolListList & confineMovement() const
Info about confining movement in all directions.
Definition: Bezier.C:157
const boolList & confineXmovement() const
Confine x movement.
Definition: Bezier.C:139
tmp< tensorField > dxdbFace(const label patchI, const label cpI, bool useChainRule=true) const
dxdb tensor for a Bezier parameterized patch
Definition: Bezier.C:279
tmp< tensorField > dndbBasedSensitivities(const label patchI, const label cpI, bool returnDimensionedNormalSens=true) const
Definition: Bezier.C:164
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
static unsigned int defaultPrecision() noexcept
Return the default precision.
Definition: IOstream.H:342
Output to file stream, using an OSstream.
Definition: OFstream.H:57
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
Ostream & printExecutionTime(OSstream &os) const
Print the elapsed ExecutionTime (cpu-time), ClockTime.
Definition: TimeIO.C:618
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Base class for adjoint solvers.
Definition: adjointSolver.H:60
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
Definition: dictionary.C:540
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Definition: dimensionSet.H:109
virtual bool write()
Write the output fields.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
Base class for incompressibleAdjoint solvers.
Base class for Field Integral-based sensitivity derivatives.
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
bool includeDistance_
Include distance variation in sens computation.
volTensorField gradDxDbMult_
grad(dx/db) multiplier
autoPtr< adjointEikonalSolver > eikonalSolver_
Adjoint eikonal equation solver.
void read()
Read options and update solver pointers if necessary.
scalarField divDxDbMult_
div(dx/db) multiplier
vectorField optionsDxDbMult_
dx/db multiplier coming from fvOptions
Abstract base class for adjoint-based sensitivities in incompressible flows.
Calculation of adjoint based sensitivities for Bezier control points using the FI appoach.
virtual void clearSensitivities()
Zero sensitivity fields and their constituents.
scalarField dSdbSens_
Term depending on delta(n dS)/delta b.
tmp< volVectorField > solveMeshMovementEqn(const label iCP, const label idir)
scalarField dVdbSens_
Term depending on delta(V)/delta b.
virtual void assembleSensitivities()
Assemble sensitivities.
scalarField optionsSens_
Term depending on fvOptions.
scalarField distanceSens_
Term depending on distance differentiation.
scalarField bcSens_
Term depending on the differenation of boundary conditions.
scalarField dndbSens_
Term depending on delta(n)/delta b.
autoPtr< boundaryVectorField > dSfdbMult_
Fields related to direct sensitivities.
const fvMesh & mesh_
Definition: sensitivity.H:69
dictionary dict_
Definition: sensitivity.H:70
splitCell * master() const
Definition: splitCell.H:113
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
Base class for creating a set of variables.
Definition: variablesSet.H:50
const word & solverName() const
Return solver name.
Definition: variablesSet.C:84
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
const volScalarField & T
dynamicFvMesh & mesh
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
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)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:83
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
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
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
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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
Type gMax(const FieldField< Field, Type > &f)
dictionary dict