Dear Community,
Greetings for the day!!
i am trying to run a two fluid 1 solid topology optimisation problem with the objective of minimising pressure dissipation energy while adhering to max volume constraints.
I’m facing an issue in my FEniCS-based optimization problem where I’m adjusting a function ( \rho_f ) across a domain to meet volume constraints. Initially, ( \rho_f ) is uniformly set to 0.5 throughout the domain. However, as the optimization progresses, ( \rho_f ) evolves to have values around 1 in half of the domain and 0 in the other half.
During each iteration, I update ( \rho_f ) using an adjoint method and MMA algorithm. However, despite significant changes in ( \rho_f ), the values of volume constraints, denoted as fval, remain constant at [0.3, 0.3].
Upon inspecting the exported values of ( \rho_f ), I confirmed that they align with the expected distribution in the final iteration. However, despite this, fval doesn’t reflect these changes. The second term in the constraint definition always evaluates to 0.2 as it’s a constant, but the first term should reflect the updates in ( \rho_f ), as it’s calculated using the assemble function, which integrates ( \rho_f ) over the domain.
This inconsistency leads to a solution where the volume constraint is not respected at all. I’m seeking insights or suggestions on why fval isn’t updating as expected and how to rectify this issue.
I am using GCMMA-MMA-Python for MMA calculation. The projection, filtering and interpolation functions are taken from the paper.
The geometry is a rectangle with two inlets on the left boundary and 2 outlets on the right boundary in the same half of the domain as the inlet. I have prescribed inlet and outlet pressures and no-slip velocity at the walls as boundary conditions.
Minimal Code:
mesh = Mesh(RectangleMesh(MPI.comm_world, Point(0.0, 0.0), Point(lx,ly), N, N, 'crossed'))
rho = Function(DensitySpace) # rho= 1 corrosponds to fluid 1, rho= 0 corrosponds to fluid 2
rho_f = Function(DensitySpace)
# Optimization Problem #
ObjFunctional1 = 0.5 * (alpha1(projection(rho_f, eta_i)) * inner(u1, u1) * dx) + ((mufluid/2) * inner(grad(u1), grad(u1)) * dx)
ObjFunctional2 = 0.5 * (alpha2(projection(rho_f, eta_i)) * inner(u2, u2) * dx) + ((mufluid/2) * inner(grad(u2), grad(u2)) * dx)
ObjFunctional = ObjFunctional1 + ObjFunctional2
NSR1 = Kcm*rhofluid * inner(dot(u1, nabla_grad(u1)), v1) * dx\
+ mufluid * inner(grad(u1), grad(v1)) * dx\
- inner(grad(p1), v1) * dx + inner(div(u1), q1) * dx + alpha1(projection(rho_f, eta_i)) * inner(u1, v1) * dx
NSR2 = Kcm*rhofluid * inner(dot(u2, nabla_grad(u2)), v2) * dx\
+ mufluid * inner(grad(u2), grad(v2)) * dx\
- inner(grad(p2), v2) * dx + inner(div(u2), q2) * dx + alpha2(projection(rho_f, eta_i)) * inner(u2, v2) * dx
vol_frac_f1 = Constant(0.2)
vol_frac_f2 = Constant(0.2)
vol_constraint_f1 = (projection(rho_f, eta_i)/Domain_area) * dx - (vol_frac_f1/Domain_area) * dx(domain = mesh) # constraint definition
vol_constraint_f2 = ((1 - projection(rho_f, eta_i))/Domain_area) * dx - (vol_frac_f2/Domain_area) * dx(domain = mesh)
ddx_vol_f1 = derivative(vol_constraint_f1, rho_f)
ddx_vol_f2 = derivative(vol_constraint_f2, rho_f)
assign(rho, interpolate(Constant(0.5), DensitySpace))
#OPTIMIZATION LOOP
while inner_count <= 500 and objectiveConvergence == False:
rho_f = pdefilter(rho, rho_f)
###### SOLVE forward problem for fluid 1
###### SOLVE forward problem for fluid 2
###### SOLVE adjoint problem for fluid 1
###### SOLVE adjoint problem for fluid 2
f0val = assemble(ObjFunctional) # Objective function value
fval[0] = assemble(vol_constraint_f1) # Volume constraint fluid 1
fval[1] = assemble(vol_constraint_f2) # Volume constraint fluid 2
unfilteredGradient.vector()[:] = assemble(ddx)[:]
filteredGradient = pdefilter(unfilteredGradient, filteredGradient)
unfiltered_volumegradient_1.vector()[:] = assemble(ddx_vol_f1)[:]
filtered_volumegradient_1 = pdefilter(unfiltered_volumegradient_1, filtered_volumegradient_1)
unfiltered_volumegradient_2.vector()[:] = assemble(ddx_vol_f2)[:]
filtered_volumegradient_2 = pdefilter(unfiltered_volumegradient_2, filtered_volumegradient_2)
df0dx[:,0] = filteredGradient.vector()
dfdx[0,:] = filtered_volumegradient_1.vector()[:]
dfdx[1,:] = filtered_volumegradient_2.vector()[:]
#Output stuff
rho.vector() = xmma[:,0].copy() # updating design variable to output from MMA
Any help would be greatly appreciated.
Thank you