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-2016 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 #include "fvOptions.H"
53 
54 template<class Type>
55 void zeroCells
56 (
57  GeometricField<Type, fvPatchField, volMesh>& vf,
58  const labelList& cells
59 )
60 {
61  forAll(cells, i)
62  {
63  vf[cells[i]] = Zero;
64  }
65 }
66 
67 
68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
69 
70 int main(int argc, char *argv[])
71 {
72  #include "postProcess.H"
73 
74  #include "setRootCase.H"
75  #include "createTime.H"
76  #include "createMesh.H"
77  #include "createControl.H"
78  #include "createFields.H"
79  #include "createFvOptions.H"
80  #include "initContinuityErrs.H"
82 
83  turbulence->validate();
84 
85  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
86 
87  Info<< "\nStarting time loop\n" << endl;
88 
89  while (simple.loop())
90  {
91  Info<< "Time = " << runTime.timeName() << nl << endl;
92 
93  laminarTransport.lookup("lambda") >> lambda;
94 
95  //alpha +=
96  // mesh.relaxationFactor("alpha")
97  // *(lambda*max(Ua & U, zeroSensitivity) - alpha);
98  alpha +=
99  mesh.fieldRelaxationFactor("alpha")
100  *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
101 
103  //zeroCells(alpha, outletCells);
104 
105  // Pressure-velocity SIMPLE corrector
106  {
107  // Momentum predictor
108 
109  tmp<fvVectorMatrix> tUEqn
110  (
111  fvm::div(phi, U)
112  + turbulence->divDevReff(U)
113  + fvm::Sp(alpha, U)
114  ==
115  fvOptions(U)
116  );
117  fvVectorMatrix& UEqn = tUEqn.ref();
118 
119  UEqn.relax();
120 
121  fvOptions.constrain(UEqn);
122 
123  solve(UEqn == -fvc::grad(p));
124 
125  fvOptions.correct(U);
126 
127  volScalarField rAU(1.0/UEqn.A());
128  volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
129  tUEqn.clear();
131  adjustPhi(phiHbyA, U, p);
132 
133  // Update the pressure BCs to ensure flux consistency
135 
136  // Non-orthogonal pressure corrector loop
137  while (simple.correctNonOrthogonal())
138  {
139  fvScalarMatrix pEqn
140  (
142  );
143 
144  pEqn.setReference(pRefCell, pRefValue);
145  pEqn.solve();
146 
147  if (simple.finalNonOrthogonalIter())
148  {
149  phi = phiHbyA - pEqn.flux();
150  }
151  }
152 
153  #include "continuityErrs.H"
154 
155  // Explicitly relax pressure for momentum corrector
156  p.relax();
157 
158  // Momentum corrector
159  U = HbyA - rAU*fvc::grad(p);
160  U.correctBoundaryConditions();
161  fvOptions.correct(U);
162  }
163 
164  // Adjoint Pressure-velocity SIMPLE corrector
165  {
166  // Adjoint Momentum predictor
167 
168  volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
169  //volVectorField adjointTransposeConvection
170  //(
171  // fvc::reconstruct
172  // (
173  // mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
174  // )
175  //);
176 
177  zeroCells(adjointTransposeConvection, inletCells);
178 
179  tmp<fvVectorMatrix> tUaEqn
180  (
181  fvm::div(-phi, Ua)
182  - adjointTransposeConvection
183  + turbulence->divDevReff(Ua)
184  + fvm::Sp(alpha, Ua)
185  ==
186  fvOptions(Ua)
187  );
188  fvVectorMatrix& UaEqn = tUaEqn.ref();
189 
190  UaEqn.relax();
191 
192  fvOptions.constrain(UaEqn);
193 
194  solve(UaEqn == -fvc::grad(pa));
195 
196  fvOptions.correct(Ua);
197 
198  volScalarField rAUa(1.0/UaEqn.A());
199  volVectorField HbyAa("HbyAa", Ua);
200  HbyAa = rAUa*UaEqn.H();
201  tUaEqn.clear();
202  surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
203  adjustPhi(phiHbyAa, Ua, pa);
204 
205  // Non-orthogonal pressure corrector loop
206  while (simple.correctNonOrthogonal())
207  {
208  fvScalarMatrix paEqn
209  (
210  fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
211  );
212 
213  paEqn.setReference(paRefCell, paRefValue);
214  paEqn.solve();
215 
216  if (simple.finalNonOrthogonalIter())
217  {
218  phia = phiHbyAa - paEqn.flux();
219  }
220  }
221 
222  #include "adjointContinuityErrs.H"
223 
224  // Explicitly relax pressure for adjoint momentum corrector
225  pa.relax();
226 
227  // Adjoint momentum corrector
228  Ua = HbyAa - rAUa*fvc::grad(pa);
229  Ua.correctBoundaryConditions();
230  fvOptions.correct(Ua);
231  }
232 
233  laminarTransport.correct();
234  turbulence->correct();
235 
236  runTime.write();
237 
238  Info<< "ExecutionTime = "
239  << runTime.elapsedCpuTime()
240  << " s\n\n" << endl;
241  }
242 
243  Info<< "End\n" << endl;
244 
245  return 0;
246 }
247 
248 
249 // ************************************************************************* //
autoPtr< compressible::turbulenceModel > turbulence
Definition: createFields.H:23
surfaceScalarField & phi
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
U
Definition: pEqn.H:83
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
void constrainPressure(volScalarField &p, const RhoType &rho, const volVectorField &U, const surfaceScalarField &phiHbyA, const RAUType &rhorAU, const MRFType &MRF)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< fvVectorMatrix > tUEqn(fvm::ddt(rho, U)+fvm::div(phi, U)+MRF.DDt(rho, U)+turbulence->divDevRhoReff(U)==fvOptions(rho, U))
zeroCells(alpha, inletCells)
const labelList & inletCells
Definition: createFields.H:95
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
HbyA
Definition: pcEqn.H:74
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:34
fv::options & fvOptions
phiHbyA
Definition: pcEqn.H:73
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
singlePhaseTransportModel laminarTransport(U, phi)
dynamicFvMesh & mesh
const cellShapeList & cells
rhoEqn solve()
List< label > labelList
A List of labels.
Definition: labelList.H:56
tmp< volVectorField > constrainHbyA(const tmp< volVectorField > &tHbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:33
static const zero Zero
Definition: zero.H:91
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0)
static const char nl
Definition: Ostream.H:262
Declare and initialise the cumulative ddjoint continuity error.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const scalar pRefValue
const label pRefCell
dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"))
const dictionary & simple
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
Calculates and prints the continuity errors.
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
messageStream Info
volScalarField & p
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Execute application functionObjects to post-process existing results.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
volScalarField rAU(1.0/UEqn.A())