DASimpleFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6  This class is modified from OpenFOAM's source code
7  applications/solvers/incompressible/simpleFoam
8 
9  OpenFOAM: The Open Source CFD Toolbox
10 
11  Copyright (C): 2011-2016 OpenFOAM Foundation
12 
13  OpenFOAM License:
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 
30 #include "DASimpleFoam.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(DASimpleFoam, 0);
38 addToRunTimeSelectionTable(DASolver, DASimpleFoam, dictionary);
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
42  char* argsAll,
43  PyObject* pyOptions)
44  : DASolver(argsAll, pyOptions),
45  simplePtr_(nullptr),
46  pPtr_(nullptr),
47  UPtr_(nullptr),
48  phiPtr_(nullptr),
49  alphaPorosityPtr_(nullptr),
50  laminarTransportPtr_(nullptr),
51  turbulencePtr_(nullptr),
52  daTurbulenceModelPtr_(nullptr),
53  daFvSourcePtr_(nullptr),
54  fvSourcePtr_(nullptr),
55  MRFPtr_(nullptr),
56  PrPtr_(nullptr),
57  PrtPtr_(nullptr),
58  TPtr_(nullptr),
59  alphatPtr_(nullptr)
60 {
61  // check whether the temperature field exists in the 0 folder
62  hasTField_ = DAUtility::isFieldReadable(meshPtr_(), "T", "volScalarField");
63 
64  // get fvSolution and fvSchemes info for fixed-point adjoint
65  const fvSolution& myFvSolution = meshPtr_->thisDb().lookupObject<fvSolution>("fvSolution");
66  if (myFvSolution.found("relaxationFactors"))
67  {
68  if (myFvSolution.subDict("relaxationFactors").found("equations"))
69  {
70  if (myFvSolution.subDict("relaxationFactors").subDict("equations").found("U"))
71  {
72  relaxUEqn_ = myFvSolution.subDict("relaxationFactors").subDict("equations").getScalar("U");
73  }
74  }
75  }
76  solverDictU_ = myFvSolution.subDict("solvers").subDict("U");
77  solverDictP_ = myFvSolution.subDict("solvers").subDict("p");
78 }
79 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
80 
82 {
83  /*
84  Description:
85  Initialize variables for DASolver
86  */
87 
88  Info << "Initializing fields for DASimpleFoam" << endl;
89  Time& runTime = runTimePtr_();
90  fvMesh& mesh = meshPtr_();
92 #include "createFieldsSimple.H"
93 
94  // read the RAS model from constant/turbulenceProperties
95  const word turbModelName(
96  IOdictionary(
97  IOobject(
98  "turbulenceProperties",
99  mesh.time().constant(),
100  mesh,
101  IOobject::MUST_READ,
102  IOobject::NO_WRITE,
103  false))
104  .subDict("RAS")
105  .lookup("RASModel"));
107 
108 #include "createAdjoint.H"
109 
110  // initialize fvSource and compute the source term
111  const dictionary& allOptions = daOptionPtr_->getAllOptions();
112  if (allOptions.subDict("fvSource").toc().size() != 0)
113  {
114  hasFvSource_ = 1;
115  Info << "Initializing DASource" << endl;
116  word sourceName = allOptions.subDict("fvSource").toc()[0];
117  word fvSourceType = allOptions.subDict("fvSource").subDict(sourceName).getWord("type");
119  fvSourceType, mesh, daOptionPtr_(), daModelPtr_(), daIndexPtr_()));
120  }
121 }
122 
124 {
125  /*
126  Description:
127  Call the primal solver to get converged state variables
128  */
129 
130 #include "createRefsSimple.H"
131 #include "createFvOptions.H"
132 
133  // call correctNut, this is equivalent to turbulence->validate();
134  daTurbulenceModelPtr_->updateIntermediateVariables();
135 
136  Info << "\nStarting time loop\n"
137  << endl;
138 
139  while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
140  {
141 
142  if (printToScreen_)
143  {
144  Info << "Time = " << runTime.timeName() << nl << endl;
145  }
146 
147  p.storePrevIter();
148 
149  // --- Pressure-velocity SIMPLE corrector
150  {
151 #include "UEqnSimple.H"
152 #include "pEqnSimple.H"
153  if (hasTField_)
154  {
155 #include "TEqnSimple.H"
156  }
157  }
158 
159  laminarTransport.correct();
161 
162  // calculate all functions
164  // calculate yPlus
166  // compute the regression model and print the feature
167  regModelFail_ = daRegressionPtr_->compute();
168  daRegressionPtr_->printInputInfo(printToScreen_);
169  // print run time
171 
172  runTime.write();
173  }
174 
175  // write the mesh to files
176  mesh.write();
177 
178  // write associated fields such as URel
179  this->writeAssociatedFields();
180 
181  Info << "End\n"
182  << endl;
183 
184  return this->checkPrimalFailure();
185 }
186 
187 // ************ the following are functions for consistent fixed-point adjoint
188 
190  Vec dFdW,
191  Vec psi)
192 {
193  // If adjoint converged, then adjConv = 0
194  // Otherwise, adjConv = 1
195  label adjConv = 1;
196 
197 #ifdef CODI_ADR
198  /*
199  Description:
200  Solve the adjoint using the fixed-point iteration method
201 
202  xvVec: the volume mesh coordinate vector
203 
204  wVec: the state variable vector
205 
206  dFdW: The dF/dW vector
207 
208  psi: The adjoint solution vector
209  */
210 
211  // Here we keep the values of the previous step psi
212  //VecZeroEntries(psi);
213 
214  label printInterval = daOptionPtr_->getOption<label>("printInterval");
215 
216  word adjEqnSolMethod = daOptionPtr_->getOption<word>("adjEqnSolMethod");
217 
218  if (adjEqnSolMethod == "fixedPoint")
219  {
220  Info << "Solving the adjoint using consistent fixed-point iteration method..."
221  << " Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
222  ;
223  label fpMaxIters = daOptionPtr_->getSubDictOption<label>("adjEqnOption", "fpMaxIters");
224  scalar fpRelTol = daOptionPtr_->getSubDictOption<scalar>("adjEqnOption", "fpRelTol");
225  scalar fpMinResTolDiff = daOptionPtr_->getSubDictOption<scalar>("adjEqnOption", "fpMinResTolDiff");
226 
227  const objectRegistry& db = meshPtr_->thisDb();
228  volVectorField& U = const_cast<volVectorField&>(db.lookupObject<volVectorField>("U"));
229  volScalarField& p = const_cast<volScalarField&>(db.lookupObject<volScalarField>("p"));
230  surfaceScalarField& phi = const_cast<surfaceScalarField&>(db.lookupObject<surfaceScalarField>("phi"));
231  volScalarField& nuTilda = const_cast<volScalarField&>(db.lookupObject<volScalarField>("nuTilda"));
232  // nuEff is only used to construct pseudoUEqn
233  volScalarField nuEff = daTurbulenceModelPtr_->nuEff();
234 
235  // We now pass dFdW directly into adjRes components at the beginning of adjoint residual calculation
236  /*
237  // Initialize field cmpts for dFdW as all-zero
238  volVectorField dFdU("dFdU", 0.0 * U);
239  volScalarField dFdP("dFdP", 0.0 * p);
240  surfaceScalarField dFdPhi("dFdPhi", 0.0 * phi);
241  volScalarField dFdNuTilda("dFdNuTilda", 0.0 * nuTilda);
242  */
243 
244  // Initialize primal residuals as all-zero
245  volVectorField URes("URes", 0.0 * U);
246  volScalarField pRes("pRes", 0.0 * p);
247  surfaceScalarField phiRes("phiRes", 0.0 * phi);
248  volScalarField nuTildaRes("nuTildaRes", 0.0 * nuTilda);
249 
250  // We now pass dFdW directly into adjRes components at the beginning of adjoint residual calculation
251  /*
252  // Pass the values in dFdW to its field cmpts
253  this->vec2Fields("vec2Field", dFdW, dFdU, dFdP, dFdPhi, dFdNuTilda);
254  */
255  /*
256  // Flip the sign for now to avoid the strange issue
257  dFdU = -dFdU;
258  dFdP = -dFdP;
259  dFdPhi = -dFdPhi;
260  dFdNuTilda = -dFdNuTilda;
261  */
262 
263  // Initialize all-zero cmpts for the adjoint vector
264  volVectorField UPsi("UPsi", 0.0 * U);
265  volScalarField pPsi("pPsi", 0.0 * p);
266  surfaceScalarField phiPsi("phiPsi", 0.0 * phi);
267  volScalarField nuTildaPsi("nuTildaPsi", 0.0 * nuTilda);
268 
269  // Pass psiVec to cmpts of the adjoint vector:
270  //this->vec2Fields("vec2Field", psi, UPsi, pPsi, phiPsi, nuTildaPsi);
271 
272  /*
273  // Change every element in the adjoint vector to one
274  forAll(UPsi, cellI)
275  {
276  UPsi[cellI][0] = 1.0;
277  UPsi[cellI][1] = 1.0;
278  UPsi[cellI][2] = 1.0;
279  }
280  forAll(pPsi, cellI)
281  {
282  pPsi[cellI] = 1.0;
283  }
284  forAll(phiPsi.primitiveFieldRef(), faceI)
285  {
286  phiPsi.primitiveFieldRef()[faceI] = 1.0;
287  }
288  forAll(phiPsi.boundaryFieldRef(), patchI)
289  {
290  forAll(phiPsi.boundaryFieldRef()[patchI], faceI)
291  {
292  phiPsi.boundaryFieldRef()[patchI][faceI] = 1.0;
293  }
294  }
295  forAll(nuTildaPsi, cellI)
296  {
297  nuTildaPsi[cellI] = 1.0;
298  }
299  */
300 
301  // Initialize all-zero cmpts for the adjoint residual
302  // Note: we tried defining below as alias of pseudo variables, it works but it actually increases mem usage..
303  volVectorField adjURes("adjURes", 0.0 * U);
304  volScalarField adjPRes("adjPRes", 0.0 * p);
305  surfaceScalarField adjPhiRes("adjPhiRes", 0.0 * phi);
306  volScalarField adjNuTildaRes("adjNuTildaRes", 0.0 * nuTilda);
307 
308  // Initialize the initial L2norm for adjURes, adjPRes, adjPhiRes, adjNuTildaRes as all-zero
309  vector initNormAdjURes = vector::zero;
310  scalar initNormAdjPRes = 0.0;
311  scalar initNormAdjPhiRes = 0.0;
312  scalar initNormAdjNuTildaRes = 0.0;
313 
314  // Below use the norm of dFdW as the initial residual norm
315  /*
316  // Pass -dfdw to adjRes:
317  this->vec2Fields("vec2Field", dFdW, adjURes, adjPRes, adjPhiRes, adjNuTildaRes);
318  adjURes = -adjURes;
319  adjPRes = -adjPRes;
320  adjPhiRes = -adjPhiRes;
321  adjNuTildaRes = -adjNuTildaRes;
322 
323  // Initialize the initial L2norm for adjURes, adjPRes, adjPhiRes, adjNuTildaRes
324  vector initNormAdjURes = L2norm(adjURes);
325  //scalar initNormAdjPRes = L2norm(adjPRes);
326  //scalar initNormAdjPhiRes = L2norm(adjPhiRes);
327  //scalar initNormAdjNuTildaRes = L2norm(adjNuTildaRes);
328 
329  // Get the norm of the total adjoint residual
330  scalar initNormTotAdjRes = sqr(initNormAdjURes[0]) + sqr(initNormAdjURes[1]) + sqr(initNormAdjURes[2]) + sqr(L2norm(adjPRes)) + sqr(L2norm(adjPhiRes)) + sqr(L2norm(adjNuTildaRes));
331  initNormTotAdjRes = sqrt(initNormTotAdjRes);
332  */
333 
334 // Construct the pseudoEqns using pseudoEqns.H
335 // Better to make pseudoUEqn, pseudoPEqn, pseudoNuTildaEqn class variables
336 // Also need to move pseudoNuTildaEqn to DASpalartAllmarasFv3 under DATurbulenceModel
337 #include "pseudoEqns.H"
338 
339  // Initialize the duplicated variables for AD
340  volScalarField p1("p1", p);
341  volScalarField p2("p2", p);
342  this->meshPtr_->setFluxRequired(p2.name());
343  surfaceScalarField phi1("phi1", phi);
344  volVectorField U1("U1", U);
345 
346  // Initialize the intermidiate variables for AD
347  volVectorField gradP("gradP", fvc::grad(p));
348  surfaceScalarField pEqnFlux("pEqnFlux", 0.0 * phi);
349  volScalarField divPhi("divPhi", fvc::div(phi));
350  surfaceScalarField fluxU("fluxU", 0.0 * phi);
351 
352  // No need to print the un-normalized primal residuals
353  /*
354  // Print primal residuals before starting the adjoint iterations
355  this->calcLduResiduals(URes, pRes, phiRes);
356  daTurbulenceModelPtr_->calcLduResidualTurb(nuTildaRes);
357  Info << "Residual for simpleFOAM after convergence: " << endl;
358  Info << "L2 norm of URes: " << L2norm(URes) << endl;
359  Info << "L2 norm of pRes: " << L2norm(pRes) << endl;
360  Info << "L2 norm of phiRes: " << L2norm(phiRes) << endl;
361  Info << "L2 norm of nuTildaRes: " << L2norm(nuTildaRes) << endl;
362  */
363 
364  label cnt = 0;
365  while (cnt < fpMaxIters)
366  {
367  if (cnt % printInterval == 0)
368  {
369  Info << "Step = " << cnt << " Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
370  }
371 
372  // ************************************************************************* //
373  // Step-1: Get D^(-T).adjURes, adjPRes, adjPhiRes, adjNuTildaRes, and update phiPsi
374  // ************************************************************************* //
375  // Step-1.1: Calculate adjURes, adjPRes, adjPhiRes, adjNuTildaRes, and update phiPsi
376 
377  // At beginning of iteration, call calcAdjointResidual once to get the adjoint residual
378  this->calcAdjointResidual(URes, pRes, phiRes, nuTildaRes, dFdW, UPsi, pPsi, phiPsi, nuTildaPsi, adjURes, adjPRes, adjPhiRes, adjNuTildaRes, cnt);
379 
380  // Below uses the norm of total adjoint residual to check convergence
381  /*
382  // Calculate L2 norm for all components of adjoint residual, then check if fpRelTol is met
383  vector normAdjURes = L2norm(adjURes);
384  scalar normAdjPRes = L2norm(adjPRes);
385  scalar normAdjPhiRes = L2norm(adjPhiRes);
386  scalar normAdjNuTildaRes = L2norm(adjNuTildaRes);
387 
388  scalar normTotAdjRes = sqr(normAdjURes[0]) + sqr(normAdjURes[1]) + sqr(normAdjURes[2]) + sqr(normAdjPRes) + sqr(normAdjPhiRes) + sqr(normAdjNuTildaRes);
389  normTotAdjRes = sqrt(normTotAdjRes);
390 
391  // Normalize total residual norm:
392  normTotAdjRes /= initNormTotAdjRes;
393 
394  if (cnt % printInterval == 0)
395  {
396  Info << "L2 norm of adjURes: " << normAdjURes[0] << " " << normAdjURes[1] << " " << normAdjURes[2] << endl;
397  Info << "L2 norm of adjPRes: " << normAdjPRes << endl;
398  Info << "L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
399  Info << "L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
400  }
401 
402  // Check if fpRelTol is met:
403  if (normTotAdjRes < fpRelTol)
404  {
405  Info << "Residual drop of " << fpRelTol << " has been achieved after " << cnt << " steps! Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;;
406  Info << "Final L2 norm of adjURes: " << normAdjURes[0] << " " << normAdjURes[1] << " " << normAdjURes[2] << endl;
407  Info << "Final L2 norm of adjPRes: " << normAdjPRes << endl;
408  Info << "Final L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
409  Info << "Final L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
410  break;
411  }
412  */
413 
414  // Calculate L2 norm for all components of adjoint residual, then check if fpRelTol is met
415  if (cnt >= 1)
416  {
417  vector normAdjURes = L2norm(adjURes);
418  scalar normAdjPRes = L2norm(adjPRes);
419  scalar normAdjPhiRes = L2norm(adjPhiRes);
420  scalar normAdjNuTildaRes = L2norm(adjNuTildaRes);
421 
422  if (cnt == 1)
423  {
424  initNormAdjURes = normAdjURes;
425  initNormAdjPRes = normAdjPRes;
426  initNormAdjPhiRes = normAdjPhiRes;
427  initNormAdjNuTildaRes = normAdjNuTildaRes;
428  }
429 
430  // Normalize residual norms:
431  normAdjURes[0] /= initNormAdjURes[0];
432  normAdjURes[1] /= initNormAdjURes[1];
433  normAdjURes[2] /= initNormAdjURes[2];
434  normAdjPRes /= initNormAdjPRes;
435  normAdjPhiRes /= initNormAdjPhiRes;
436  normAdjNuTildaRes /= initNormAdjNuTildaRes;
437 
438  if (cnt % printInterval == 0)
439  {
440  Info << "Normalized L2 norm of adjURes: " << normAdjURes[0] << " " << normAdjURes[1] << " " << normAdjURes[2] << endl;
441  Info << "Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
442  Info << "Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
443  Info << "Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
444  }
445 
446  // Check if fpRelTol is met:
447  if (normAdjURes[0] < fpRelTol && normAdjURes[1] < fpRelTol && normAdjURes[2] < fpRelTol && normAdjPRes < fpRelTol && normAdjPhiRes < fpRelTol && normAdjNuTildaRes < fpRelTol)
448  {
449  Info << "Residual drop of " << fpRelTol << " has been achieved after " << cnt << " steps! Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
450  Info << "Adjoint is converged and succussful.." << endl;
451  Info << "Normalized L2 norm of adjURes: " << normAdjURes[0] << " " << normAdjURes[1] << " " << normAdjURes[2] << endl;
452  Info << "Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
453  Info << "Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
454  Info << "Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
455  adjConv = 0;
456  break;
457  }
458  }
459 
460  // Update phiPsi
461  phiPsi += adjPhiRes;
462 
463  // ************************************************************************* //
464  // Step-1.2: Get D^(-T) adjURes, overwrite adjURes
465 
466  forAll(adjURes, cellI)
467  {
468  adjURes[cellI] /= UDiag[cellI];
469  }
470 
471  // ************************************************************************* //
472  // Step-2: Get alphaP.(F_p^grad)^T.(D^(-T).adjURes) - alphaP.adjPRes + (F_p^flux.F_p^aug)^T.adjPhiRes
473  // ************************************************************************* //
474  // Step-2.1: Start with (F_p^grad)^T.(D^(-T).adjURes) - adjPRes
475 
476  // Initialize step_2 as -adjPRes
477  //volScalarField step_2("step_2", -adjPRes);
478  // Use alias to save some mem usage
479  volScalarField& step_2 = adjPRes;
480  step_2 = -step_2;
481 
482  // Do (F_p^grad)^T.(D^(-T).adjURes) with ADR
483  // Here, F_p^grad is the operator of fvc::grad(p) * mesh().V() on p
484 
485  if (cnt == 0)
486  {
487  this->globalADTape_.setActive();
488 
489  // Now register p1 as input:
490  forAll(p1, cellI)
491  {
492  this->globalADTape_.registerInput(p1[cellI]);
493  }
494 
495  // Get tape position before calling function(s)
496  gradPStart_ = this->globalADTape_.getPosition();
497 
498  // Correct boundaries to link the intermediate results
499  p1.correctBoundaryConditions();
500  U.correctBoundaryConditions();
501 
502  // Calculate fvc::grad(p) * mesh().V()
503  gradP = fvc::grad(p1);
504  forAll(gradP, cellI)
505  {
506  gradP.primitiveFieldRef()[cellI] *= U.mesh().V()[cellI];
507  }
508 
509  // Get tape position before calling function(s)
510  gradPEnd_ = this->globalADTape_.getPosition();
511 
512  // stop recording
513  this->globalADTape_.setPassive();
514  }
515 
516  // set the AD seed to the output var
517  // Seed with (D^(-T).adjURes), which is the overwritten adjURes.
518  forAll(gradP, cellI)
519  {
520  gradP[cellI][0].setGradient(adjURes[cellI][0].getValue());
521  gradP[cellI][1].setGradient(adjURes[cellI][1].getValue());
522  gradP[cellI][2].setGradient(adjURes[cellI][2].getValue());
523  }
524 
525  // evaluate the tape to compute the derivative of the seeded output wrt all the input
526  this->globalADTape_.evaluate(gradPEnd_, gradPStart_);
527  forAll(p1, cellI)
528  {
529  step_2[cellI] += p1[cellI].getGradient();
530  }
531 
532  // Clear adjoints for future Jacobian calculations
533  this->globalADTape_.clearAdjoints();
534 
535  // ************************************************************************* //
536  // Step-2.2: Get alphaP*step_2.1 = alphaP.((F_p^grad)^T.(D^(-T).adjURes) - adjPRes)
537  scalar alphaP = p.mesh().fieldRelaxationFactor(p.name());
538  step_2 *= alphaP;
539 
540  //Info<<"Show step_2 after *alphaP: "<< step_2 << endl;
541 
542  // ************************************************************************* //
543  // Step-2.3: Get step_2.2 + (F_p^flux.F_p^aug)^T.adjPhiRes = alphaP.((F_p^grad)^T.(D(-T).adjURes) - adjPRes) + (F_p^flux.F_p^aug)^T.adjPhiRes
544  // Here, (F_p^flux.F_p^aug) is the operator pEqn.flux() on p, (F_p^flux.F_p^aug)^T.adjPhiRes is calculated with ADR
545 
546  if (cnt == 0)
547  {
548  // First, get pEqn1, note that we only need the l.h.s. of pEqn1
549  fvScalarMatrix pEqn1(
550  fvm::laplacian(rAU, p2));
551  //Info<<"Show pEqn1.flux(): "<< pEqn1.flux() << endl;
552 
553  // Re-active the tape, no need to initiate it again
554  this->globalADTape_.setActive();
555 
556  // Now register p2 as input:
557  forAll(p2, cellI)
558  {
559  this->globalADTape_.registerInput(p2[cellI]);
560  }
561 
562  // Get tape position before calling function(s)
563  pEqnFluxStart_ = this->globalADTape_.getPosition();
564 
565  // Correct boundaries to link the intermediate results
566  p2.correctBoundaryConditions();
567  U.correctBoundaryConditions();
568 
569  // Calculate pEqn1.flux():
570  pEqnFlux = pEqn1.flux();
571 
572  // Get tape position after calling function(s)
573  pEqnFluxEnd_ = this->globalADTape_.getPosition();
574 
575  // stop recording
576  this->globalADTape_.setPassive();
577  }
578 
579  // set the AD seed to the output var pEqnFlux
580  // Seed with adjPhiRes
581  // Seed internal of pEqnFlux:
582  forAll(pEqnFlux.primitiveFieldRef(), faceI)
583  {
584  pEqnFlux.primitiveFieldRef()[faceI].setGradient(adjPhiRes.primitiveFieldRef()[faceI].getValue());
585  }
586  // Seed boundary of pEqnFlux:
587  forAll(pEqnFlux.boundaryFieldRef(), patchI)
588  {
589  forAll(pEqnFlux.boundaryFieldRef()[patchI], faceI)
590  {
591  pEqnFlux.boundaryFieldRef()[patchI][faceI].setGradient(adjPhiRes.boundaryFieldRef()[patchI][faceI].getValue());
592  }
593  }
594 
595  // evaluate the tape to compute the derivative of the seeded output wrt all the input
596  this->globalADTape_.evaluate(pEqnFluxEnd_, pEqnFluxStart_);
597  forAll(p2, cellI)
598  {
599  step_2[cellI] += p2[cellI].getGradient();
600  }
601 
602  // Clear adjoints for future Jacobian calculations
603  this->globalADTape_.clearAdjoints();
604  // reset the tape for future recording
605  //tape.reset();
606 
607  //Info<<"Show final form of step_2: "<< step_2 << endl;
608 
609  // ************************************************************************* //
610  // Step-3: Solve and get step_3 = A_p^(-T).step_2,
611  // which is A_p^(-T).(alphaP.((F_p^grad)^T.(D(-T).adjURes) - adjPRes) + (F_p^flux.F_p^aug)^T.adjPhiRes)
612  // then update pPsi += step_3
613 
614  // Solve, using the inverse transpose product function
615  //invTranProd_pEqn(rAU, p, step_2, pseudoP);
616  //Info<<"Show pseudoP: "<< pseudoP << endl;
617 
618  // *********************************** //
619  // Below is equivalent to: invTranProd_pEqn(rAU, p, step_2, pseudoP);
620  // Here, pSource is a redundant alias, which is an "input argument" that takes the value of step_2 into solvePseudopEqn.H
621  volScalarField& pSource = step_2;
622 #include "solvePseudoPEqn.H"
623  // *********************************** //
624 
625  // Get step_3 as the alias of pseudoP
626  volScalarField& step_3 = pseudoP;
627  //Info<<"Show step_3: "<< step_3 << endl;
628 
629  // Update pPsi
630  pPsi += step_3;
631  //Info<<"Show pPsi: "<< pPsi << endl;
632 
633  // ************************************************************************* //
634  // Step-4: Get to the point before solve with A_U^(-T)
635  // which is Q^T.{D^(-T).adjURes - D^(-T)[E^T.F^T.M_p^(-T).alphaP.B^T.D^(-T).adjURes - E^T.F^T.M_p^(-T).alphaP.adjPRes + (E^T.F^T.M_p^(-T).G^T - E^T).adjPhiRes]}
636  // ************************************************************************* //
637  // Step-4.1: Get (F_phi^div)^T.step_3 - adjPhiRes = (F_phi^div)^T.(A_p^(-T).(alphaP.((F_p^grad)^T.(D(-T).adjURes) - adjPRes) + (F_p^flux.F_p^aug)^T.adjPhiRes)) - adjPhiRes
638  // (F_phi^div)^T.step_3 is calculated with ADR
639  // Here, (F_phi^div) is the operator fvc::div(phiHbyA) on phiHbyA
640  // We don't want to calculate phiHbyA, so let's use phi instead
641  // So, (F_phi^div) is the operator fvc::div(phi) on phi
642 
643  // Initiate a phi-like temp variable as - adjPhiRes
644  //surfaceScalarField phiTemp = -adjPhiRes;
645  // Use alia to same some mem usage
646  surfaceScalarField& phiTemp = adjPhiRes;
647  phiTemp = -phiTemp;
648 
649  if (cnt == 0)
650  {
651  // Re-active the tape, no need to initiate it again
652  this->globalADTape_.setActive();
653 
654  // Register phi1 as input
655  // Note that both the internal field of phi1 and the boundary fields need to be registered
656  // Register internal of phi1:
657  forAll(phi1.primitiveFieldRef(), faceI)
658  {
659  this->globalADTape_.registerInput(phi1.primitiveFieldRef()[faceI]);
660  }
661  // Register boundary of phi:
662  forAll(phi1.boundaryFieldRef(), patchI)
663  {
664  forAll(phi1.boundaryFieldRef()[patchI], faceI)
665  {
666  this->globalADTape_.registerInput(phi1.boundaryFieldRef()[patchI][faceI]);
667  }
668  }
669 
670  // Get tape position before calling function(s)
671  divPhiStart_ = this->globalADTape_.getPosition();
672 
673  // This is not needed, but whatever... XD
674  // Correct boundaries to link the intermediate results
675  //U.correctBoundaryConditions();
676  //p.correctBoundaryConditions();
677 
678  // Calculate fvc::div(phi), then scale with cell volumes mesh().V()
679  divPhi = fvc::div(phi1);
680  forAll(divPhi, cellI)
681  {
682  divPhi.primitiveFieldRef()[cellI] *= p.mesh().V()[cellI];
683  }
684 
685  // Get tape position after calling function(s)
686  divPhiEnd_ = this->globalADTape_.getPosition();
687 
688  // stop recording
689  this->globalADTape_.setPassive();
690  }
691 
692  // set the AD seed to the output var divPhi
693  // Seed with step_3
694  forAll(divPhi, cellI)
695  {
696  divPhi[cellI].setGradient(step_3[cellI].getValue());
697  }
698 
699  // evaluate the tape to compute the derivative of the seeded output wrt all the input
700  this->globalADTape_.evaluate(divPhiEnd_, divPhiStart_);
701  forAll(phi1.primitiveFieldRef(), faceI)
702  {
703  phiTemp.primitiveFieldRef()[faceI] += phi1.primitiveFieldRef()[faceI].getGradient();
704  }
705  forAll(phi1.boundaryFieldRef(), patchI)
706  {
707  forAll(phi1.boundaryFieldRef()[patchI], faceI)
708  {
709  phiTemp.boundaryFieldRef()[patchI][faceI] += phi1.boundaryFieldRef()[patchI][faceI].getGradient();
710  }
711  }
712 
713  // Clear adjoints for future Jacobian calculations
714  this->globalADTape_.clearAdjoints();
715  // reset the tape for future recording
716  //tape.reset();
717 
718  //Info<<"Show step-4.1: "<< phiTemp << endl;
719 
720  // ************************************************************************* //
721  // Step-4.2: Get (D^(-T).adjURes) - D^(-T).(F_U^flux.F_U^aug)^T.phiTemp
722  // Here, (F_U^flux.F_U^aug) is the operator fvc::flux(HbyA) on HbyA
723  // We don't want to calculate HbyA, so we use U instead
724  // HbyA and U have the same boundary settings, so the operator as a function stays the same
725  // So, (F_U^flux.F_U^aug) is the operator fvc::flux(U) on U
726 
727  // Update UPsi, for the 1st time
728  // We have to do this before creating step_4 as the alias of adjURes
729  forAll(UPsi, cellI)
730  {
731  UPsi[cellI] -= adjURes[cellI]; // Note that the overwritten adjURes is actually D^(-T).adjURes
732  }
733 
734  // Initiate a U-like temporary variable step_4 as D^(-T).adjURes, which is the overwritten adjURes
735  //volVectorField step_4("step_4", adjURes);
736  // Use alias to save some mem usage
737  volVectorField& step_4 = adjURes;
738 
739  if (cnt == 0)
740  {
741  // Re-active the tape, no need to initiate it again
742  this->globalADTape_.setActive();
743 
744  // Register U1 as input, note that U1 has 3 components
745  forAll(U1, cellI)
746  {
747  this->globalADTape_.registerInput(U1[cellI][0]);
748  this->globalADTape_.registerInput(U1[cellI][1]);
749  this->globalADTape_.registerInput(U1[cellI][2]);
750  }
751 
752  // Get tape position before calling function(s)
753  fluxUStart_ = this->globalADTape_.getPosition();
754 
755  // Correct boundaries to link the intermediate results
756  U1.correctBoundaryConditions();
757  p.correctBoundaryConditions();
758 
759  // Calculate fvc::flux(U)
760  fluxU = fvc::flux(U1);
761 
762  // Get tape position after calling function(s)
763  fluxUEnd_ = this->globalADTape_.getPosition();
764 
765  // stop recording
766  this->globalADTape_.setPassive();
767  }
768 
769  // set the AD seed to the output var fluxU
770  // Seed with phiTemp
771  // Seed internal of fluxU:
772  forAll(fluxU.primitiveFieldRef(), faceI)
773  {
774  fluxU.primitiveFieldRef()[faceI].setGradient(phiTemp.primitiveFieldRef()[faceI].getValue());
775  }
776  // Seed boundary of fluxU:
777  forAll(fluxU.boundaryFieldRef(), patchI)
778  {
779  forAll(fluxU.boundaryFieldRef()[patchI], faceI)
780  {
781  fluxU.boundaryFieldRef()[patchI][faceI].setGradient(phiTemp.boundaryFieldRef()[patchI][faceI].getValue());
782  }
783  }
784 
785  // evaluate the tape to compute the derivative of the seeded output wrt all the input
786  this->globalADTape_.evaluate(fluxUEnd_, fluxUStart_);
787  forAll(U1, cellI)
788  {
789  for (label cmpt = 0; cmpt < 3; cmpt++)
790  {
791  step_4[cellI][cmpt] -= U1[cellI][cmpt].getGradient() / UDiag[cellI];
792  }
793  }
794 
795  // Clear adjoints for future Jacobian calculations
796  this->globalADTape_.clearAdjoints();
797  // reset the tape for future recording
798  //tape.reset();
799 
800  //Info<<"Show step_4: "<< step_4 << endl;
801 
802  // ************************************************************************* //
803  // Step-4.3: Get Q^T.step_4, this is the last sub-step of step-4
804  // We use pseudoUEqn, whose upper and lower are already swapped
805  step_4.primitiveFieldRef() = -pseudoUEqn.lduMatrix::H(step_4);
806  //Info<<"Show step_4: "<< step_4 << endl;
807 
808  // ************************************************************************* //
809  // Step-5: Solve and get step_5 = A_U^(-T).step_4
810  // Then update UPsi += (step_5 - D^(-T).adjURes)
811 
812  // Solve, using the inverse transpose product function
813  //invTranProd_UEqn(U, phi, nuEff, p, step_4, pseudoU);
814  //Info<<"Show pseudoU: "<< pseudoU << endl;
815 
816  // *********************************** //
817  // Below is equivalent to: invTranProd_UEqn(U, phi, nuEff, p, step_4, pseudoU);
818  // Here, USource is a redundant alias, which is an "input argument" that takes the value of step_4 into solvePseudoUEqn.H
819  volVectorField& USource = step_4;
820 #include "solvePseudoUEqn.H"
821  // *********************************** //
822 
823  // Get step_5 as the alias of pseudoU
824  volVectorField& step_5 = pseudoU;
825  //Info<<"Show step_5: "<< step_5 << endl;
826 
827  // Update UPsi, for the 2nd time
828  forAll(UPsi, cellI)
829  {
830  UPsi[cellI] += step_5[cellI];
831  }
832  //Info<<"Show UPsi: "<< UPsi << endl;
833 
834  // ************************************************************************* //
835  // Step-6: (SA turbulence) Solve and get pseudoNuTilda = A_nuTilda^(-T).adjNuTildaRes
836  // Then update nuTildaPsi -= pseudoNuTilda
837 
838  // Solve, using the inverse transpose product function
839  //invTranProd_nuTildaEqn(U, phi, nuTilda, y, nu, sigmaNut, kappa, Cb1, Cb2, Cw1, Cw2, Cw3, Cv1, Cv2, adjNuTildaRes, pseudoNuTilda);
840 
841  // *********************************** //
842  // Below is equivalent to: invTranProd_nuTildaEqn(U, phi, nuTilda, y, nu, sigmaNut, kappa, Cb1, Cb2, Cw1, Cw2, Cw3, Cv1, Cv2, adjNuTildaRes, pseudoNuTilda);
843  // Here, nuTildaSource is a redundant alias, which is an "input argument" that takes the value of adjNuTildaRes into solvePseudoUEqn.H
844 
845  volScalarField& nuTildaSource = adjNuTildaRes;
846 
847  daTurbulenceModelPtr_->rhsSolvePseudoNuTildaEqn(nuTildaSource);
848  //#include "solvePseudonuTildaEqn.H"
849  // *********************************** //
850 
851  // Update nuTildaPsi
852  forAll(nuTildaPsi, cellI)
853  {
854  nuTildaPsi[cellI] -= pseudoNuTilda[cellI];
855  }
856 
857  cnt++;
858  }
859 
860  // If fpRelTol not met, check if the relaxed adjoint relTol is met
861  if (adjConv == 1)
862  {
863  Info << "Prescribed adjoint convergence of " << fpRelTol << " has not been met after " << cnt << " steps! Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
864  Info << "Now checking if the relaxed convergence is met..." << endl;
865  this->calcAdjointResidual(URes, pRes, phiRes, nuTildaRes, dFdW, UPsi, pPsi, phiPsi, nuTildaPsi, adjURes, adjPRes, adjPhiRes, adjNuTildaRes, cnt);
866 
867  vector normAdjURes = L2norm(adjURes);
868  scalar normAdjPRes = L2norm(adjPRes);
869  scalar normAdjPhiRes = L2norm(adjPhiRes);
870  scalar normAdjNuTildaRes = L2norm(adjNuTildaRes);
871 
872  // Normalize residual norms:
873  normAdjURes[0] /= initNormAdjURes[0];
874  normAdjURes[1] /= initNormAdjURes[1];
875  normAdjURes[2] /= initNormAdjURes[2];
876  normAdjPRes /= initNormAdjPRes;
877  normAdjPhiRes /= initNormAdjPhiRes;
878  normAdjNuTildaRes /= initNormAdjNuTildaRes;
879 
880  scalar relaxedFpRelTol = fpRelTol * fpMinResTolDiff;
881  if (normAdjURes[0] < relaxedFpRelTol && normAdjURes[1] < relaxedFpRelTol && normAdjURes[2] < relaxedFpRelTol && normAdjPRes < relaxedFpRelTol && normAdjPhiRes < relaxedFpRelTol && normAdjNuTildaRes < relaxedFpRelTol)
882  {
883  Info << "Relaxed residual drop of " << relaxedFpRelTol << " has been achieved! Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
884  Info << "Adjoint is still considered succussful.." << endl;
885  adjConv = 0;
886  }
887  else
888  {
889  Info << "Relaxed residual drop of " << relaxedFpRelTol << " has not been met! Execution Time: " << meshPtr_->time().elapsedCpuTime() << " s" << endl;
890  Info << "Adjoint has failed to converge.." << endl;
891  }
892  Info << "Normalized L2 norm of adjURes: " << normAdjURes[0] << " " << normAdjURes[1] << " " << normAdjURes[2] << endl;
893  Info << "Normalized L2 norm of adjPRes: " << normAdjPRes << endl;
894  Info << "Normalized L2 norm of adjPhiRes: " << normAdjPhiRes << endl;
895  Info << "Normalized L2 norm of adjNuTildaRes: " << normAdjNuTildaRes << endl;
896  }
897 
898  // converged, assign the field Psi to psiVec
899  this->vec2Fields("field2Vec", psi, UPsi, pPsi, phiPsi, nuTildaPsi);
900  }
901  else
902  {
903  FatalErrorIn("adjEqnSolMethod not valid") << exit(FatalError);
904  }
905 
906 #endif
907  return adjConv;
908 }
909 
911  const word mode,
912  Vec cVec,
913  volVectorField& UField,
914  volScalarField& pField,
915  surfaceScalarField& phiField,
916  volScalarField& nuTildaField)
917 {
918 #ifdef CODI_ADR
919  PetscScalar* cVecArray;
920  if (mode == "vec2Field")
921  {
922  VecGetArray(cVec, &cVecArray);
923 
924  // U
925  forAll(meshPtr_->cells(), cellI)
926  {
927  for (label comp = 0; comp < 3; comp++)
928  {
929  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("U", cellI, comp);
930  UField[cellI][comp] = cVecArray[adjLocalIdx];
931  }
932  }
933  // p
934  forAll(meshPtr_->cells(), cellI)
935  {
936  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("p", cellI);
937  pField[cellI] = cVecArray[adjLocalIdx];
938  }
939  // phi
940  forAll(meshPtr_->faces(), faceI)
941  {
942  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("phi", faceI);
943 
944  if (faceI < daIndexPtr_->nLocalInternalFaces)
945  {
946  phiField[faceI] = cVecArray[adjLocalIdx];
947  }
948  else
949  {
950  label relIdx = faceI - daIndexPtr_->nLocalInternalFaces;
951  label patchIdx = daIndexPtr_->bFacePatchI[relIdx];
952  label faceIdx = daIndexPtr_->bFaceFaceI[relIdx];
953  phiField.boundaryFieldRef()[patchIdx][faceIdx] = cVecArray[adjLocalIdx];
954  }
955  }
956  // nuTilda
957  forAll(meshPtr_->cells(), cellI)
958  {
959  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("nuTilda", cellI);
960  nuTildaField[cellI] = cVecArray[adjLocalIdx];
961  }
962 
963  VecRestoreArray(cVec, &cVecArray);
964  }
965  else if (mode == "field2Vec")
966  {
967  VecGetArray(cVec, &cVecArray);
968 
969  // U
970  forAll(meshPtr_->cells(), cellI)
971  {
972  for (label comp = 0; comp < 3; comp++)
973  {
974  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("U", cellI, comp);
975  cVecArray[adjLocalIdx] = UField[cellI][comp].value();
976  }
977  }
978  // p
979  forAll(meshPtr_->cells(), cellI)
980  {
981  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("p", cellI);
982  cVecArray[adjLocalIdx] = pField[cellI].value();
983  }
984  // phi
985  forAll(meshPtr_->faces(), faceI)
986  {
987  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("phi", faceI);
988 
989  if (faceI < daIndexPtr_->nLocalInternalFaces)
990  {
991  cVecArray[adjLocalIdx] = phiField[faceI].value();
992  }
993  else
994  {
995  label relIdx = faceI - daIndexPtr_->nLocalInternalFaces;
996  label patchIdx = daIndexPtr_->bFacePatchI[relIdx];
997  label faceIdx = daIndexPtr_->bFaceFaceI[relIdx];
998  cVecArray[adjLocalIdx] = phiField.boundaryFieldRef()[patchIdx][faceIdx].value();
999  }
1000  }
1001  // nuTilda
1002  forAll(meshPtr_->cells(), cellI)
1003  {
1004  label adjLocalIdx = daIndexPtr_->getLocalAdjointStateIndex("nuTilda", cellI);
1005  cVecArray[adjLocalIdx] = nuTildaField[cellI].value();
1006  }
1007 
1008  VecRestoreArray(cVec, &cVecArray);
1009  }
1010  else
1011  {
1012  FatalErrorIn("mode not valid") << exit(FatalError);
1013  }
1014 #endif
1015 }
1016 
1018  const volVectorField& mySource,
1019  volVectorField& pseudoU)
1020 {
1021  /*
1022  Description:
1023  Inverse transpose product, MU^(-T)
1024  Based on inverseProduct_UEqn from simpleFoamPrimal, but swaping upper() and lower()
1025  We won't ADR this function, so we can treat most of the arguments as const
1026  */
1027 
1028  /*
1029  const objectRegistry& db = meshPtr_->thisDb();
1030  const surfaceScalarField& phi = db.lookupObject<surfaceScalarField>("phi");
1031  volScalarField nuEff = daTurbulenceModelPtr_->nuEff();
1032 
1033  // Get the pseudoUEqn,
1034  // the most important thing here is to make sure the l.h.s. matches that of UEqn.
1035  fvVectorMatrix pseudoUEqn(
1036  fvm::div(phi, pseudoU, "div(phi,U)")
1037  - fvm::laplacian(nuEff, pseudoU)
1038  - fvc::div(nuEff * dev2(T(fvc::grad(pseudoU))), "div((nuEff*dev2(T(grad(U)))))"));
1039  pseudoUEqn.relax(relaxUEqn_);
1040 
1041  // Swap upper() and lower()
1042  List<scalar> temp = pseudoUEqn.upper();
1043  pseudoUEqn.upper() = pseudoUEqn.lower();
1044  pseudoUEqn.lower() = temp;
1045 
1046  // Overwrite the r.h.s.
1047  pseudoUEqn.source() = mySource;
1048 
1049  // Make sure that boundary contribution to source is zero,
1050  // Alternatively, we can deduct source by boundary contribution, so that it would cancel out during solve.
1051  forAll(pseudoU.boundaryField(), patchI)
1052  {
1053  const fvPatch& pp = pseudoU.boundaryField()[patchI].patch();
1054  forAll(pp, faceI)
1055  {
1056  label cellI = pp.faceCells()[faceI];
1057  pseudoUEqn.source()[cellI] -= pseudoUEqn.boundaryCoeffs()[patchI][faceI];
1058  }
1059  }
1060 
1061  // Before solve, force xEqn.psi() to be solved into all zero
1062  forAll(pseudoU.primitiveFieldRef(), cellI)
1063  {
1064  pseudoU.primitiveFieldRef()[cellI][0] = 0;
1065  pseudoU.primitiveFieldRef()[cellI][1] = 0;
1066  pseudoU.primitiveFieldRef()[cellI][2] = 0;
1067  }
1068 
1069  pseudoUEqn.solve(solverDictU_);
1070 */
1071 }
1072 
1074  const volScalarField& mySource,
1075  volScalarField& pseudoP)
1076 {
1077  /*
1078  Description:
1079  Inverse transpose product, Mp^(-T)
1080  Based on inverseProduct_pEqn from simpleFoamPrimal, but swaping upper() and lower()
1081  We won't ADR this function, so we can treat most of the arguments as const
1082  */
1083 
1084  /*
1085  const objectRegistry& db = meshPtr_->thisDb();
1086  const volVectorField& U = db.lookupObject<volVectorField>("U");
1087  const surfaceScalarField& phi = db.lookupObject<surfaceScalarField>("phi");
1088  volScalarField nuEff = daTurbulenceModelPtr_->nuEff();
1089 
1090  // Construct UEqn first
1091  fvVectorMatrix UEqn(
1092  fvm::div(phi, U)
1093  - fvm::laplacian(nuEff, U)
1094  - fvc::div(nuEff * dev2(T(fvc::grad(U)))));
1095  // Without this, pRes would be way off.
1096  UEqn.relax();
1097 
1098  // create a scalar field with 1/A, reverse of A() of U
1099  volScalarField rAU(1.0 / UEqn.A());
1100 
1101  // Get the pseudoPEqn,
1102  // the most important thing here is to make sure the l.h.s. matches that of pEqn.
1103  fvScalarMatrix pseudoPEqn(fvm::laplacian(rAU, pseudoP));
1104 
1105  // Swap upper() and lower()
1106  List<scalar> temp = pseudoPEqn.upper();
1107  pseudoPEqn.upper() = pseudoPEqn.lower();
1108  pseudoPEqn.lower() = temp;
1109 
1110  // Overwrite the r.h.s.
1111  pseudoPEqn.source() = mySource;
1112 
1113  // pEqn.setReference(pRefCell, pRefValue);
1114  // Here, pRefCell is a label, and pRefValue is a scalar
1115  // In actual implementation, they need to passed into this function.
1116  pseudoPEqn.setReference(0, 0.0);
1117 
1118  // Make sure that boundary contribution to source is zero,
1119  // Alternatively, we can deduct source by boundary contribution, so that it would cancel out during solve.
1120  forAll(pseudoP.boundaryField(), patchI)
1121  {
1122  const fvPatch& pp = pseudoP.boundaryField()[patchI].patch();
1123  forAll(pp, faceI)
1124  {
1125  label cellI = pp.faceCells()[faceI];
1126  pseudoPEqn.source()[cellI] -= pseudoPEqn.boundaryCoeffs()[patchI][faceI];
1127  }
1128  }
1129 
1130  // Before solve, force xEqn.psi() to be solved into all zero
1131  forAll(pseudoP.primitiveFieldRef(), cellI)
1132  {
1133  pseudoP.primitiveFieldRef()[cellI] = 0;
1134  }
1135 
1136  pseudoPEqn.solve(solverDictP_);
1137 */
1138 }
1139 
1141  volVectorField& URes,
1142  volScalarField& pRes,
1143  surfaceScalarField& phiRes)
1144 {
1145  const objectRegistry& db = meshPtr_->thisDb();
1146  const volVectorField& U = db.lookupObject<volVectorField>("U");
1147  const volScalarField& p = db.lookupObject<volScalarField>("p");
1148  const surfaceScalarField& phi = db.lookupObject<surfaceScalarField>("phi");
1149  volScalarField nuEff = daTurbulenceModelPtr_->nuEff();
1150 
1151  fvVectorMatrix UEqn(
1152  fvm::div(phi, U)
1153  - fvm::laplacian(nuEff, U)
1154  - fvc::div(nuEff * dev2(T(fvc::grad(U))))); //This term is needed in res though...
1155 
1156  List<vector>& USource = UEqn.source();
1157  // Note we cannot use UEqn.D() here, because boundary contribution to diag have 3 components, and they can be different.
1158  // Thus we use UEqn.diag() here, and we correct both source and diag later.
1159  List<scalar>& UDiag = UEqn.diag();
1160 
1161  // Get fvc::grad(p), so that it can be added to r.h.s.
1162  volVectorField gradP(fvc::grad(p));
1163 
1164  // Initiate URes, with no boundary contribution
1165  for (label i = 0; i < U.size(); i++)
1166  {
1167  URes[i] = UDiag[i] * U[i] - USource[i] + U.mesh().V()[i] * gradP[i];
1168  }
1169  URes.primitiveFieldRef() -= UEqn.lduMatrix::H(U);
1170 
1171  // Add boundary contribution to source and diag
1172  forAll(U.boundaryField(), patchI)
1173  {
1174  const fvPatch& pp = U.boundaryField()[patchI].patch();
1175  forAll(pp, faceI)
1176  {
1177  // Both ways of getting cellI work
1178  // Below is the previous way of getting the address
1179  label cellI = pp.faceCells()[faceI];
1180  // Below is using lduAddr().patchAddr(patchi)
1181  //label cellI = UEqn.lduAddr().patchAddr(patchI)[faceI];
1182  for (label cmpt = 0; cmpt < 3; cmpt++)
1183  {
1184  URes[cellI][cmpt] += UEqn.internalCoeffs()[patchI][faceI][cmpt] * U[cellI][cmpt];
1185  }
1186  //Info << "UEqn.internalCoeffs()[" << patchI << "][" << faceI <<"]= " << UEqn.internalCoeffs()[patchI][faceI] <<endl;
1187  URes[cellI] -= UEqn.boundaryCoeffs()[patchI][faceI];
1188  }
1189  }
1190 
1191  // Below is not necessary, but it doesn't hurt
1192  URes.correctBoundaryConditions();
1193 
1194  UEqn.relax(); // Without this, pRes would be way off.
1195 
1196  volScalarField rAU(1.0 / UEqn.A()); // create a scalar field with 1/A, reverse of A() of U
1197  volVectorField HbyA("HbyA", U); // initialize a vector field with U and pass it to HbyA
1198  HbyA = rAU * UEqn.H(); // basically, HbyA = 1/A * H, H_by_A, need to verify source code though...
1199  surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); // get the flux of HbyA, phi_H_by_A
1200 
1201  fvScalarMatrix pEqn(
1202  fvm::laplacian(rAU, p) == fvc::div(phiHbyA));
1203 
1204  List<scalar>& pSource = pEqn.source();
1205  List<scalar>& pDiag = pEqn.diag();
1206 
1207  // Initiate pRes, with no boundary contribution
1208  for (label i = 0; i < p.size(); i++)
1209  {
1210  pRes[i] = pDiag[i] * p[i] - pSource[i];
1211  }
1212  pRes.primitiveFieldRef() -= pEqn.lduMatrix::H(p);
1213 
1214  // Boundary correction
1215  forAll(p.boundaryField(), patchI)
1216  {
1217  const fvPatch& pp = p.boundaryField()[patchI].patch();
1218  forAll(pp, faceI)
1219  {
1220  // Both ways of getting cellI work
1221  // Below is the previous way of getting the address
1222  label cellI = pp.faceCells()[faceI];
1223  // Below is using lduAddr().patchAddr(patchi)
1224  //label cellI = pEqn.lduAddr().patchAddr(patchI)[faceI];
1225  //myDiag[cellI] += TEqn.internalCoeffs()[patchI][faceI];
1226  pRes[cellI] += pEqn.internalCoeffs()[patchI][faceI] * p[cellI];
1227  pRes[cellI] -= pEqn.boundaryCoeffs()[patchI][faceI];
1228  }
1229  }
1230 
1231  // Below is not necessary, but it doesn't hurt
1232  pRes.correctBoundaryConditions();
1233 
1234  // Then do phiRes
1235  // Note: DAFoam also uses this formula for phiRes
1236  phiRes = phiHbyA - pEqn.flux() - phi;
1237 }
1238 
1240  volVectorField& URes,
1241  volScalarField& pRes,
1242  surfaceScalarField& phiRes,
1243  volScalarField& nuTildaRes,
1244  Vec dFdW,
1245  volVectorField& UPsi,
1246  volScalarField& pPsi,
1247  surfaceScalarField& phiPsi,
1248  volScalarField& nuTildaPsi,
1249  volVectorField& adjURes,
1250  volScalarField& adjPRes,
1251  surfaceScalarField& adjPhiRes,
1252  volScalarField& adjNuTildaRes,
1253  label& cnt)
1254 {
1255 #ifdef CODI_ADR
1256  volVectorField& U = const_cast<volVectorField&>(meshPtr_->thisDb().lookupObject<volVectorField>("U"));
1257  volScalarField& p = const_cast<volScalarField&>(meshPtr_->thisDb().lookupObject<volScalarField>("p"));
1258  volScalarField& nuTilda = const_cast<volScalarField&>(meshPtr_->thisDb().lookupObject<volScalarField>("nuTilda"));
1259  surfaceScalarField& phi = const_cast<surfaceScalarField&>(meshPtr_->thisDb().lookupObject<surfaceScalarField>("phi"));
1260 
1261  // Pass -dfdw to adjRes:
1262  this->vec2Fields("vec2Field", dFdW, adjURes, adjPRes, adjPhiRes, adjNuTildaRes);
1263  adjURes = -adjURes;
1264  adjPRes = -adjPRes;
1265  adjPhiRes = -adjPhiRes;
1266  adjNuTildaRes = -adjNuTildaRes;
1267 
1268  if (cnt == 0)
1269  {
1270  this->globalADTape_.reset();
1271  this->globalADTape_.setActive();
1272 
1273  // register all (3+1) state variables as input
1274  // Start with U, note that U has 3 components
1275  forAll(U, cellI)
1276  {
1277  this->globalADTape_.registerInput(U[cellI][0]);
1278  this->globalADTape_.registerInput(U[cellI][1]);
1279  this->globalADTape_.registerInput(U[cellI][2]);
1280  }
1281  // Now register p as input:
1282  forAll(p, cellI)
1283  {
1284  this->globalADTape_.registerInput(p[cellI]);
1285  }
1286  // Then, register phi as input
1287  // Note that both the internal field of phi and the boundary fields need to be registered
1288  // Register internal of phi:
1289  forAll(phi.primitiveFieldRef(), faceI)
1290  {
1291  this->globalADTape_.registerInput(phi.primitiveFieldRef()[faceI]);
1292  }
1293  // Register boundary of phi:
1294  forAll(phi.boundaryFieldRef(), patchI)
1295  {
1296  forAll(phi.boundaryFieldRef()[patchI], faceI)
1297  {
1298  this->globalADTape_.registerInput(phi.boundaryFieldRef()[patchI][faceI]);
1299  }
1300  }
1301  // And then, register turbulence variable nuTilda as input:
1302  forAll(nuTilda, cellI)
1303  {
1304  this->globalADTape_.registerInput(nuTilda[cellI]);
1305  }
1306 
1307  // Get tape position before calling function(s)
1308  adjResStart_ = this->globalADTape_.getPosition();
1309 
1310  // Correct boundaries to link the intermediate results
1311  U.correctBoundaryConditions();
1312  p.correctBoundaryConditions();
1313  nuTilda.correctBoundaryConditions();
1314 
1315  // Construct nuEff before calling lduCalcAllRes
1316  daTurbulenceModelPtr_->updateIntermediateVariables();
1317 
1318  // Call the residual functions
1319  this->calcLduResiduals(URes, pRes, phiRes);
1320  daTurbulenceModelPtr_->calcLduResidualTurb(nuTildaRes);
1321 
1322  // No registerOutput needed in positional tapes
1323  /*
1324  // register output
1325  forAll(URes, cellI)
1326  {
1327  tape.registerOutput(URes[cellI][0]);
1328  tape.registerOutput(URes[cellI][1]);
1329  tape.registerOutput(URes[cellI][2]);
1330  }
1331  forAll(pRes, cellI)
1332  {
1333  tape.registerOutput(pRes[cellI]);
1334  }
1335  forAll(phiRes.primitiveFieldRef(), faceI)
1336  {
1337  tape.registerOutput(phiRes[faceI]);
1338  }
1339  // Seed boundary of phiRes:
1340  forAll(phiRes.boundaryFieldRef(), patchI)
1341  {
1342  forAll(phiRes.boundaryFieldRef()[patchI], faceI)
1343  {
1344  tape.registerOutput(phiRes.boundaryFieldRef()[patchI][faceI]);
1345  }
1346  }
1347  forAll(nuTildaRes, cellI)
1348  {
1349  tape.registerOutput(nuTildaRes[cellI]);
1350  }
1351  */
1352 
1353  // Get tape position after calling function(s)
1354  adjResEnd_ = this->globalADTape_.getPosition();
1355 
1356  // stop recording
1357  this->globalADTape_.setPassive();
1358  }
1359 
1360  // set the AD seed to the output var
1361  // Start with URes, note that URes has 3 components
1362  forAll(URes, cellI)
1363  {
1364  URes[cellI][0].setGradient(UPsi[cellI][0].getValue());
1365  URes[cellI][1].setGradient(UPsi[cellI][1].getValue());
1366  URes[cellI][2].setGradient(UPsi[cellI][2].getValue());
1367  }
1368  // Now seed pRes:
1369  forAll(pRes, cellI)
1370  {
1371  pRes[cellI].setGradient(pPsi[cellI].getValue());
1372  }
1373  // Then, seed phiRes:
1374  // Seed internal of phiRes:
1375  forAll(phiRes.primitiveFieldRef(), faceI)
1376  {
1377  phiRes.primitiveFieldRef()[faceI].setGradient(phiPsi.primitiveFieldRef()[faceI].getValue());
1378  }
1379  // Seed boundary of phiRes:
1380  forAll(phiRes.boundaryFieldRef(), patchI)
1381  {
1382  forAll(phiRes.boundaryFieldRef()[patchI], faceI)
1383  {
1384  phiRes.boundaryFieldRef()[patchI][faceI].setGradient(phiPsi.boundaryFieldRef()[patchI][faceI].getValue());
1385  }
1386  }
1387  // And then, seed nuTildaRes:
1388  forAll(nuTildaRes, cellI)
1389  {
1390  nuTildaRes[cellI].setGradient(nuTildaPsi[cellI].getValue());
1391  }
1392 
1393  // evaluate the tape to compute the derivative of the seeded output wrt all the input
1394  this->globalADTape_.evaluate(adjResEnd_, adjResStart_);
1395  forAll(U, cellI)
1396  {
1397  adjURes[cellI][0] += U[cellI][0].getGradient();
1398  adjURes[cellI][1] += U[cellI][1].getGradient();
1399  adjURes[cellI][2] += U[cellI][2].getGradient();
1400  }
1401  forAll(p, cellI)
1402  {
1403  adjPRes[cellI] += p[cellI].getGradient();
1404  }
1405  forAll(phi.primitiveFieldRef(), faceI)
1406  {
1407  adjPhiRes.primitiveFieldRef()[faceI] += phi.primitiveFieldRef()[faceI].getGradient();
1408  }
1409  forAll(phi.boundaryFieldRef(), patchI)
1410  {
1411  forAll(phi.boundaryFieldRef()[patchI], faceI)
1412  {
1413  adjPhiRes.boundaryFieldRef()[patchI][faceI] += phi.boundaryFieldRef()[patchI][faceI].getGradient();
1414  }
1415  }
1416  forAll(nuTilda, cellI)
1417  {
1418  adjNuTildaRes[cellI] += nuTilda[cellI].getGradient();
1419  }
1420 
1421  // Clear adjoints for future Jacobian calculations
1422  this->globalADTape_.clearAdjoints();
1423 #endif
1424 }
1425 
1426 // L2 norm of a scalar list (or the primitive field of a volScalarField)
1427 // Scale off cell volumes
1428 scalar DASimpleFoam::L2norm(const volScalarField& v)
1429 {
1430  scalar L2normV = 0.0;
1431 
1432  forAll(v, cellI)
1433  {
1434  L2normV += sqr(v[cellI] / meshPtr_->V()[cellI]);
1435  }
1436  L2normV = sqrt(L2normV);
1437 
1438  return L2normV;
1439 }
1440 
1441 // Scale off cell volumes
1442 vector DASimpleFoam::L2norm(const volVectorField& U)
1443 {
1444  vector L2normU = vector::zero;
1445 
1446  forAll(U, cellI)
1447  {
1448  L2normU[0] += sqr(U[cellI][0] / meshPtr_->V()[cellI]);
1449  L2normU[1] += sqr(U[cellI][1] / meshPtr_->V()[cellI]);
1450  L2normU[2] += sqr(U[cellI][2] / meshPtr_->V()[cellI]);
1451  }
1452  L2normU[0] = sqrt(L2normU[0]);
1453  L2normU[1] = sqrt(L2normU[1]);
1454  L2normU[2] = sqrt(L2normU[2]);
1455 
1456  return L2normU;
1457 }
1458 
1459 // L2 norm of a surfaceScalarField
1460 scalar DASimpleFoam::L2norm(const surfaceScalarField& Phi)
1461 {
1462  scalar L2normPhi = 0.0;
1463 
1464  forAll(Phi.primitiveField(), faceI)
1465  {
1466  L2normPhi += sqr(Phi.primitiveField()[faceI]);
1467  }
1468  forAll(Phi.boundaryField(), patchI)
1469  {
1470  forAll(Phi.boundaryField()[patchI], faceI)
1471  {
1472  L2normPhi += sqr(Phi.boundaryField()[patchI][faceI]);
1473  }
1474  }
1475  L2normPhi = sqrt(L2normPhi);
1476 
1477  return L2normPhi;
1478 }
1479 
1480 // This function is for swapping the upper and lower of a xxEqn
1481 void DASimpleFoam::swap(List<scalar>& a, List<scalar>& b)
1482 {
1483  List<scalar> temp = a;
1484  a = b;
1485  b = temp;
1486 }
1487 
1488 } // End namespace Foam
1489 
1490 // ************************************************************************* //
Foam::DASimpleFoam::DASimpleFoam
DASimpleFoam(char *argsAll, PyObject *pyOptions)
Definition: DASimpleFoam.C:41
Foam::DASimpleFoam::relaxUEqn_
scalar relaxUEqn_
fvSolution parameters for fixed-point adjoint
Definition: DASimpleFoam.H:133
Foam::DASimpleFoam::swap
void swap(List< scalar > &a, List< scalar > &b)
Definition: DASimpleFoam.C:1481
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Foam::DASolver::checkPrimalFailure
label checkPrimalFailure()
check whether the primal fails based on residual and regression fail flag
Definition: DASolver.C:2472
Foam::DASolver::loop
label loop(Time &runTime)
return whether to loop the primal solution, similar to runTime::loop() except we don't do file IO
Definition: DASolver.C:147
Foam::DASolver::printToScreen_
label printToScreen_
whether to print primal information to the screen
Definition: DASolver.H:150
pseudoP
volScalarField pseudoP("pseudoP", p)
Foam::DASolver::daRegressionPtr_
autoPtr< DARegression > daRegressionPtr_
DARegression pointer.
Definition: DASolver.H:111
Foam::DASolver::daModelPtr_
autoPtr< DAModel > daModelPtr_
DAModel pointer.
Definition: DASolver.H:84
Foam::DASimpleFoam::invTranProdUEqn
void invTranProdUEqn(const volVectorField &mySource, volVectorField &pseudoU)
Inverse transpose product, MU^(-T)
Definition: DASimpleFoam.C:1017
Foam::DASimpleFoam::solvePrimal
virtual label solvePrimal()
solve the primal equations
Definition: DASimpleFoam.C:123
Foam::DASolver::daIndexPtr_
autoPtr< DAIndex > daIndexPtr_
DAIndex pointer.
Definition: DASolver.H:87
Foam::DASimpleFoam::solverDictU_
dictionary solverDictU_
Definition: DASimpleFoam.H:134
Foam::DASimpleFoam::daTurbulenceModelPtr_
autoPtr< DATurbulenceModel > daTurbulenceModelPtr_
DATurbulenceModel pointer.
Definition: DASimpleFoam.H:79
Foam::DASolver
Definition: DASolver.H:54
pseudoU
volVectorField pseudoU("pseudoU", U)
Foam::DASimpleFoam::hasTField_
label hasTField_
whether to have the temperature field
Definition: DASimpleFoam.H:115
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
UEqn
fvVectorMatrix & UEqn
Definition: UEqnPimpleDyM.H:15
pseudoNuTilda
volScalarField & pseudoNuTilda
Definition: pseudoEqns.H:34
Foam::DASolver::daOptionPtr_
autoPtr< DAOption > daOptionPtr_
DAOption pointer.
Definition: DASolver.H:81
DASimpleFoam.H
Foam::DATurbulenceModel::New
static autoPtr< DATurbulenceModel > New(const word modelType, const fvMesh &mesh, const DAOption &daOption)
Definition: DATurbulenceModel.C:161
Foam::DASimpleFoam::hasFvSource_
label hasFvSource_
whether to have fvSource term
Definition: DASimpleFoam.H:91
createAdjoint.H
Foam::DASimpleFoam::initSolver
virtual void initSolver()
initialize fields and variables
Definition: DASimpleFoam.C:81
HbyA
volVectorField & HbyA
Definition: pEqnPimpleDyM.H:17
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DASimpleFoam::solverDictP_
dictionary solverDictP_
Definition: DASimpleFoam.H:135
phi1
phi1
Definition: p1EqnIrkPimple.H:19
Foam::DASimpleFoam::daFvSourcePtr_
autoPtr< DAFvSource > daFvSourcePtr_
DASource pointer.
Definition: DASimpleFoam.H:82
phiHbyA
phiHbyA
Definition: pEqnRhoSimpleC.H:36
createFieldsSimple.H
p
volScalarField & p
Definition: createRefsPimpleDyM.H:6
Foam::DASimpleFoam::runFPAdj
virtual label runFPAdj(Vec dFdW, Vec psi)
solve the adjoint equation using the fixed-point iteration method
Definition: DASimpleFoam.C:189
createSimpleControlPython.H
Foam::DAFvSource::New
static autoPtr< DAFvSource > New(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAFvSource.C:47
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
phi
surfaceScalarField & phi
Definition: createRefsPimpleDyM.H:8
Foam::DASimpleFoam::L2norm
scalar L2norm(const volScalarField &v)
Definition: DASimpleFoam.C:1428
Foam
Definition: checkGeometry.C:32
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
Foam::DASolver::writeAssociatedFields
void writeAssociatedFields()
write associated fields such as relative velocity
Definition: DASolver.C:2526
TEqnSimple.H
UDiag
List< scalar > UDiag
Definition: pseudoEqns.H:18
allOptions
const dictionary & allOptions
Definition: createRefsPimpleDyM.H:14
Foam::DASolver::regModelFail_
label regModelFail_
whether the regModel compute fails
Definition: DASolver.H:153
pseudoEqns.H
laminarTransport
singlePhaseTransportModel & laminarTransport
Definition: createRefsPimpleDyM.H:9
p1
p1
Definition: p1EqnIrkPimple.H:24
solvePseudoUEqn.H
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DASimpleFoam::calcLduResiduals
void calcLduResiduals(volVectorField &URes, volScalarField &pRes, surfaceScalarField &phiRes)
calculate all the flow residuals except for the turbulence one
Definition: DASimpleFoam.C:1140
Foam::DASimpleFoam::invTranProdPEqn
void invTranProdPEqn(const volScalarField &mySource, volScalarField &pseudoP)
Inverse transpose product, Mp^(-T)
Definition: DASimpleFoam.C:1073
Foam::DASolver::runTimePtr_
autoPtr< Time > runTimePtr_
runTime pointer
Definition: DASolver.H:75
Foam::DASolver::meshPtr_
autoPtr< fvMesh > meshPtr_
fvMesh pointer
Definition: DASolver.H:78
solvePseudoPEqn.H
rAU
volScalarField rAU(1.0/UEqn.A())
Foam::DAUtility::isFieldReadable
static label isFieldReadable(const fvMesh &mesh, const word fieldName, const word fieldType)
Definition: DAUtility.C:817
createRefsSimple.H
UEqnSimple.H
Foam::DASimpleFoam::vec2Fields
void vec2Fields(const word mode, Vec cVec, volVectorField &UField, volScalarField &pField, surfaceScalarField &phiField, volScalarField &nuTildaField)
assign coupled petsc vec to fields or vice versa
Definition: DASimpleFoam.C:910
U1
U1
Definition: p1EqnIrkPimple.H:26
Foam::DASimpleFoam::calcAdjointResidual
void calcAdjointResidual(volVectorField &URes, volScalarField &pRes, surfaceScalarField &phiRes, volScalarField &nuTildaRes, Vec dFdW, volVectorField &UPsi, volScalarField &pPsi, surfaceScalarField &phiPsi, volScalarField &nuTildaPsi, volVectorField &adjURes, volScalarField &adjPRes, surfaceScalarField &adjPhiRes, volScalarField &adjNuTildaRes, label &cnt)
Definition: DASimpleFoam.C:1239
runTime
Time & runTime
Definition: createRefsHeatTransfer.H:1
pseudoUEqn
fvVectorMatrix pseudoUEqn(fvm::div(phi, pseudoU, divScheme) - fvm::laplacian(nuEff, pseudoU))
Foam::DASolver::printElapsedTime
void printElapsedTime(const Time &runTime, const label printToScreen)
Definition: DASolver.H:256
pEqnSimple.H
psi
const volScalarField & psi
Definition: createRefsRhoPimple.H:14
dafoam_plot3dtransform.mode
mode
Definition: dafoam_plot3dtransform.py:21
Foam::DASolver::calcAllFunctions
void calcAllFunctions(label print=0)
calculate the values of all objective functions and print them to screen
Definition: DASolver.C:204
p2
p2
Definition: p2EqnIrkPimple.H:24