I’m trying to solve the double pipe problem from Borrval and Petersson (2003), using the same methods as in A detailed introduction to density-based topology optimisation of fluid flow problems with implementation in MATLAB.
My code is able to solve the optimisation problem with similar results for the geometry and object value. However when solving the final forward problem, it results in a checkerboard pressure plot.
When solving problem with Taylor-Hood elements, the problem disappears. And i have been able to solve the stationary problem of the initial problem (no topology optimsation), using first order elements.
Therefore i suspect that it has something to do with the implementation of the SPUG/PSPG stabilisation in combination the Brinkman penalty.
from dolfin import *
from dolfin_adjoint import *
import numpy as np
import matplotlib.pyplot as plt
from MMA import mmasub # Import MMA functions
from pyadjoint.placeholder import Placeholder
set_log_level(LogLevel.ERROR) # Only print message if error message
# turning reorder_dofs_serial to avoid scrambling of solution when extracted at the cost of some computing power
# https://fenicsproject.org/qa/1522/recent-change-in-numbering-convention-for-dofs/
parameters["reorder_dofs_serial"] = False
## Physical Parameters
Uin = 1; rho = 1e-3; mu = 1 # Inlet velocity, density and viscosity
lx = 1; ly = 1 # Length and height of domain
## Domain Size and Meshing
nely = 102; nelx = int(nely*lx/ly) # Number of elements in the x and y direction
y_l = ly/6 # Inlet length
## Meshing and functionspaces
mesh = RectangleMesh(Point(0.0, 0.0), Point(lx, ly), nelx, nely, diagonal="crossed") # Creating a rectangle 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
V0 = FunctionSpace(mesh, "DG", 0) # Discontinous functionspace
## Allowable fluid volume fraction
volfrac = 1/3; xinit = volfrac
voltot = lx*ly
## Brinkman penalisation
alphamax = 2.5*mu/(0.01**2); alphamin = 2.5*mu/(100*2)
## Continuation strategy
ainit = 2.5*mu/(0.1**2) # Initial brinkman penalty
qinit = (-xinit*(alphamax-alphamin) - ainit+alphamax)/(xinit*(ainit-alphamin)) # Initial intpolation factor
qavec = qinit/np.array([[1], [2], [10], [20]]); qanum = len(qavec)-1
qastep = 0
qa = Constant(np.squeeze(qavec[qastep]))
Placeholder(qa)
## Optimisation loop setting
niter = 5 # Max number of iterations
tol_opt = 1e-3 # Convergence tolreance, based on the absolute error
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, qa):
"""
Solves the forward problem (solving PDE) using the RAMP method
Inputs: gamma is the density variable field
Outputs: Variables used for computing the objectiive function
"""
# ----- 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, qa): # RAMP method interpolating
return alphamin + (alphamax - alphamin)*(1-gamma)/(1+qa*gamma)
## 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, 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, qastep), u)*dx2
)
## SUPG/PSPG stabilisation
F_SUPG = tau*(rho*inner(grad(u)*w, grad(u)*u) + inner(grad(u)*w, grad(p)) + RAMP(gamma, qa)*inner(grad(u)*w, u))*dx2 #
F_PSPG = tau/rho*(inner(grad(q), grad(u)*u) + inner(grad(q), grad(p)) + RAMP(gamma, qa)*inner(grad(q), 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
solver.solve() # Solves the problem
u, p = up.split()
PHI = (0.5 * inner(RAMP(gamma, qa) * u, u)*dx2 + 0.5* mu * inner(grad(u), grad(u)) * dx2) # Calculates the dissipated energy
return PHI, u, p # return variables for computing f0
def Optimization(gamma, rf_f0, rf_h):
"""
Performs gradient based topology optimazation using the Methode of moving assymptotes
based on a MMA and GCMMA script by Krister Svanberg (https://people.kth.se/~krille/mmagcmma.pdf)
which is translated to Python Arjee Deetman (https://github.com/arjendeetman/GCMMA-MMA-Python)
Inputs: gamma is the density variable field, rf_f0 is the reduced functional for f0, rf_h is the reduced functional for h, f0_history is
Outputs: gamma
"""
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
## Initiate iteration counters and previous objective
iter = 0
f0_kprev = 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 old design values and update gamma
xold2 = xold1 # Design variables, two iterations ago
xold1 = xval # Design variables, one iterations ago
# Update density field
gamma_k.vector()[:] = xmma.flatten() # Overwrite gamma with the new gamma from the MMA solver
# Performance of the current design
if iter > 1:
f0_kprev = f0_k
f0_k = rf_f0(gamma_k) # Objective function
h_k = rf_h(gamma_k) # Constraint function
return gamma_k
if __name__ == '__main__':
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
# ----- INITIALIZE DENSITIES -----
gamma = interpolate(Constant(volfrac),V0) # Unfiltered design variable
PHI, u, p = forward(gamma, qastep) # Solve forward problem (physics problem)
# Objective and Constraint function
h = assemble(gamma*dx2)/(volfrac*voltot) - 1
f0 = assemble(PHI)
control = Control(gamma)
# Reduced functionals
rf_f0 = ReducedFunctional(f0, control)
rf_h = ReducedFunctional(h, control)
gamma = Optimization(gamma, rf_f0, rf_h)
PHI, u, p = forward(gamma, qastep) # Solve forward problem (physics problem)
plt.figure(4)
plt.colorbar(plot(p))
plt.title("Pressure at final design")
plt.show()