adjointShapeOptimizationFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 "setRootCaseLists.H"
75  #include "createTime.H"
76  #include "createMesh.H"
77  #include "createControl.H"
78  #include "createFields.H"
79  #include "initContinuityErrs.H"
81 
82  turbulence->validate();
83 
84  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
85 
86  Info<< "\nStarting time loop\n" << endl;
87 
88  while (simple.loop(runTime))
89  {
90  Info<< "Time = " << runTime.timeName() << nl << endl;
91 
92  laminarTransport.lookup("lambda") >> lambda;
93 
94  // alpha +=
95  // mesh.relaxationFactor("alpha")
96  // *(lambda*max(Ua & U, zeroSensitivity) - alpha);
97  alpha +=
98  mesh.fieldRelaxationFactor("alpha")
99  *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
100 
102  // zeroCells(alpha, outletCells);
103 
104  // Pressure-velocity SIMPLE corrector
105  {
106  // Momentum predictor
107 
108  tmp<fvVectorMatrix> tUEqn
109  (
110  fvm::div(phi, U)
111  + turbulence->divDevReff(U)
112  + fvm::Sp(alpha, U)
113  ==
114  fvOptions(U)
115  );
116  fvVectorMatrix& UEqn = tUEqn.ref();
117 
118  UEqn.relax();
119 
120  fvOptions.constrain(UEqn);
121 
122  solve(UEqn == -fvc::grad(p));
123 
124  fvOptions.correct(U);
125 
126  volScalarField rAU(1.0/UEqn.A());
127  volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
128  tUEqn.clear();
130  adjustPhi(phiHbyA, U, p);
131 
132  // Update the pressure BCs to ensure flux consistency
134 
135  // Non-orthogonal pressure corrector loop
136  while (simple.correctNonOrthogonal())
137  {
138  fvScalarMatrix pEqn
139  (
141  );
142 
143  pEqn.setReference(pRefCell, pRefValue);
144  pEqn.solve();
145 
146  if (simple.finalNonOrthogonalIter())
147  {
148  phi = phiHbyA - pEqn.flux();
149  }
150  }
151 
152  #include "continuityErrs.H"
153 
154  // Explicitly relax pressure for momentum corrector
155  p.relax();
156 
157  // Momentum corrector
158  U = HbyA - rAU*fvc::grad(p);
159  U.correctBoundaryConditions();
160  fvOptions.correct(U);
161  }
162 
163  // Adjoint Pressure-velocity SIMPLE corrector
164  {
165  // Adjoint Momentum predictor
166 
167  volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
168  // volVectorField adjointTransposeConvection
169  //(
170  // fvc::reconstruct
171  // (
172  // mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
173  // )
174  //);
175 
176  zeroCells(adjointTransposeConvection, inletCells);
177 
178  tmp<fvVectorMatrix> tUaEqn
179  (
180  fvm::div(-phi, Ua)
181  - adjointTransposeConvection
182  + turbulence->divDevReff(Ua)
183  + fvm::Sp(alpha, Ua)
184  ==
185  fvOptions(Ua)
186  );
187  fvVectorMatrix& UaEqn = tUaEqn.ref();
188 
189  UaEqn.relax();
190 
191  fvOptions.constrain(UaEqn);
192 
193  solve(UaEqn == -fvc::grad(pa));
194 
195  fvOptions.correct(Ua);
196 
197  volScalarField rAUa(1.0/UaEqn.A());
198  volVectorField HbyAa("HbyAa", Ua);
199  HbyAa = rAUa*UaEqn.H();
200  tUaEqn.clear();
201  surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
202  adjustPhi(phiHbyAa, Ua, pa);
203 
204  // Non-orthogonal pressure corrector loop
205  while (simple.correctNonOrthogonal())
206  {
207  fvScalarMatrix paEqn
208  (
209  fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
210  );
211 
212  paEqn.setReference(paRefCell, paRefValue);
213  paEqn.solve();
214 
215  if (simple.finalNonOrthogonalIter())
216  {
217  phia = phiHbyAa - paEqn.flux();
218  }
219  }
220 
221  #include "adjointContinuityErrs.H"
222 
223  // Explicitly relax pressure for adjoint momentum corrector
224  pa.relax();
225 
226  // Adjoint momentum corrector
227  Ua = HbyAa - rAUa*fvc::grad(pa);
228  Ua.correctBoundaryConditions();
229  fvOptions.correct(Ua);
230  }
231 
232  laminarTransport.correct();
233  turbulence->correct();
234 
235  runTime.write();
236 
237  Info<< "ExecutionTime = "
238  << runTime.elapsedCpuTime()
239  << " s\n\n" << endl;
240  }
241 
242  Info<< "End\n" << endl;
243 
244  return 0;
245 }
246 
247 
248 // ************************************************************************* //
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
fv::options & fvOptions
surfaceScalarField & phi
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)
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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
tmp< volScalarField > rAU
Definition: initCorrectPhi.H:1
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
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:97
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0)
static const char nl
Definition: Ostream.H:265
Declare and initialise the cumulative ddjoint continuity error.
Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar("pos", dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar("neg", dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"))
U
Definition: pEqn.H:72
label pRefCell
Definition: createFields.H:106
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: [].
scalar pRefValue
Definition: createFields.H:107
simpleControl simple(mesh)
zeroField Sp
Definition: alphaSuSp.H:2