This page was generated from unit-7-optimization/03a_Shape_Derivative_Laplace_SemiAuto.ipynb.
7.4 PDE-Constrained Shape Optimization in NGSolve¶
- We want to solve the PDE-constrained shape optimization problem :nbsphinx-math:`begin{equation}
underset{Omegasubset mathsf{D}}{mbox{min}} ; J(u) := int_Omega |u-u_d|^q ; dx, quad qge 2
- end{equation}` subject to that \((\Omega,u)\) satisfy :nbsphinx-math:`begin{equation}
int_Omega nabla u cdot nabla v ; dx = int_Omega f v ; dx ; quad text{ for all } v in H_0^1(Omega),
end{equation}` where \(\Omega \subset \mathbb R^2\) for given \(u_d, f \in C^1(\mathbb R^2)\).
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry
from ngsolve.solvers import *
[2]:
geo = SplineGeometry()
geo.AddCircle(c=(0.5,0.5), r=0.5, bc = 'circle')
#geo.AddRectangle( (0.2, 0.2), (0.8, 0.8), bcs = ('b','r','t','l'))
mesh = Mesh(geo.GenerateMesh(maxh = 0.1))
[3]:
#given data of our problem (chosen such that \Omega^* = [0,1]^2 is the optimal shape)
f = CoefficientFunction(2*y*(1-y)+2*x*(1-x))
ud = x*(1-x)*y*(1-y)
grad_f = CoefficientFunction( (f.Diff(x), f.Diff(y) ) )
grad_ud = CoefficientFunction( (ud.Diff(x), ud.Diff(y) ) )
[4]:
fes = H1(mesh, order=2, dirichlet=".*")
gfu = GridFunction(fes)
scene_u = Draw (gfu, mesh, "state")
gfp = GridFunction(fes)
scene_p = Draw (gfp, mesh, "adjoint")
Note that (for linear problems) the operator on the left hand side of the adjoint equation is just the transpose of the state operator.
Automatic Shape Differentiation¶
- The formula for the shape derivative was derived as the partial derivative of the perturbed Lagrangian (brought back to the original domain): :nbsphinx-math:`begin{align*}
- dmathcal J(Omega; X) = frac{partial}{partial t} left( int_Omega |u - u_d^t|^q ,mbox{det}(F_t) mbox{d} x
int_{Omega} (F_t^{-top}nabla u) cdot (F_t^{-top} nabla p), mbox{det}(F_t) , dx - int_{Omega}f^t p ,mbox{det}(F_t) ,dx right)|_{t=0}
end{align*}` where
\(T_t(x)=x+tX(x)=y\)
\(F_t = DT_t = I+t DX\)
\(u_d^t = u_d \circ T_t\)
\(f^t = f \circ T_t\)
The integrand depends on the parameter \(t\) only via \(F_t\) and \(T_t\). We define the Lagrangian in this form, involving a parameter t that has the value 0.
[5]:
VEC = H1(mesh, order=2, dim=2)
PHI, X = VEC.TnT()
#define symbolic parameter t=0 in order to allow differentiation w.r.t. t
t = Parameter(0)
F = Id(2) + t * grad(X).trans
def Equation(u,v):
return ( (Inv(F.trans)*grad(u))*(Inv(F.trans)*grad(v))-f*v)*Det(F)
q=6
def CostAuto(u):
return (u-ud)**q*Det(F)
def CostAuto2(u):
return (u-ud)**q
Lagrangian = CostAuto(gfu) + Equation(gfu,gfp)
State equation¶
Equation can also be used to define the bilinear form. The following defines left and right hand side of the PDE in a “BilinearForm”:
[6]:
u, v = fes.TnT()
aAuto = BilinearForm(fes, symmetric=True)
aAuto += Equation(u,v)*dx
Now the PDE can be conveniently solved by calling Newton’s method (which terminates after one iteration since the PDE is linear)
[7]:
gfu.vec[:]=0
Newton(aAuto,gfu,freedofs=fes.FreeDofs())
scene_u.Redraw()
Newton iteration 0
err = 0.13016086186788947
Newton iteration 1
err = 9.481807124502225e-17
Adjoint equation¶
- We set up the adjoint equation :nbsphinx-math:`begin{align*}
mbox{Find } p in H_0^1(Omega): int_Omega nabla w cdot nabla p , mbox dx = - partial_u J(u)(w) quad text{ for all } w in H_0^1(Omega)
- end{align*}` where \(u\) is the solution to the state equation. For \(J(u) = \int_\Omega |u-u_d|^2 \mbox dx\), we get :nbsphinx-math:`begin{align*}
partial_u J(u)(w) = 2 int_Omega (u-u_d)w ,mbox dx.
end{align*}` However, we can also use the Diff(…) command:
[8]:
p, w = fes.TnT()
fadjoint = LinearForm(fes)
fadjoint += -1*(CostAuto(gfu)*dx).Diff(gfu,w)
[9]:
def SolveAdjointEquation():
rhs = gfp.vec.CreateVector()
rhs.data = fadjoint.vec - aAuto.mat.T * gfp.vec
update = gfp.vec.CreateVector()
update.data = aAuto.mat.Inverse(fes.FreeDofs()).T * rhs
gfp.vec.data += update
[10]:
fadjoint.Assemble()
SolveAdjointEquation()
scene_p.Redraw()
Automatic Shape Differentiation¶
- Denoting the integrand by \(G^{u,p}\), the shape derivative is given by :nbsphinx-math:`begin{align*}
dmathcal J(Omega; X) =& left( left. frac{partial G^{u,p}}{partial t} + frac{d G^{u,p}}{dy} cdot frac{d T_t}{dt}right)rightrvert_{t=0} \ =& left. frac{partial G^{u,p}}{partial t}rightrvert_{t=0} + frac{d G^{u,p}}{dy} cdot X
end{align*}`
[11]:
dJOmegaAuto = LinearForm(VEC)
dJOmegaAuto += Lagrangian.Diff(t) * dx
dJOmegaAuto += (Lagrangian.Diff(x, X[0]) + Lagrangian.Diff(y, X[1])) * dx
[12]:
b = BilinearForm(VEC)
b += InnerProduct(grad(X),grad(PHI))*dx + InnerProduct(X,PHI)*dx
gfX = GridFunction(VEC)
# gfset denotes the deformation of the original domain and will be updated during the shape optimization
gfset = GridFunction(VEC)
gfset.Set((0,0))
scene=Draw(gfset,mesh,"gfset")
[13]:
def SolveDeformationEquationAuto():
rhs = gfX.vec.CreateVector()
rhs.data = dJOmegaAuto.vec - b.mat * gfX.vec
update = gfX.vec.CreateVector()
update.data = b.mat.Inverse(VEC.FreeDofs()) * rhs
gfX.vec.data += update
[14]:
b.Assemble()
dJOmegaAuto.Assemble()
SolveDeformationEquationAuto()
Draw(-gfX, mesh, "-gfX")
[14]:
BaseWebGuiScene
[15]:
print('Cost at initial design', Integrate (CostAuto2(gfu), mesh))
gfset.Set((0,0))
scale = 0.5 / Norm(gfX.vec)
gfset.vec.data -= scale * gfX.vec
Cost at initial design 6.457553863133461e-13
[16]:
mesh.SetDeformation(gfset)
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs())
print('Cost at new design', Integrate (CostAuto2(gfu), mesh))
scene.Redraw()
Newton iteration 0
err = 0.1571174420960136
Newton iteration 1
err = 1.1609319199315723e-16
Cost at new design 1.5689143159982939e-13
Thus, the user has to enter the PDE (in its transformed form) only once.
Finally, let us again run the full algorithm:
[17]:
#reset to and solve for initial configuration
scene_u = Draw(gfu)
gfset.Set((0,0))
mesh.SetDeformation(gfset)
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs())
LineSearch = True
iter_max = 600
Jold = Integrate(CostAuto2(gfu), mesh)
converged = False
for k in range(iter_max):
print('cost at iteration', k, ': ', Jold)
mesh.SetDeformation(gfset)
scene_u.Redraw()
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs(), printing = False)
fadjoint.Assemble()
SolveAdjointEquation()
b.Assemble()
dJOmegaAuto.Assemble()
SolveDeformationEquationAuto()
scale = 0.01 / Norm(gfX.vec)
gfsetOld = gfset
gfset.vec.data -= scale * gfX.vec
Jnew = Integrate(CostAuto2(gfu), mesh)
if LineSearch:
while Jnew > Jold and scale > 1e-12:
#input('a')
scale = scale / 2
if scale <= 1e-12:
converged = True
break
gfset.vec.data = gfsetOld.vec - scale * gfX.vec
mesh.SetDeformation(gfset)
gfu.vec[:]=0
Newton(aAuto, gfu, fes.FreeDofs(), printing = False)
Jnew = Integrate(CostAuto2(gfu), mesh)
if converged==True:
break
Jold = Jnew
Redraw(blocking=True)
Newton iteration 0
err = 0.13016086186788947
Newton iteration 1
err = 9.481807124502225e-17
cost at iteration 0 : 6.457553863133461e-13
cost at iteration 1 : 5.983887299218896e-13
cost at iteration 2 : 5.335221619386397e-13
cost at iteration 3 : 4.750742838717911e-13
cost at iteration 4 : 4.224887358578915e-13
cost at iteration 5 : 3.7524777180142396e-13
cost at iteration 6 : 3.3287057919500564e-13
cost at iteration 7 : 2.9491159086510645e-13
cost at iteration 8 : 2.6095879603665375e-13
cost at iteration 9 : 2.3063205708646567e-13
cost at iteration 10 : 2.035814377293271e-13
cost at iteration 11 : 1.7948554798915158e-13
cost at iteration 12 : 1.5804991105732506e-13
cost at iteration 13 : 1.390053569757802e-13
cost at iteration 14 : 1.2210644796799573e-13
cost at iteration 15 : 1.0712994014784748e-13
cost at iteration 16 : 9.387328623151561e-14
cost at iteration 17 : 8.215318371514817e-14
cost at iteration 18 : 7.180417268866989e-14
cost at iteration 19 : 6.267728691993569e-14
cost at iteration 20 : 5.463876088363814e-14
cost at iteration 21 : 4.7568793746481066e-14
cost at iteration 22 : 4.136036852582171e-14
cost at iteration 23 : 3.591812007049436e-14
cost at iteration 24 : 3.115723823748801e-14
cost at iteration 25 : 2.7002381389221666e-14
cost at iteration 26 : 2.3386558633849078e-14
cost at iteration 27 : 2.0249916108851958e-14
cost at iteration 28 : 1.753833523841301e-14
cost at iteration 29 : 1.5201732827926245e-14
cost at iteration 30 : 1.3191992545794986e-14
cost at iteration 31 : 1.1460695749849676e-14
cost at iteration 32 : 9.957576406582612e-15
cost at iteration 33 : 8.632213301301494e-15
cost at iteration 34 : 7.442386197370637e-15
cost at iteration 35 : 6.3661987177662546e-15
cost at iteration 36 : 5.402363949733956e-15
cost at iteration 37 : 4.554359969711412e-15
cost at iteration 38 : 3.8192092703343004e-15
cost at iteration 39 : 3.1880716117992034e-15
cost at iteration 40 : 2.6497757284759458e-15
cost at iteration 41 : 2.192980672152958e-15
cost at iteration 42 : 1.807073039235905e-15
cost at iteration 43 : 1.4824524085938693e-15
cost at iteration 44 : 1.2105670389576063e-15
cost at iteration 45 : 9.8386230723172e-16
cost at iteration 46 : 7.956991898751652e-16
cost at iteration 47 : 6.402708251352465e-16
cost at iteration 48 : 5.125175552215949e-16
cost at iteration 49 : 4.080526005409263e-16
cost at iteration 50 : 3.2308490353689827e-16
cost at iteration 51 : 2.5436044958738657e-16
cost at iteration 52 : 1.9908817075715602e-16
cost at iteration 53 : 1.5490715358823715e-16
cost at iteration 54 : 1.1979726809862288e-16
cost at iteration 55 : 9.212807207889791e-17
cost at iteration 56 : 7.056068995185006e-17
cost at iteration 57 : 5.48598300413397e-17
cost at iteration 58 : 4.5620525712213855e-17
cost at iteration 59 : 4.0057472118834696e-17
cost at iteration 60 : 3.4702532171283616e-17
cost at iteration 61 : 3.037407306318569e-17
cost at iteration 62 : 2.6423128089062752e-17
cost at iteration 63 : 2.31418986707482e-17
cost at iteration 64 : 2.017942953154694e-17
cost at iteration 65 : 1.7699199067396037e-17
cost at iteration 66 : 1.5463241816146383e-17
cost at iteration 67 : 1.3585240797098516e-17
cost at iteration 68 : 1.1890357269593688e-17
cost at iteration 69 : 1.0464222098274457e-17
cost at iteration 70 : 9.174736905044622e-18
cost at iteration 71 : 8.088233597576724e-18
cost at iteration 72 : 7.10371850808119e-18
cost at iteration 73 : 6.273188400999334e-18
cost at iteration 74 : 5.518883715501759e-18
cost at iteration 75 : 4.881859065694522e-18
cost at iteration 76 : 4.301900055104616e-18
cost at iteration 77 : 3.8116153869111924e-18
cost at iteration 78 : 3.364123504576092e-18
cost at iteration 79 : 2.9854703797564188e-18
cost at iteration 80 : 2.6389590134674287e-18
cost at iteration 81 : 2.34550544821156e-18
cost at iteration 82 : 2.076233220445727e-18
cost at iteration 83 : 1.848022846562556e-18
cost at iteration 84 : 1.638036846255227e-18
cost at iteration 85 : 1.4599592095860127e-18
cost at iteration 86 : 1.295643557967857e-18
cost at iteration 87 : 1.1562237757843933e-18
cost at iteration 88 : 1.0272199524750795e-18
cost at iteration 89 : 9.177172945532354e-19
cost at iteration 90 : 8.161189943924632e-19
cost at iteration 91 : 7.298538412439378e-19
cost at iteration 92 : 6.496052580099691e-19
cost at iteration 93 : 5.814554792305476e-19
cost at iteration 94 : 5.179007907358631e-19
cost at iteration 95 : 4.639242007514461e-19
cost at iteration 96 : 4.1346972484571927e-19
cost at iteration 97 : 3.70620597959826e-19
cost at iteration 98 : 3.3048106493599525e-19
cost at iteration 99 : 2.963969430877318e-19
cost at iteration 100 : 2.6440463520760735e-19
cost at iteration 101 : 2.372457142368346e-19
cost at iteration 102 : 2.1170660861943507e-19
cost at iteration 103 : 1.9003446481002304e-19
cost at iteration 104 : 1.696195168733629e-19
cost at iteration 105 : 1.5230527370652611e-19
cost at iteration 106 : 1.359679725857099e-19
cost at iteration 107 : 1.2212251338461037e-19
cost at iteration 108 : 1.090360928276733e-19
cost at iteration 109 : 9.795681007809638e-20
cost at iteration 110 : 8.746619338228156e-20
cost at iteration 111 : 7.85961791016252e-20
cost at iteration 112 : 7.018100679592305e-20
cost at iteration 113 : 6.307764016567679e-20
cost at iteration 114 : 5.632367302210743e-20
cost at iteration 115 : 5.063433734347229e-20
cost at iteration 116 : 4.52112226787407e-20
cost at iteration 117 : 4.065445310144534e-20
cost at iteration 118 : 3.629835222580623e-20
cost at iteration 119 : 3.2649134532160286e-20
cost at iteration 120 : 2.914908379616561e-20
cost at iteration 121 : 2.622733558449903e-20
cost at iteration 122 : 2.3414490728912542e-20
cost at iteration 123 : 2.107599031769807e-20
cost at iteration 124 : 1.8815110740573947e-20
cost at iteration 125 : 1.694431678404717e-20
cost at iteration 126 : 1.512700123870333e-20
cost at iteration 127 : 1.3631341148759687e-20
cost at iteration 128 : 1.2170641758056992e-20
cost at iteration 129 : 1.0975949022874114e-20
cost at iteration 130 : 9.802075491719958e-21
cost at iteration 131 : 8.848932110333834e-21
cost at iteration 132 : 7.90582021728165e-21
cost at iteration 133 : 7.146616890740038e-21
cost at iteration 134 : 6.389181576826966e-21
cost at iteration 135 : 5.7857503347494835e-21
cost at iteration 136 : 5.177680113458935e-21
cost at iteration 137 : 4.699385935481947e-21
cost at iteration 138 : 4.21136681280123e-21
cost at iteration 139 : 3.833569483477622e-21
cost at iteration 140 : 3.441855861922309e-21
cost at iteration 141 : 3.1446727955247575e-21
cost at iteration 142 : 2.8299495204246234e-21
cost at iteration 143 : 2.597266196702026e-21
cost at iteration 144 : 2.3437657764262963e-21
cost at iteration 145 : 2.162453776838353e-21
cost at iteration 146 : 1.9572996450087074e-21
cost at iteration 147 : 1.816613241063529e-21
cost at iteration 148 : 1.6493516287711283e-21
cost at iteration 149 : 1.5404818931789125e-21
cost at iteration 150 : 1.4027492066767962e-21
cost at iteration 151 : 1.318514625569354e-21
cost at iteration 152 : 1.2037771004462694e-21
cost at iteration 153 : 1.1384090745979975e-21
cost at iteration 154 : 1.0417190669940193e-21
cost at iteration 155 : 9.906684523066128e-22
cost at iteration 156 : 9.084219266008396e-22
cost at iteration 157 : 8.681174510907018e-22
cost at iteration 158 : 7.978679870892573e-22
cost at iteration 159 : 7.654003911116294e-22
cost at iteration 160 : 7.058844769001284e-22
cost at iteration 161 : 6.78607123801438e-22
cost at iteration 162 : 6.303711206591396e-22
cost at iteration 163 : 6.053698260107228e-22
cost at iteration 164 : 5.713387159445664e-22
cost at iteration 165 : 5.438440714730507e-22
cost at iteration 166 : 5.217583948720182e-22
cost at iteration 167 : 4.909771009672683e-22
cost at iteration 168 : 4.757063605745574e-22
cost at iteration 169 : 4.452313250306558e-22
cost at iteration 170 : 4.345881717273167e-22
cost at iteration 171 : 4.0491712080306655e-22
cost at iteration 172 : 3.982523475627112e-22
cost at iteration 173 : 3.690175858589996e-22
cost at iteration 174 : 3.6668841420365605e-22
cost at iteration 175 : 3.3726716886980163e-22
[ ]: