setForwardADSeeds.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v3
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #ifdef CODI_AD_FORWARD
9 
10 Info << endl
11  << "Forward-Mode AD is activated, setting the AD seeds..." << endl;
12 
13 if (daOptionPtr_->getAllOptions().subDict("useAD").getWord("mode") == "forward")
14 {
15 
16  word designVarName = daOptionPtr_->getAllOptions().subDict("useAD").getWord("dvName");
17 
18  dictionary designVarDict = daOptionPtr_->getAllOptions().subDict("designVar");
19  // get the subDict for this dvName
20  dictionary dvSubDict = designVarDict.subDict(designVarName);
21  word type = dvSubDict.getWord("designVarType");
22 
23  if (type == "FFD")
24  {
25  // NOTE: seedIndex is not needed for FFD here, instead, we set FFD2XvSeedVec
26  // in pyDAFoam.py and then read it here
27  // for Xv derivs
28  pointField meshPointsForwardAD = mesh.points();
29 
30  Info << "Setting seeds based on FFD2XvSeedVec: " << endl;
31  PetscScalar* seedArray;
32  VecGetArray(FFD2XvSeedVec_, &seedArray);
33  forAll(meshPointsForwardAD, i)
34  {
35  for (label j = 0; j < 3; j++)
36  {
37  label localIdx = daIndexPtr_->getLocalXvIndex(i, j);
38  PetscScalar seedVal = seedArray[localIdx];
39  meshPointsForwardAD[i][j].setGradient(seedVal);
40  }
41  }
42  VecRestoreArray(FFD2XvSeedVec_, &seedArray);
43 
44  mesh.movePoints(meshPointsForwardAD);
45  }
46  else if (type == "Field")
47  {
48  label seedIndex = daOptionPtr_->getAllOptions().subDict("useAD").getLabel("seedIndex");
49 
50  word fieldName = dvSubDict.getWord("fieldName");
51  word fieldType = dvSubDict.getWord("fieldType");
52 
53  if (fieldType == "scalar")
54  {
55  volScalarField& state = const_cast<volScalarField&>(
56  meshPtr_->thisDb().lookupObject<volScalarField>(fieldName));
57 
58  forAll(state, cellI)
59  {
60  label glxIdx = daIndexPtr_->getGlobalCellIndex(cellI);
61  if (glxIdx == seedIndex || seedIndex == -1)
62  {
63  state[cellI].setGradient(1.0);
64  }
65  }
66  }
67  else if (fieldType == "vector")
68  {
69  volVectorField& state = const_cast<volVectorField&>(
70  meshPtr_->thisDb().lookupObject<volVectorField>(fieldName));
71 
72  forAll(state, cellI)
73  {
74  for (label i = 0; i < 3; i++)
75  {
76  label glxIdx = daIndexPtr_->getLocalCellVectorIndex(cellI, i);
77  if (glxIdx == seedIndex || seedIndex == -1)
78  {
79  state[cellI][i].setGradient(1.0);
80  }
81  }
82  }
83  }
84  else
85  {
86  FatalErrorIn("") << "fieldType not valid. Options: scalar or vector"
87  << abort(FatalError);
88  }
89  }
90  else if (type == "ACTD")
91  {
92  // for ACTD derivs
93  scalarList actDVListForwardAD(9);
94 
95  label seedIndex = daOptionPtr_->getAllOptions().subDict("useAD").getLabel("seedIndex");
96 
97  DAFvSource& fvSource = const_cast<DAFvSource&>(
98  meshPtr_->thisDb().lookupObject<DAFvSource>("DAFvSource"));
99 
100  word diskName = dvSubDict.getWord("actuatorName");
101  dictionary fvSourceSubDict = daOptionPtr_->getAllOptions().subDict("fvSource");
102  word source = fvSourceSubDict.subDict(diskName).getWord("source");
103  if (source == "cylinderAnnulusSmooth")
104  {
105  // get the design variable vals
106  for (label i = 0; i < 9; i++)
107  {
108  if (i == seedIndex)
109  {
110  actDVListForwardAD[i] = fvSource.getActuatorDVs(diskName, i);
111  actDVListForwardAD[i].setGradient(1.0);
112  fvSource.setActuatorDVs(diskName, i, actDVListForwardAD[i]);
113  }
114  }
115  // the actuatorDVs are updated, now we need to recompute fvSource
116  // this is not needed for the residual partials because fvSource
117  // will be automatically calculated in the UEqn, but for the
118  // obj partials, we need to manually recompute fvSource
119  fvSource.updateFvSource();
120  }
121  else
122  {
123  FatalErrorIn("") << "only support cylinderAnnulusSmooth!"
124  << abort(FatalError);
125  }
126  }
127  else if (type == "AOA")
128  {
129  // for AOA derivs
130  scalar aoaForwardAD;
131 
132  // get info from dvSubDict. This needs to be defined in the pyDAFoam
133  // name of the boundary patch
134  wordList patches;
135  dvSubDict.readEntry<wordList>("patches", patches);
136  // the streamwise axis of aoa, aoa = tan( U_normal/U_flow )
137  word flowAxis = dvSubDict.getWord("flowAxis");
138  word normalAxis = dvSubDict.getWord("normalAxis");
139 
140  HashTable<label> axisIndices;
141  axisIndices.set("x", 0);
142  axisIndices.set("y", 1);
143  axisIndices.set("z", 2);
144  label flowAxisIndex = axisIndices[flowAxis];
145  label normalAxisIndex = axisIndices[normalAxis];
146 
147  volVectorField& U = const_cast<volVectorField&>(
148  meshPtr_->thisDb().lookupObject<volVectorField>("U"));
149 
150  // now we need to get the Ux, Uy values from the inout patches
151  scalar Ux0 = -1e16, Uy0 = -1e16;
152  forAll(patches, idxI)
153  {
154  word patchName = patches[idxI];
155  label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
156  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
157  {
158  if (U.boundaryField()[patchI].type() == "fixedValue")
159  {
160  Uy0 = U.boundaryField()[patchI][0][normalAxisIndex];
161  Ux0 = U.boundaryField()[patchI][0][flowAxisIndex];
162  break;
163  }
164  else if (U.boundaryField()[patchI].type() == "inletOutlet")
165  {
166  mixedFvPatchField<vector>& inletOutletPatch =
167  refCast<mixedFvPatchField<vector>>(U.boundaryFieldRef()[patchI]);
168  Uy0 = inletOutletPatch.refValue()[0][normalAxisIndex];
169  Ux0 = inletOutletPatch.refValue()[0][flowAxisIndex];
170  break;
171  }
172  else
173  {
174  FatalErrorIn("") << "boundaryType: " << U.boundaryFieldRef()[patchI].type()
175  << " not supported!"
176  << "Available options are: fixedValue, inletOutlet"
177  << abort(FatalError);
178  }
179  }
180  }
181  // need to reduce the U value across all processors, this is because some of
182  // the processors might not own the prescribed patches so their U value will be still -1e16, but
183  // when calling the following reduce function, they will get the correct U from other processors
184  reduce(Ux0, maxOp<scalar>());
185  reduce(Uy0, maxOp<scalar>());
186  aoaForwardAD = atan(Uy0 / Ux0);
187 
188  // convert to degree
189  aoaForwardAD *= 180.0 / constant::mathematical::pi.getValue();
190 
191  Info << "Setting the seed for AOA: " << aoaForwardAD << endl;
192  aoaForwardAD.setGradient(1.0);
193 
194  // convert to back to radian
195  aoaForwardAD *= constant::mathematical::pi.getValue() / 180.0;
196 
197  // assign the angle of attack to the static variable in DAUtility. It will be used in determing
198  // the forceDir in DAObjeFuncForce.C
199  DAUtility::angleOfAttackRadForwardAD = aoaForwardAD;
200 
201  // set far field Ux, Uy
202  forAll(patches, idxI)
203  {
204  word patchName = patches[idxI];
205  label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
206 
207  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
208  {
209  scalar UMag = sqrt(Ux0 * Ux0 + Uy0 * Uy0);
210  scalar UxNew = UMag * cos(aoaForwardAD);
211  scalar UyNew = UMag * sin(aoaForwardAD);
212 
213  if (U.boundaryField()[patchI].type() == "fixedValue")
214  {
215  forAll(U.boundaryField()[patchI], faceI)
216  {
217  U.boundaryFieldRef()[patchI][faceI][flowAxisIndex] = UxNew;
218  U.boundaryFieldRef()[patchI][faceI][normalAxisIndex] = UyNew;
219  }
220  }
221  else if (U.boundaryField()[patchI].type() == "inletOutlet")
222  {
223  mixedFvPatchField<vector>& inletOutletPatch =
224  refCast<mixedFvPatchField<vector>>(U.boundaryFieldRef()[patchI]);
225 
226  forAll(U.boundaryField()[patchI], faceI)
227  {
228  inletOutletPatch.refValue()[faceI][flowAxisIndex] = UxNew;
229  inletOutletPatch.refValue()[faceI][normalAxisIndex] = UyNew;
230  }
231  }
232  }
233  }
234  }
235  else if (type == "BC")
236  {
237  scalar BCForwardAD = -1e16;
238 
239  if (designVarName == "MRF")
240  {
241  const IOMRFZoneListDF& MRF = meshPtr_->thisDb().lookupObject<IOMRFZoneListDF>("MRFProperties");
242 
243  // first, we get the current value of omega and assign it to BC
244  scalar& omega = const_cast<scalar&>(MRF.getOmegaRef());
245  BCForwardAD = omega;
246 
247  Info << "Setting the seed for BCForwardAD: " << BCForwardAD << endl;
248  BCForwardAD.setGradient(1.0);
249  // now set BC to omega
250  omega = BCForwardAD;
251  }
252  else
253  {
254 
255  // get info from dvSubDict. This needs to be defined in the pyDAFoam
256  // name of the variable for changing the boundary condition
257  word varName = dvSubDict.getWord("variable");
258  // name of the boundary patch
259  wordList patches;
260  dvSubDict.readEntry<wordList>("patches", patches);
261  // the compoent of a vector variable, ignore when it is a scalar
262  label comp = dvSubDict.getLabel("comp");
263 
264  // Now get the BC value
265  forAll(patches, idxI)
266  {
267  word patchName = patches[idxI];
268  label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
269  if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
270  {
271  volVectorField& state(const_cast<volVectorField&>(
272  meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
273  // for decomposed domain, don't set BC if the patch is empty
274  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
275  {
276  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
277  {
278  BCForwardAD = state.boundaryFieldRef()[patchI][0][comp];
279  }
280  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
281  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
282  {
283  mixedFvPatchField<vector>& inletOutletPatch =
284  refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
285  BCForwardAD = inletOutletPatch.refValue()[0][comp];
286  }
287  }
288  }
289  else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
290  {
291  volScalarField& state(const_cast<volScalarField&>(
292  meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
293  // for decomposed domain, don't set BC if the patch is empty
294  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
295  {
296  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
297  {
298  BCForwardAD = state.boundaryFieldRef()[patchI][0];
299  }
300  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
301  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
302  {
303  mixedFvPatchField<scalar>& inletOutletPatch =
304  refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
305  BCForwardAD = inletOutletPatch.refValue()[0];
306  }
307  }
308  }
309  }
310  // need to reduce the BC value across all processors, this is because some of
311  // the processors might not own the prescribed patches so their BC value will be still -1e16, but
312  // when calling the following reduce function, they will get the correct BC from other processors
313  reduce(BCForwardAD, maxOp<scalar>());
314 
315  Info << "Setting the seed for BCForwardAD: " << BCForwardAD << endl;
316  BCForwardAD.setGradient(1.0);
317 
318  // ******* now set BC ******
319  forAll(patches, idxI)
320  {
321  word patchName = patches[idxI];
322  label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
323  if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
324  {
325  volVectorField& state(const_cast<volVectorField&>(
326  meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
327  // for decomposed domain, don't set BC if the patch is empty
328  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
329  {
330  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
331  {
332  forAll(state.boundaryFieldRef()[patchI], faceI)
333  {
334  state.boundaryFieldRef()[patchI][faceI][comp] = BCForwardAD;
335  }
336  }
337  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
338  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
339  {
340  mixedFvPatchField<vector>& inletOutletPatch =
341  refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
342  vector val = inletOutletPatch.refValue()[0];
343  val[comp] = BCForwardAD;
344  inletOutletPatch.refValue() = val;
345  }
346  }
347  }
348  else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
349  {
350  volScalarField& state(const_cast<volScalarField&>(
351  meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
352  // for decomposed domain, don't set BC if the patch is empty
353  if (meshPtr_->boundaryMesh()[patchI].size() > 0)
354  {
355  if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
356  {
357  forAll(state.boundaryFieldRef()[patchI], faceI)
358  {
359  state.boundaryFieldRef()[patchI][faceI] = BCForwardAD;
360  }
361  }
362  else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
363  || state.boundaryFieldRef()[patchI].type() == "outletInlet")
364  {
365  mixedFvPatchField<scalar>& inletOutletPatch =
366  refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
367  inletOutletPatch.refValue() = BCForwardAD;
368  }
369  }
370  }
371  }
372  }
373  }
374 }
375 
376 #endif
U
U
Definition: pEqnPimpleDyM.H:60
forAll
forAll(pseudoP.boundaryField(), patchI)
Definition: solvePseudoPEqn.H:10
MRF
IOMRFZoneListDF & MRF
Definition: createRefsRhoSimple.H:18
fvSource
volScalarField & fvSource
Definition: createRefsHeatTransfer.H:7
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
source
pseudoPEqn source()

Generated by   doxygen 1.8.17