Dear everyone,
I’m trying to solve the double pipe problem from Borrval and Petersson (2003) , using dolfin_adjoint and MMA (the code for MMA can be found here).
When solving the forward problem before starting optimisation loop the resulting presure plot is as expected without fault as seen (left), however when running the optimisation loop it results in the plot to the right.
The full code results in the correct final design, with the only error being the pressure plot. MY MWE are:
from fenics import *
from fenics_adjoint import *
import numpy as np
import matplotlib.pyplot as plt
from MMA import mmasub # Import MMA functions
fenics.parameters["reorder_dofs_serial"] = False
## Newton solver parameters
nlmax = 25 # Maximum number of iterations
nltol = 1e-6 # Newton solver tolerance
## Physical Parameters
Uin = 1; rho = 1e-3; mu = 1 # Inlet velocity, density and viscosity
## Domain Size and Meshing
ly = 1; lx = 1
nely = 102; nelx = int(nely*lx/ly) # Number of elements in the x and y direction
def Problem(lx: float, ly: float, nelx: float, nely: float, Uin: float):
y_l = ly/6
mesh = UnitSquareMesh(nelx, nely, diagonal="crossed") # Creating a unitsquare mesh
P1 = VectorElement('P', mesh.ufl_cell(), 1)
P2 = FiniteElement('P', mesh.ufl_cell(), 1)
element = MixedElement([P1, P2])
V = FunctionSpace(mesh, element)
## Creating the inlet class
class FluidInlet_Lower(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0) and (x[1] > ly/6 and x[1] < (2*ly)/6)
class FluidInlet_Upper(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0) and x[1] > (4*ly)/6 and x[1] < (5*ly)/6
## Creating the outlet class
class FluidOutlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], lx) and ((x[1] > ly/6 and x[1] < (2*ly)/6) or (x[1] > (4*ly)/6 and x[1] < (5*ly)/6))
## Creating a class contaning all outer walls
class AllWalls(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
## Creates a domain containing all facets
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
## Initialising of classes
allwalls = AllWalls()
fluidinlet_lower = FluidInlet_Lower()
fluidinlet_upper = FluidInlet_Upper()
fluidoutlet = FluidOutlet()
## Marking of classes
allwalls.mark(facets, 4) # Marks all outer walls with 4
fluidinlet_lower.mark(facets, 1) # Overwrites and marks the lower inlet with 1
fluidinlet_upper.mark(facets, 2) # Overwrites and marks the upper inlet with 2
fluidoutlet.mark(facets, 3) # Overwrites and marks outlet with 3
u_in_exp_lower = Expression('Uin*(-4*pow((x[1]/y_l-1), 2)+4*(x[1]/y_l-1))', Uin=Uin, y_l=y_l, degree=2)
u_in_exp_upper = Expression('Uin*(-4*pow(((x[1]-3*y_l)/y_l-1), 2)+4*((x[1]-3*y_l)/y_l-1))', Uin=Uin, y_l=y_l, degree=2)
## Boundary conditions
BcV_Inlet_Lower = DirichletBC(V.sub(0).sub(0), u_in_exp_lower, facets, 1)
BcV_Inlet_Upper = DirichletBC(V.sub(0).sub(0), u_in_exp_upper, facets, 2)
BcV_Inlet1 = DirichletBC(V.sub(0).sub(1), Constant(0.0), facets, 1)
BcV_Inlet2 = DirichletBC(V.sub(0).sub(1), Constant(0.0), facets, 2)
BcV_Outlet = DirichletBC(V.sub(0).sub(1), Constant(0.0), facets, 3)
BcV_Wall = DirichletBC(V.sub(0), Constant((0.0, 0.0)), facets, 4)
BcP_Outlet = DirichletBC(V.sub(1), Constant(0.0), facets, 3)
BcV = [BcV_Outlet, BcV_Wall, BcV_Inlet_Lower, BcV_Inlet_Upper, BcV_Inlet1, BcV_Inlet2]
BcP = [BcP_Outlet]
bcs = BcV + BcP # Collecting all boundary conditions
return V, mesh, bcs
V, mesh, bcs = Problem(lx, ly, nelx, nely, Uin) # Function for defining the problem, result in the functionspace V, mesh and dirichelt boundary conditions
V0 = FunctionSpace(mesh, "DG", 0)
## Allowable fluid volume fraction
volfrac = 1/3
## Brinkman penalisation
alphamax = 2.5*mu/(0.01**2); alphamin = 2.5*mu/(100*2)
ainit = 2.5*mu/(0.1**2) # Initial brinkman penalty
qa = (-volfrac*(alphamax-alphamin) - ainit+alphamax)/(volfrac*(ainit-alphamin)) # Initial intpolation factor
## Optimisation loop setting
niter = 1 # Maximum number of iterations
mvlim = 0.2 # Maximum change of gamma_k per iteration
dx2 = dx(metadata={"quadrature_degree":2*3}) # Increases mathmatical degree to increase stability
def forward(gamma_tilde, qa):
# ----- DEFINING PROBLEM -----
# Defining the trial and test functions
w, q = TestFunctions(V) # Creating test function for the velocity and pressure functionspace
up = Function(V) # Creating a conbined function for velocity and pressure
# Material interpolation
def RAMP(gamma_tilde, qa): # RAMP method interpolating
return alphamin + (alphamax - alphamin)*(1-gamma_tilde)/(1+qa*gamma_tilde)
## Assign initial conditions
e_u0 = Expression(('Uin', '0.0'), degree=1, Uin=Uin) # Initial velocity expression
e_p0 = Expression('0.0', degree=1) # Initial pressure expression
u0 = interpolate(e_u0, V.sub(0).collapse()) # Interpolates the initial values over functionspace
p0 = interpolate(e_p0, V.sub(1).collapse())# Interpolates the initial values over functionspace
assign(up, [u0, p0]) # Assigns the initial values to the functionspace
u, p = split(up) # Splits into two seperate functionspaces
## Stabilisation parameter
vnorm = sqrt(dot(u, u)) # Calculates the velocity norm
tau_1 = h/(2*vnorm) # Calculates tau_1 SUPG/PSPG stability factor
tau_3 = (rho*h**2)/(12*mu) # Calculates tau_3 SUPG/PSPG sstability factor
tau_4 = rho/RAMP(gamma_tilde, qa) # Calculates tau_4 SUPG/PSPG sstability factor
tau = (tau_1**(-2) + tau_3**(-2) + tau_4**(-2))**(-1/2) # Calculating the collected SUPG/PSPG parameter
## Navier-Stokes equation
F = (rho*inner(grad(u)* u, w)*dx2 + # Convection term
mu*inner(grad(u), grad(w))*dx2 - # Viscous term
inner(p, div(w))*dx2 + # Pressure term
div(u)*q*dx2 + # Divergence term
inner(w*RAMP(gamma_tilde, qa), u)*dx2
)
## SUPG/PSPG stabilisation
F_SUPG = inner(tau*grad(u)*w, (rho*grad(u)*u + grad(p)+ RAMP(gamma_tilde, qa)*u))*dx2
F_PSPG = inner(tau/rho*grad(q), (grad(u)*u + grad(p) + RAMP(gamma_tilde, qa)*u))*dx2
F += F_SUPG + F_PSPG
J = derivative(F, up) # Calculates the Jacobian
problem = NonlinearVariationalProblem(F, up, bcs, J) # Sets up the problem
solver = NonlinearVariationalSolver(problem) # Saves the solver
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['relaxation_parameter'] = 1.
prm['newton_solver']['relative_tolerance'] = nltol
prm['newton_solver']['absolute_tolerance'] = 1e-12
prm['newton_solver']['maximum_iterations'] = nlmax
prm['newton_solver']['error_on_nonconvergence'] = True
prm['newton_solver']['linear_solver'] = 'mumps
solver.solve() # Solves the problem
u, p = up.split(deepcopy=True)
PHI = (0.5 * inner(RAMP(gamma_tilde, qa) * u, u)*dx2 + 0.5* mu * inner(grad(u), grad(u)) * dx2) # Calculates the dissipated energy
return PHI, p
def Optimization(gamma, rf_f0, rf_h):
with pyadjoint.stop_annotating() as _:
# ------ INITIALASING LOOP -----
## MMA Topology Optimization settings
# Initialize MMA parameters
m = 1 # Number of constraints
n = gamma.vector()[:].size # Number of design variables
xval = gamma.vector()[:].reshape(-1, 1) # Initial design variables
xold1 = xval # Design variables, one iterations ago, initial
xold2 = xval # Design variables, two iterations ago, initial
low = 0
upp = 0
NumberOfConstraints = 1 # Only constraint: Volume constraint
a0 = 1
ai = 0*np.ones((NumberOfConstraints, 1))
c = 1000*np.ones((NumberOfConstraints, 1))
d = 1*np.ones((NumberOfConstraints, 1))
## Load inital designs and preallocate sensitivities
gamma_k = gamma
df0dgamma = Function(V0, name="Object function Sensitivity")
dhdgamma = Function(V0, name="Constraint function Sensitivity")
## Performance of inital design
f0_k = rf_f0(gamma_k) # Objective function
h_k = rf_h(gamma_k) # Constraint function
Mnd = assemble(4*gamma_k*(1-gamma_k)*dx2)/(lx*ly)*100 #Measure of non-discreteness
## Initiate iteration counters and previous objective
iter = 0
## OPTIMIZATION LOOP
while (iter < niter):
iter += 1 # Update iteration counter
# Evaluate sensitivities
df0dgamma.assign(rf_f0.derivative())
dhdgamma.assign(rf_h.derivative())
# Scaling functions and storing them in matrices for the MMA solver
xval = gamma_k.vector()[:].reshape(-1, 1) # Saving design variable in xval
if iter == 1: # Scaling with the initial compliance
scale = np.abs(f0_k/10)
f0val = 1 + f0_k/scale # Scaling f0 at xval
df0dgammaarr = df0dgamma.vector()[:].reshape(-1, 1)/scale # nx1 matrix of sensitivities of f0 at xval
hval = h_k # Array of constraint functions at xval
dhdgammaarr = dhdgamma.vector()[:].reshape(m, n) # nx1 matrix of sensitivities of h at xval
# Setting min and max limits for the xval of next iteration
xmin = np.maximum(0, gamma_k.vector()[:].reshape(-1, 1)-mvlim)
xmax = np.minimum(1, gamma_k.vector()[:].reshape(-1, 1)+mvlim)
# Call MMA
xmma, ymma, zmma, lam, xsi, eta, mu1, zet, s, low, upp = \
mmasub(m, n, iter, xval, xmin, xmax, xold1, xold2,
f0val, df0dgammaarr, hval, dhdgammaarr, low, upp, a0, ai, c, d, move=1)
# Update density field
gamma_k.vector()[:] = xmma.flatten() # Overwrite gamma with the new gamma from the MMA solver
_, p = forward(gamma_k, qa)
plt.figure(1)
plt.clf()
plt.colorbar(plot(p))
plt.title(f"Pressure at iter {iter}")
plt.show()
rho = Constant(rho) # Creating the density as a constant
mu = Constant(mu) # Creating the viscosity as a constant
h = CellDiameter(mesh) # Diameter of each element
gamma = interpolate(Constant(volfrac),V0) # Unfiltered design variable
PHI, p = forward(gamma, qa) # Solving the forward problem for initial values
plt.colorbar(plot(p)) # Plot of the initial pressure feild
plt.title('Pressure plot before optimisation loop')
plt.show()
# Objective and Constraint function
h = assemble(gamma*dx2)/(volfrac) - 1 # Volume constraint
f0 = assemble(PHI) # Objective function of dissipated energy
control = Control(gamma) # Control function
# Reduced functionals
rf_f0 = ReducedFunctional(f0, control)
rf_h = ReducedFunctional(h, control)
Optimization(gamma, rf_f0, rf_h)