adjointShapeOptimizationFoam.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 | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  ajointShapeOptimizationFoam
26 
27 Description
28  Steady-state solver for incompressible, turbulent flow of non-Newtonian
29  fluids with optimisation of duct shape by applying "blockage" in regions
30  causing pressure loss as estimated using an adjoint formulation.
31 
32  References:
33  \verbatim
34  "Implementation of a continuous adjoint for topology optimization of
35  ducted flows"
36  C. Othmer,
37  E. de Villiers,
38  H.G. Weller
39  AIAA-2007-3947
40  http://pdf.aiaa.org/preview/CDReadyMCFD07_1379/PV2007_3947.pdf
41  \endverbatim
42 
43  Note that this solver optimises for total pressure loss whereas the
44  above paper describes the method for optimising power-loss.
45 
46 \*---------------------------------------------------------------------------*/
47 
48 #include "fvCFD.H"
51 #include "simpleControl.H"
52 
53 template<class Type>
54 void zeroCells
55 (
56  GeometricField<Type, fvPatchField, volMesh>& vf,
57  const labelList& cells
58 )
59 {
60  forAll(cells, i)
61  {
62  vf[cells[i]] = pTraits<Type>::zero;
63  }
64 }
65 
66 
67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
68 
69 int main(int argc, char *argv[])
70 {
71  #include "setRootCase.H"
72 
73  #include "createTime.H"
74  #include "createMesh.H"
75 
76  simpleControl simple(mesh);
77 
78  #include "createFields.H"
79  #include "initContinuityErrs.H"
81 
82  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
83 
84  Info<< "\nStarting time loop\n" << endl;
85 
86  while (simple.loop())
87  {
88  Info<< "Time = " << runTime.timeName() << nl << endl;
89 
90  laminarTransport.lookup("lambda") >> lambda;
91 
92  //alpha +=
93  // mesh.relaxationFactor("alpha")
94  // *(lambda*max(Ua & U, zeroSensitivity) - alpha);
95  alpha +=
96  mesh.fieldRelaxationFactor("alpha")
97  *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
98 
100  //zeroCells(alpha, outletCells);
101 
102  // Pressure-velocity SIMPLE corrector
103  {
104  // Momentum predictor
105 
106  tmp<fvVectorMatrix> UEqn
107  (
108  fvm::div(phi, U)
109  + turbulence->divDevReff(U)
110  + fvm::Sp(alpha, U)
111  );
112 
113  UEqn().relax();
114 
115  solve(UEqn() == -fvc::grad(p));
116 
117  volScalarField rAU(1.0/UEqn().A());
118  volVectorField HbyA("HbyA", U);
119  HbyA = rAU*UEqn().H();
120  UEqn.clear();
122  (
123  "phiHbyA",
124  fvc::interpolate(HbyA) & mesh.Sf()
125  );
126  adjustPhi(phiHbyA, U, p);
127 
128  // Non-orthogonal pressure corrector loop
129  while (simple.correctNonOrthogonal())
130  {
131  fvScalarMatrix pEqn
132  (
134  );
135 
136  pEqn.setReference(pRefCell, pRefValue);
137  pEqn.solve();
138 
139  if (simple.finalNonOrthogonalIter())
140  {
141  phi = phiHbyA - pEqn.flux();
142  }
143  }
144 
145  #include "continuityErrs.H"
146 
147  // Explicitly relax pressure for momentum corrector
148  p.relax();
149 
150  // Momentum corrector
151  U = HbyA - rAU*fvc::grad(p);
152  U.correctBoundaryConditions();
153  }
154 
155  // Adjoint Pressure-velocity SIMPLE corrector
156  {
157  // Adjoint Momentum predictor
158 
159  volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
160  //volVectorField adjointTransposeConvection
161  //(
162  // fvc::reconstruct
163  // (
164  // mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
165  // )
166  //);
167 
168  zeroCells(adjointTransposeConvection, inletCells);
169 
170  tmp<fvVectorMatrix> UaEqn
171  (
172  fvm::div(-phi, Ua)
173  - adjointTransposeConvection
174  + turbulence->divDevReff(Ua)
175  + fvm::Sp(alpha, Ua)
176  );
177 
178  UaEqn().relax();
179 
180  solve(UaEqn() == -fvc::grad(pa));
181 
182  volScalarField rAUa(1.0/UaEqn().A());
183  volVectorField HbyAa("HbyAa", Ua);
184  HbyAa = rAUa*UaEqn().H();
185  UaEqn.clear();
186  surfaceScalarField phiHbyAa
187  (
188  "phiHbyAa",
189  fvc::interpolate(HbyAa) & mesh.Sf()
190  );
191  adjustPhi(phiHbyAa, Ua, pa);
192 
193  // Non-orthogonal pressure corrector loop
194  while (simple.correctNonOrthogonal())
195  {
196  fvScalarMatrix paEqn
197  (
198  fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
199  );
200 
201  paEqn.setReference(paRefCell, paRefValue);
202  paEqn.solve();
203 
204  if (simple.finalNonOrthogonalIter())
205  {
206  phia = phiHbyAa - paEqn.flux();
207  }
208  }
209 
210  #include "adjointContinuityErrs.H"
211 
212  // Explicitly relax pressure for adjoint momentum corrector
213  pa.relax();
214 
215  // Adjoint momentum corrector
216  Ua = HbyAa - rAUa*fvc::grad(pa);
217  Ua.correctBoundaryConditions();
218  }
219 
220  laminarTransport.correct();
221  turbulence->correct();
222 
223  runTime.write();
224 
225  Info<< "ExecutionTime = "
226  << runTime.elapsedCpuTime()
227  << " s\n\n" << endl;
228  }
229 
230  Info<< "End\n" << endl;
231 
232  return 0;
233 }
234 
235 
236 // ************************************************************************* //
Calculates and prints the continuity errors.
singlePhaseTransportModel laminarTransport(U, phi)
rhoEqn solve()
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:34
phiHbyA
Definition: pcEqn.H:74
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
autoPtr< compressible::turbulenceModel > turbulence
Definition: createFields.H:23
dynamicFvMesh & mesh
const scalar pRefValue
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
fvVectorMatrix UEqn(fvm::ddt(rho, U)+fvm::div(phi, U)+turbulence->divDevRhoReff(U)==fvOptions(rho, U))
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
volScalarField & p
Definition: createFields.H:51
Declare and initialise the cumulative ddjoint continuity error.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"))
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0)
surfaceScalarField & phi
const cellShapeList & cells
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
volScalarField rAU(1.0/UEqn.A())
List< label > labelList
A List of labels.
Definition: labelList.H:56
const labelList & inletCells
Definition: createFields.H:95
const label pRefCell
const dictionary & simple
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
zeroCells(alpha, inletCells)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
int main(int argc, char *argv[])
Definition: postCalc.C:54
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
U
Definition: pEqn.H:82
HbyA
Definition: pEqn.H:7