# Description of AdjointShapeOptimizationFoam and how to implement new cost functions Ulf Nilsson

```Overview
Introduction
Theory
Implementation
Modifications
Results
and how to implement new cost functions
Ulf Nilsson
A project work in the course
CFD with OpenSource software,
taught by
H˚
akan Nilsson.
Chalmers University of Technology,
Gothenburg, Sweden
December 10, 2013
Ulf Nilsson
December 10, 2013
1 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Overview
Ulf Nilsson
I
Introduction
I
Theory
I
Implementation
I
Modifications
I
Results
December 10, 2013
2 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Introduction
Optimization in CFD
I
I
Increased importance
Optimization algorithm
-
I
Simplex
Evolution Strategies
...
Demanding
- Number of variables
- Constraints
Ulf Nilsson
December 10, 2013
3 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Introduction
Ulf Nilsson
I
I
Tool to compute sensitivity w.r.t. design variables
I
Well behaved objective function
I
Large number of design variables
I
Constraints, number of functions
December 10, 2013
4 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Objective function
minimize J = J(v, p, α)
(1)
such that
(v · ∇)v + ∇p − ∇ · (2νD(v)) + αv
= 0,
(2)
∇·v
= 0.
(3)
I
Gerneral cost function
I
I
Porosity α
I
Topology optimization
Ulf Nilsson
December 10, 2013
5 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Lagrangian relaxation
Z
minimize L := J +
(u, q)RdΩ.
(4)
Ω
I
Lagrange multipliers, (u, q)
I
A penalty to violate the constraints
I
Adjoint velocity, u, and pressure, q
I
No physical meaning
Ulf Nilsson
December 10, 2013
6 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Variations w.r.t. flow and design variables
δL = δα L + δv L + δp L,
(5)
δv L + δp L = 0.
I
I
I
I
(6)
Total variation.
Choose adjoint variables so that the variation with respect to flow
variables vanishes.
Equation of the adjoint variables and an expression for the sensitivity
of L with respect to the porosity of each cell, derivations out of the
scope of the project.
Approximations needed to arrive at the adjoint equations.
- Forzen turbulence.
- Ducted flow specialization.
I
Integration by parts to arrive at volume contribution and boundary
contribution.
Ulf Nilsson
December 10, 2013
7 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Sensitivity
∂L
= ui · vi Vi ,
∂αi
(7)
I
Expression for cell i
I
Vi the cell volume.
I
Equation of the adjoint variables and an expression for the sensitivity
of L with respect to the porosity of each cell
Ulf Nilsson
December 10, 2013
8 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
−2D(u)v = −∇q + ∇ · (2νD(u)) − αu −
∇·u=
I
∂JΩ
,
∂v
∂JΩ
.
∂p
(8)
(9)
Contribution from volumetric part of the cost function, JΩ .
Ulf Nilsson
December 10, 2013
9 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Wall and inlet boundary conditions
ut = 0,
(10)
∂JΩ
,
un = −
∂p
n · ∇q = 0.
(11)
(12)
Outlet boundary conditions
q = u · v + un vn + ν(n · ∇)un +
0 = vn ut + ν(n · ∇)ut +
∂JΓ
.
∂vt
∂JΓ
,
∂vn
(13)
(14)
December 10, 2013
Ulf Nilsson
10 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
v
p
Wall
No-slip
Inlet
Prescribed value
Outlet
p=0
Table: Boundary conditions for the primal quantities v and p.
I
Only valid for specific primal boundary condition.
I
Contribution from both the volume and boundary part of the cost
function, JΩ and JΓ respectively.
I
Eq. 13 and 14 used to calculate the adjoint pressure and tangential
part of the adjoint velocity at the outlet.
December 10, 2013
Ulf Nilsson
11 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Summary of the theory
I
Introducing a Lagrange relaxation, with multipliers u and q.
I
for the sensititivity.
I
Lengthy derivations.
I
Approximations and requirements to ensure a valid method.
I
If the cost function satisfies Jω = 0 the dependence of J enters the
solver only through the boundary conditions.
I
The optimization problem reduces to solving two flow equations, the
I
When the sensitivity is calculated an algorithm is needed to update
the porosity.
December 10, 2013
Ulf Nilsson
12 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Steepest descent algorithm
Steepest descent
pk = −∇f (xk ),
(15)
αn+1 = αn − ui · vi Vi δ.
(16)
I
General search direction in the opposite direction of the gradient of
the objective function.
I
Considering a linear problem a minimum will be found for small
enough steps in the search direction.
I
Using the sensitivity, ∂L/∂α to update.
I
The step length δ. An optimization problem of its own.
I
Under relaxation and limits in the final implementation.
December 10, 2013
Ulf Nilsson
13 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Solution procedure
Main-loop
Solve the
primal system
(v, p)
“Frozen tubulence”
Update the
porosity
Solve the
Steepest descent
(u, q)
Convergence
assessment
Converged
End time or
End
convergence condition
Figure 1: Block scheme of the solution procedure used in adjointShapeOptimizationFoam.
December 10, 2013
Ulf Nilsson
14 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
\$FOAM SOLVERS/incompressible/
Let us study the implementation together, trying to identify the underlying
theory we have covered so far. In the order it is implemented in the code.
I
Porosity update
- Implementation of the steepest descent using under relaxation.
- Correct?
I
The primal pressure-velocity SIMPLE corrector
- How is the sink term implemented?
I
- Is the current implementation done for objective functions with
volumteric contribution (JΩ 6= 0)?
December 10, 2013
Ulf Nilsson
15 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Porosity update
94
95
96
a l p h a +=
mesh . f i e l d R e l a x a t i o n F a c t o r ( " alpha " )
∗ ( min ( max ( a l p h a + lambda ∗ ( Ua & U) , z e r o A l p h a ) , alphaMax ) − a l p h a ) ;
αn+1 = αn (1 − γ) + γmin (max ((αn + ui · vi Vi δ) .0) , αmax ) ,
(17)
December 10, 2013
Ulf Nilsson
16 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Porosity update
However the there seems to be an incorrect sign in the steepest descent
implementation. The correct eq. and implementation should be
αn+1 = αn (1 − γ) + γmin (max ((αn − ui · vi Vi δ) , 0) , αmax ) .
94
95
96
(18)
a l p h a +=
mesh . f i e l d R e l a x a t i o n F a c t o r ( " alpha " )
∗ ( min ( max ( a l p h a − lambda ∗ ( Ua & U) , z e r o A l p h a ) , alphaMax ) − a l p h a ) ;
December 10, 2013
Ulf Nilsson
17 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Primal SIMPLE loop
103
104
105
106
107
108
109
110
111
112
113
114
// Momentum p r e d i c t o r
tmp<f v V e c t o r M a t r i x > UEqn
(
fvm : : d i v ( p h i , U)
+ t u r b u l e n c e −>d i v D e v R e f f (U)
+ fvm : : Sp ( a l p h a , U)
);
UEqn ( ) . r e l a x ( ) ;
s o l v e ( UEqn ( ) == −f v c : : g r a d ( p ) ) ;
December 10, 2013
Ulf Nilsson
18 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Primal SIMPLE loop
I
fvm::Sp(alpha, U) Implicit implementation of the extra
source term.
December 10, 2013
Ulf Nilsson
19 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
156 // A d j o i n t Momentum p r e d i c t o r
157
158
v o l V e c t o r F i e l d a d j o i n t T r a n s p o s e C o n v e c t i o n ( ( f v c : : g r a d ( Ua ) & U ) ) ;
159
// v o l V e c t o r F i e l d a d j o i n t T r a n s p o s e C o n v e c t i o n
160
// (
161
// f v c : : r e c o n s t r u c t
162
// (
163
//
mesh . magSf ( ) ∗ ( f v c : : sn Grad ( Ua ) & f v c : : i n t e r p o l a t e (U) )
164
// )
165
// ) ;
166
167
zeroCells ( adjointTransposeConvection , i n l e t C e l l s );
168
169
tmp<f v V e c t o r M a t r i x > UaEqn
170
(
171
fvm : : d i v (− p h i , Ua )
172
173
+ t u r b u l e n c e −>d i v D e v R e f f ( Ua )
174
+ fvm : : Sp ( a l p h a , Ua )
175
);
176
177
UaEqn ( ) . r e l a x ( ) ;
178
179
s o l v e ( UaEqn ( ) == −f v c : : g r a d ( pa ) ) ;
December 10, 2013
Ulf Nilsson
20 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
I
Similar to the primal system.
I
function, i.e. current implementation treats cost functions of
the type JΓ = 0
December 10, 2013
Ulf Nilsson
21 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
changed indicates additional problems with the current solver. Since the
actual solver seems to agree with the theory, what remains to be examined
is the boundary conditions.
First the boundary condition for a specific cost function need to be
calculated.
December 10, 2013
Ulf Nilsson
22 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Power dissipation
Z
J := −
1
(p + v 2 )v · n dΓ.
2
(19)
Γ
I
Only boundary contribution.
I
Energy loss, due to pressure drop and kinetic energy.
December 10, 2013
Ulf Nilsson
23 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Inlet boundary conditions
ut = 0,
un = −
(20)
∂JΩ
,
∂p
n · ∇q = 0.


ut




un




n · ∇q
=(
0,
=
0
vn
(21)
(22)
at wall,
at inlet,
(23)
= 0.
December 10, 2013
Ulf Nilsson
24 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Inlet boundary conditions
I
Use existing conditions
- Prescribed value
- No slip
I
Could be useful to implement a inlet boundary condition for
the adjoint velocity, if the velocity profile of the primal velocity
at the inlet is not easily described.
December 10, 2013
Ulf Nilsson
25 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Outlet boundary conditions
q = u · v + un vn + ν(n · ∇)un +
0 = vn ut + ν(n · ∇)ut +
(
q
0
∂JΓ
,
∂vn
(24)
∂JΓ
.
∂vt
(25)
= u · v + un vn + ν(n · ∇)un − 21 v 2 − vn2 ,
= vn (ut − vt ) + ν(n · ∇)ut .
(26)
December 10, 2013
Ulf Nilsson
26 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Outlet boundary conditions
I
Use eq. 26 to prescribe a value to the adjoint pressure
I
Use eq. 26 to prescribe a value to the tangential component
ut =
ν
−∆
uneigh,t − vp,n vp,t
,
ν
vp,n + ∆
(27)
December 10, 2013
Ulf Nilsson
27 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Implementation of pressure condition
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
// ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Member F u n c t i o n s
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ //
v o i d Foam : : a d j o i n t O u t l e t P r e s s u r e F v P a t c h S c a l a r F i e l d : : u p d a t e C o e f f s ( )
{
i f ( updated ( ) )
{
return ;
}
const f v s P a t c h F i e l d <s c a l a r >& p h i p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phi" ) ;
const f v s P a t c h F i e l d <s c a l a r >& p h i a p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phia" ) ;
const f v P a t c h F i e l d <v e c t o r >& Up =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("U" ) ;
const f v P a t c h F i e l d <v e c t o r >& Uap =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("Ua" ) ;
o p e r a t o r ==(( p h i a p / p a t c h ( ) . magSf ( ) − 1 . 0 ) ∗ p h i p / p a t c h ( ) . magSf ( ) + ( Up & Uap ) ) ;
fixedValueFvPatchScalarField : : updateCoeffs ( ) ;
}
December 10, 2013
Ulf Nilsson
28 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Implementation of velocity condition
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
// ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Member F u n c t i o n s
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ //
// Update t h e c o e f f i c i e n t s a s s o c i a t e d w i t h t h e p a t c h f i e l d
v o i d Foam : : a d j o i n t O u t l e t V e l o c i t y F v P a t c h V e c t o r F i e l d : : u p d a t e C o e f f s ( )
{
i f ( updated ( ) )
{
return ;
}
const f v s P a t c h F i e l d <s c a l a r >& p h i a p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phia" ) ;
const f v P a t c h F i e l d <v e c t o r >& Up =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("U" ) ;
s c a l a r F i e l d Un ( mag ( p a t c h ( ) . n f ( ) & Up ) ) ;
v e c t o r F i e l d UtHat ( ( Up − p a t c h ( ) . n f ( ) ∗ Un ) / ( Un + SMALL ) ) ;
v e c t o r F i e l d Uan ( p a t c h ( ) . n f ( ) ∗ ( p a t c h ( ) . n f ( ) & p a t c h I n t e r n a l F i e l d ( ) ) ) ;
v e c t o r F i e l d : : o p e r a t o r =( p h i a p ∗ p a t c h ( ) . S f ( ) / s q r ( p a t c h ( ) . magSf ( ) ) + UtHat ) ;
// v e c t o r F i e l d : : o p e r a t o r =(Uan + UtHat ) ;
fixedValueFvPatchVectorField : : updateCoeffs ( ) ;
}
December 10, 2013
Ulf Nilsson
29 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Analyzing the implementation
Listing 5 gives an implementation of the pressure condition according to
q = (un − 1)vn + u · v.
(28)
1
q = u · v + un vn + ν(n · ∇)un − v 2 − vn2 ,
2
(29)
compare to
Listing 6 gives an implementation of the velocity condition according to
up,t =
vp − vp,n
.
up,n + SMALL
(30)
compare to
ut =
ν
−∆
uneigh,t − vp,n vp,t
,
ν
vp,n + ∆
(31)
December 10, 2013
Ulf Nilsson
30 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Analyzing the theory
I
Not equal to the theory.
I
Another cost function, total pressure loss.
I
Approximations or derivations done makes the cost function
hard to identify.
December 10, 2013
Ulf Nilsson
31 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Analyzing the theory
I
Changing the sign of the steepest descent implementation.
I
Implementing the cost functions according to derived
equations.
December 10, 2013
Ulf Nilsson
32 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
I
I
Similar to the existing implementation
- Velocity in the neighbouring node
- Effective viscosity
const s c a l a r F i e l d & d e l t a i n v = p a t c h ( ) . d e l t a C o e f f s ( ) ; // d i s t a n c e ˆ( −1) between p a t c h and n e i g h b o u r i n g node .
// c r e a t e t h e o b j e c t n e e d e d t o g e t t h e v i s c o u s i t y :
const i n c o m p r e s s i b l e : : RASModel& r a s M o d e l =
db ( ) . l o o k u p O b j e c t <i n c o m p r e s s i b l e : : RASModel>(" RASProperties " ) ;
// n u e f f from t h e t u r u l e n c e model , r a s M o d e l :
s c a l a r F i e l d n u e f f = rasModel . nuEff ( ) ( ) . boundaryField ( ) [ patch ( ) . index ( ) ] ;
// N e i g h b o u r i n g node ’ s v e l o c i t y ( n o r m a l component ) :
s c a l a r F i e l d U a n e i g h n = ( Uap . p a t c h I n t e r n a l F i e l d ( ) & p a t c h ( ) . n f ( ) ) ;
December 10, 2013
Ulf Nilsson
33 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
I
Vector fields instead of scalar fields
// p a t c h I n t e r n a l F i e l d t o g e t t h e a d j o i n t v e l o c i t y o f n e i g h b o u r i n g node .
v e c t o r F i e l d U a n e i g h = Uap . p a t c h I n t e r n a l F i e l d ( ) ;
// I t s t a n g e n t i a l and n o r m a l components
v e c t o r F i e l d Uaneigh n = ( Uaneigh & normal )∗ normal ;
v e c t o r F i e l d Uaneigh t = Uaneigh − Uaneigh n ;
December 10, 2013
Ulf Nilsson
34 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Case settings
The boundary conditions of the primal and adjoint variables are set
according to the table below. More detailed information and m4 script of
the “box” example case can be found in the provided files.
u
q
Wall
Prescribed value, 0
Inlet
Prescribed value, un = vn
Outlet
December 10, 2013
Ulf Nilsson
35 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily porosity field
December 10, 2013
Ulf Nilsson
36 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily value of the objective function
0.006
”valueObjFunc.xy”
0.005
J
0.004
0.003
0.002
0.001
0
0
500
1000
1500
2000
2500
3000
Iteration step
December 10, 2013
Ulf Nilsson
37 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily value of the objective function, without
porosity update
0.006
”valueObjFunc.xy”
0.005
J
0.004
0.003
0.002
0.001
0
0
500
1000
1500
2000
2500
3000
Iteration step
December 10, 2013
Ulf Nilsson
38 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
Pipe bend example
December 10, 2013
Ulf Nilsson
39 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
Thank you for listening!
Questions
December 10, 2013
Ulf Nilsson
40 /
40
```