adjointShapeOptimisationFoam.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-2023 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  adjointShapeOptimisationFoam
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 optimisation 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 "argList.H"
49 #include "timeSelector.H"
50 #include "viscosityModel.H"
52 #include "simpleControl.H"
53 #include "pressureReference.H"
54 #include "findRefCell.H"
55 #include "constrainPressure.H"
56 #include "constrainHbyA.H"
57 #include "adjustPhi.H"
58 #include "fvModels.H"
59 #include "fvConstraints.H"
60 
61 #include "fvcDdt.H"
62 #include "fvcGrad.H"
63 #include "fvcFlux.H"
64 
65 #include "fvmDdt.H"
66 #include "fvmDiv.H"
67 #include "fvmLaplacian.H"
68 
69 using namespace Foam;
70 
71 
72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73 
74 template<class Type>
75 inline void zeroCells
76 (
77  VolField<Type>& vf,
78  const labelList& cells
79 )
80 {
81  forAll(cells, i)
82  {
83  vf[cells[i]] = Zero;
84  }
85 }
86 
87 
88 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
89 
90 int main(int argc, char *argv[])
91 {
92  #include "postProcess.H"
93 
94  #include "setRootCase.H"
95  #include "createTime.H"
96  #include "createMesh.H"
97  #include "createControl.H"
98  #include "createFields.H"
99  #include "initContinuityErrs.H"
100  #include "initAdjointContinuityErrs.H"
101 
102  turbulence->validate();
103 
104  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
105 
106  Info<< "\nStarting time loop\n" << endl;
107 
108  while (simple.loop(runTime))
109  {
110  Info<< "Time = " << runTime.userTimeName() << nl << endl;
111 
112  viscosity->lookup("lambda") >> lambda;
113 
114  // alpha +=
115  // mesh.relaxationFactor("alpha")
116  // *(lambda*max(Ua & U, zeroSensitivity) - alpha);
117  alpha +=
118  mesh.solution().fieldRelaxationFactor("alpha")
119  *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
120 
122  // zeroCells(alpha, outletCells);
123 
124  turbulence->predict();
125 
126  // Pressure-velocity SIMPLE corrector
127  {
128  fvModels.correct();
129 
130  // Momentum predictor
131 
133  (
134  fvm::div(phi, U)
135  + turbulence->divDevSigma(U)
136  + fvm::Sp(alpha, U)
137  ==
138  fvModels.source(U)
139  );
140  fvVectorMatrix& UEqn = tUEqn.ref();
141 
142  UEqn.relax();
143 
145 
146  solve(UEqn == -fvc::grad(p));
147 
149 
150  volScalarField rAU(1.0/UEqn.A());
152  tUEqn.clear();
154  adjustPhi(phiHbyA, U, p);
155 
156  // Update the pressure BCs to ensure flux consistency
158 
159  // Non-orthogonal pressure corrector loop
160  while (simple.correctNonOrthogonal())
161  {
162  fvScalarMatrix pEqn
163  (
165  );
166 
167  pEqn.setReference
168  (
171  );
172 
173  pEqn.solve();
174 
176  {
177  phi = phiHbyA - pEqn.flux();
178  }
179  }
180 
181  #include "continuityErrs.H"
182 
183  // Explicitly relax pressure for momentum corrector
184  p.relax();
185 
186  // Momentum corrector
187  U = HbyA - rAU*fvc::grad(p);
188  U.correctBoundaryConditions();
190  }
191 
192  // Adjoint Pressure-velocity SIMPLE corrector
193  {
194  // Adjoint Momentum predictor
195 
196  volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
197  // volVectorField adjointTransposeConvection
198  //(
199  // fvc::reconstruct
200  // (
201  // mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
202  // )
203  //);
204 
205  zeroCells(adjointTransposeConvection, inletCells);
206 
207  tmp<fvVectorMatrix> tUaEqn
208  (
209  fvm::div(-phi, Ua)
210  - adjointTransposeConvection
211  + turbulence->divDevSigma(Ua)
212  + fvm::Sp(alpha, Ua)
213  ==
214  fvModels.source(Ua)
215  );
216  fvVectorMatrix& UaEqn = tUaEqn.ref();
217 
218  UaEqn.relax();
219 
220  fvConstraints.constrain(UaEqn);
221 
222  solve(UaEqn == -fvc::grad(pa));
223 
225 
226  volScalarField rAUa(1.0/UaEqn.A());
227  volVectorField HbyAa("HbyAa", Ua);
228  HbyAa = rAUa*UaEqn.H();
229  tUaEqn.clear();
230  surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
231  adjustPhi(phiHbyAa, Ua, pa);
232 
233  // Non-orthogonal pressure corrector loop
234  while (simple.correctNonOrthogonal())
235  {
236  fvScalarMatrix paEqn
237  (
238  fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
239  );
240 
241  paEqn.setReference(paRefCell, paRefValue);
242  paEqn.solve();
243 
245  {
246  phia = phiHbyAa - paEqn.flux();
247  }
248  }
249 
250  #include "adjointContinuityErrs.H"
251 
252  // Explicitly relax pressure for adjoint momentum corrector
253  pa.relax();
254 
255  // Adjoint momentum corrector
256  Ua = HbyAa - rAUa*fvc::grad(pa);
257  Ua.correctBoundaryConditions();
259  }
260 
261  viscosity->correct();
262  turbulence->correct();
263 
264  runTime.write();
265 
266  Info<< "ExecutionTime = "
267  << runTime.elapsedCpuTime()
268  << " s\n\n" << endl;
269  }
270 
271  Info<< "End\n" << endl;
272 
273  return 0;
274 }
275 
276 
277 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Calculates and prints the continuity errors.
int main(int argc, char *argv[])
void zeroCells(VolField< Type > &vf, const labelList &cells)
For cases which do no have a pressure boundary adjust the balance of fluxes to obey continuity....
Generic GeometricField class.
void relax(const scalar alpha)
Relax field (for steady-state solution).
Finite volume constraints.
Definition: fvConstraints.H:67
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:603
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
tmp< volScalarField > A() const
Return the central coefficient.
Definition: fvMatrix.C:808
tmp< VolField< Type > > H() const
Return the H operation source.
Definition: fvMatrix.C:868
tmp< SurfaceField< Type > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:964
void setReference(const label celli, const Type &value, const bool forceReference=false)
Set reference level for solution.
Definition: fvMatrix.C:563
Finite volume models.
Definition: fvModels.H:65
tmp< fvMatrix< Type > > source(const VolField< Type > &field) const
Return source for an equation.
virtual void correct()
Correct the fvModels.
Definition: fvModels.C:380
bool correctNonOrthogonal(const bool finalIter=false)
Non-orthogonal corrector loop.
bool finalNonOrthogonalIter() const
Flag to indicate the last non-orthogonal iteration.
Provides controls for the pressure reference in closed-volume simulations.
scalar refValue() const
Return the pressure reference level.
label refCell() const
Return the cell in which the reference pressure is set.
bool loop(Time &time)
Time loop loop.
Definition: simpleControl.C:80
A class for managing temporary objects.
Definition: tmp.H:55
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Abstract base class for all fluid physical properties.
Definition: viscosity.H:50
tmp< fvVectorMatrix > tUEqn(fvm::div(phi, U)+MRF.DDt(rho, U)+turbulence->divDevTau(U)==fvModels.source(rho, U))
fvVectorMatrix & UEqn
Definition: UEqn.H:11
Calculates and prints the continuity errors.
simpleControl simple(mesh)
Find the reference cell nearest (in index) to the given cell but which is not on a cyclic,...
Calculate the first temporal derivative.
Calculate the face-flux of the given field.
Calculate the gradient of the given field.
Calculate the matrix for the first temporal derivative.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
const cellShapeList & cells
Declare and initialise the cumulative adjoint continuity error.
Declare and initialise the cumulative continuity error.
volScalarField rAU(1.0/UEqn.A())
U
Definition: pEqn.H:72
volVectorField & HbyA
Definition: pEqn.H:13
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho) *fvc::flux(HbyA))
dimensionedScalar alphaMax(viscosity->lookup("alphaMax"))
const labelList & inletCells
Definition: createFields.H:94
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0)
dimensionedScalar lambda(viscosity->lookup("lambda"))
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
tmp< SurfaceField< typename innerProduct< vector, Type >::type > > flux(const VolField< Type > &vf)
Return the face-flux field obtained from the given volVectorField.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
Namespace for OpenFOAM.
void constrainPressure(volScalarField &p, const RhoType &rho, const volVectorField &U, const surfaceScalarField &phiHbyA, const RAUType &rhorAU, const MRFType &MRF)
static const zero Zero
Definition: zero.H:97
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:35
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
messageStream Info
void constrainHbyA(volVectorField &HbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:34
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
static const char nl
Definition: Ostream.H:266
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
Execute application functionObjects to post-process existing results.
autoPtr< incompressible::momentumTransportModel > turbulence(incompressible::momentumTransportModel::New(U, phi, viscosity))
volScalarField & p