Modifying a value during optimisation

[Nlanglet]

Hello, In order to carry out topology optimisation simulations in fluid mechanics using the density based method, I need to modify a q value over the iterations, used in the k function: def k(alpha): return k_max + (k_min - k_max) alpha (1. + q) / (alpha + q) However, I can’t get the optimisation to take account of the change in the value of q. I get exactly the same result without changing the value of q. I also tried to define k as def(k,q) but without success.

Here is part of the code:

# Initial setup for topology optimization
alpha = interpolate(Constant(f_V-0.01), A) # Initial guess for the design variable
set_working_tape(Tape())              # Clear all annotations and restart the adjoint model

# Solve the simulation
u = solve_forward_problem(alpha)

# Visualization files
alpha_pvd_file = File("%s/alpha_iterations.pvd" %(output_folder)); alpha_viz = Function(A, name = "AlphaVisualisation")
dj_pvd_file = File("%s/dj_iterations.pvd" %(output_folder)); dj_viz = Function(A, name = "dJVisualisation")

# Callback during topology optimization
global current_iteration
current_iteration = 0
def derivative_cb_post(j, dj, current_alpha):
    global current_iteration
    if current_iteration >=10:
        q = 0.1
    print("\n [Iteration: %d] J = %1.7e, q = %1.2f\n" %(current_iteration, j, q))
    obj.append(j)
    ité.append(current_iteration)
    current_iteration += 1
    # Save for visualization
    alpha_viz.assign(current_alpha); alpha_pvd_file << alpha_viz
    dj_viz.assign(dj); dj_pvd_file << dj_viz

if current_iteration >=10:
    q = 0.1

# Objective function
u_split = split(u); v = u_split[0]
J = assemble((1/2.*(mu_)*inner(grad(v) + grad(v).T, grad(v) + grad(v).T) + inner(k(alpha) * v, v))*dx)
#J = assemble(w_ud*dot(v - u_target, v - u_target) * dx(6) + inner(k(alpha) * v, v)*dx)

print(" Current objective function value: %1.7e" %(J))

# Set the topology optimization problem and solver
alpha_C = Control(alpha)
Jhat = ReducedFunctional(J, alpha_C, derivative_cb_post = derivative_cb_post)
problem_min = MinimizationProblem(Jhat, bounds = (0.0, 1.0), constraints = [UFLInequalityConstraint((f_V - alpha)*dx, alpha_C)])

solver_opt = IPOPTSolver(problem_min, parameters = {'maximum_iterations': 50 })

# Perform topology optimization
alpha_opt = solver_opt.solve()
alpha.assign(alpha_opt); alpha_viz.assign(alpha)
alpha_pvd_file << alpha_viz

Thank you !

[dham]

There are a few things to note.

First, the code is incomplete. Without the code from solve_forward_problem there is no way we can tell what you are doing.

However, the assignments q=0.1 are just local variable bindings. They have nothing to do with the problem that the topology optimisation solves.

An effective strategy would have to change something on the tape, but without the problem code, I can't guess what that would be.

[Nlanglet]

Thank you for your response, I admit I have difficulty understanding what exactly the Tape function of pyadjoint does. Here is the function that create the mesh and then the function that solves the direct problem, the rest of the code consists of defining the boundary conditions.

# Create the 2D mesh and plot it
N_mesh = 40
delta_x = x_max - x_min; delta_y = y_max - y_min
mesh = RectangleMesh(Point(x_min, y_min), Point(x_max, y_max), int(N_mesh*delta_x/delta_y), N_mesh, diagonal = "crossed")
File('%s/mesh.pvd' %(output_folder)) << mesh

# Function spaces -> MINI element (2D)
V1_element =  FiniteElement('Lagrange', mesh.ufl_cell(), 1)
B_element = FiniteElement('Bubble', mesh.ufl_cell(), 3)
V_element = VectorElement(NodalEnrichedElement(V1_element, B_element)) # Velocity
P_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1) # Pressure
if flow_regime == 'laminar':
    U_element = MixedElement([V_element, P_element])
U = FunctionSpace(mesh, U_element) # Mixed function space
A_element = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
A = FunctionSpace(mesh, A_element) # Design variable function space (nodal)

and

# Function to solve the forward problem
def solve_forward_problem(alpha):
    # Set the state vector and test functions
    u = Function(U); u.rename("StateVariable", "StateVariable")
    u_split = split(u); v = u_split[0]; p = u_split[1]
    w = TestFunction(U)
    w_split = split(w); w_v = w_split[0]; w_p = w_split[1]

    # Additional definitions
    n = FacetNormal(mesh)
    nu_ = mu_/rho_

    # Weak form of the pressure-velocity formulation
    I_ = as_tensor(np.eye(2)); T = -p*I_ + (mu_)*(grad(v) + grad(v).T); v_mat = v
    F = div(v) * w_p*dx + inner( grad(w_v), T )*dx + rho_ * inner( dot(grad(v), v), w_v )*dx + inner(k(alpha) * v_mat, w_v)*dx

    # Boundary conditions
    if Probtype == "PipeBend":
        bcs = [DirichletBC(U.sub(0), wall_velocity_value, boundary_markers, marker_numbers['wall']),
            DirichletBC(U.sub(0), inlet_velocity_expression, boundary_markers, marker_numbers['inlet'])]

    if Probtype == "DoublePipe":
        bcs = [DirichletBC(U.sub(0), wall_velocity_value, boundary_markers, marker_numbers['wall']),
            DirichletBC(U.sub(0), inlet_velocity_expression1, boundary_markers, marker_numbers['inlet1']),
            DirichletBC(U.sub(0), inlet_velocity_expression2, boundary_markers, marker_numbers['inlet2']),
            DirichletBC(U.sub(0), inlet_velocity_expression1, boundary_markers, marker_numbers['outlet1']),
            DirichletBC(U.sub(0), inlet_velocity_expression2, boundary_markers, marker_numbers['outlet2'])]

    if Probtype == "Diffuser2D":
        bcs = [DirichletBC(U.sub(0), inlet_velocity_expression, boundary_markers, marker_numbers['inlet']),
            DirichletBC(U.sub(0).sub(1), Constant((0.0)), boundary_markers, marker_numbers['outlet']),
            DirichletBC(U.sub(1), Constant(0.0), boundary_markers, marker_numbers['outlet']),
            DirichletBC(U.sub(0), wall_velocity_value, boundary_markers, marker_numbers['wall'])]

    # Solve the simulation
    dF = derivative(F, u); problem = NonlinearVariationalProblem(F, u, bcs, dF)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    return u

Thank you !