# Double Pipe TopOp Problem Using Linear Elements

**URL:** https://fenicsproject.discourse.group/t/double-pipe-topop-problem-using-linear-elements/12091
**Category:** dolfin-adjoint
**Created:** [August 31, 2023, 1:12pm UTC](https://fenicsproject.discourse.group/t/double-pipe-topop-problem-using-linear-elements/12091 "2023-08-31T13:12:24Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![Morten](https://avatars.discourse-cdn.com/v4/letter/m/cdc98d/32.png) [@Morten](https://fenicsproject.discourse.group/u/Morten)
#### Post date: [August 31, 2023, 1:12pm UTC](https://fenicsproject.discourse.group/t/double-pipe-topop-problem-using-linear-elements/12091/1 "2023-08-31T13:12:24Z")

</div>

I’m trying to solve the double pipe problem from [Borrval and Petersson (2003)](https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.426), using the same methods as in [A detailed introduction to density-based topology optimisation of fluid flow problems with implementation in MATLAB](https://link.springer.com/article/10.1007/s00158-022-03420-9).

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.

```auto
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()

```