Topology optimization loop cannot update the new design variable

Pasuta · March 6, 2024, 7:10am

I am trying to do topology optimization where design variable is eps but i cannot update the new eps for calculation in the next iteration as error showing , this is my full code Calculate sensitivity equation of reaction-diffusion systems with dolfinx
and this is my eps updating equation
Screenshot 2024-03-06 134931

where eps(new) and eps(old) are the porosity of the current and previous iterations, could somebody suggest me or give me a simple example

# Function to compute the objective
def J_eps(m_n):
    eps.x.array[:] = m_n
    solver.solve(c)
    J_old = fem.assemble_scalar(J_form)

    return J_old

# Function to compute the gradient of the reaction rate w.r.t. eps
def grad_J(m_n,dJde=dJde):
    eps.x.array[:] = m_n

    # Direct solution
    solver.solve(c)

    # Adjoint solution
    lamda = problem_adj.solve()

    # Reassemble dRde
    dRde.zeroEntries()
    fem.petsc.assemble_matrix(dRde, dRde_form)
    dRde.assemble()

    # Compute gradient
    # zero gradient first
    with dJde.localForm() as dJde_loc:
        dJde_loc.set(0)

    fem.petsc.assemble_vector(dJde, dJde_form)
    dJde.assemble()
    
    B = dolfinx.fem.petsc.assemble_vector(dJde_form)
    B.assemble()
    
    dRde.mult(lamda.vector, B)
    dJde += B
    dJde.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

   
    return dJde

# START ITERATION ---------------------------------
mn = 0.5
loop= 100
J_eps_s = []
for i in range(loop):
    i= i+1

    #solve objective function
    J_eps(mn)

    #solve gradient
    grad_J(mn,dJde=dJde)

    #update eps
    eps.x.array[:] = mn
    eps_new = eps-10**17*dJde.array
    mn.vector()[:] = eps_new.vector()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[16], line 18
     16 eps.x.array[:] = mn
     17 eps_new = eps+10**17*dJde.array
---> 18 mn.vector()[:] = eps_new.vector()

AttributeError: 'numpy.ndarray' object has no attribute 'vector'

dokken · March 6, 2024, 3:42pm

Please note that you are mixing PETSc-Vectors, its local form (a view into memory) and general arrays (dolfinx.la.Vector) in your code.

An easy fix to your problem is simply calling

   eps.x.array[:] = mn
    eps_new = eps.vector.array-10**17*dJde.array
    mn.vector()[:] = eps_new

as you then extract the array in the first step.

Please note that making a reproducible (runnable) minimal example is key to getting accurate help.

Pasuta · March 8, 2024, 3:36pm

I tried to do as you suggested but this happens many times, it has [ * ] as in the picture like it becomes unresponsive and I couldn’t see the result if it’s worked or not . Do you have any idea about this??

mn = fem.Function(V)
mn.x.array[:] = 0.5 
loop= 100

for i in range(loop):
    i = i+1
    
    #solve objective function
    J_eps(mn)

    #solve gradient
    grad_J(mn,dJde=dJde)

    #update eps
    eps.x.array[:] = mn
    eps_new = eps.vector.array-10**17*dJde.array
    mn.vector()[:] = eps_new

dokken · March 9, 2024, 1:04pm

I cannot give you any further help, as you have not supplied a minimal reproducible example.

bilenemek · March 9, 2024, 3:51pm

Please have a look at this project for topology optimization, it is working pretty well!

http://www.dolfin-adjoint.org/en/latest/

Pasuta · March 11, 2024, 9:34am

please look at this when i put f (just randomly, try to see why code goes wrong) instead of eps, error is shown below

f.x.array[:] = mn

mn = fem.Function(V)
mn.x.array[:] = 0.5 
loop= 100

for i in range(loop):
    i = i+1

    #solve objective function
    #J_eps(mn)

    #solve gradient
    #grad_J(mn,dJde=dJde)

    #update eps
    f.x.array[:] = mn
    eps_new = eps.vector.array-10**17*dJde.array
    mn.vector()[:] = eps_new

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 15
      6 i = i+1
      8 #solve objective function
      9 #J_eps(mn)
     10 
   (...)
     13 
     14 #update eps
---> 15 f.x.array[:] = mn
     16 eps_new = eps.vector.array-10**17*dJde.array
     17 mn.vector()[:] = eps_new

NameError: name 'f' is not defined

but when i put eps.x.array[:] = mn there are no results or error occurred as i said before, it’s like something or code (maybe JupyterLab) doesn’t respond.

mn = fem.Function(V)
mn.x.array[:] = 0.5 
loop= 100

for i in range(loop):
    i = i+1

    #solve objective function
    #J_eps(mn)

    #solve gradient
    #grad_J(mn,dJde=dJde)

    #update eps
    eps.x.array[:] = mn
    eps_new = eps.vector.array-10**17*dJde.array
    mn.vector()[:] = eps_new

Pasuta · March 11, 2024, 1:59pm

it seems like the code isn’t available, I can’t preview

Pasuta · March 12, 2024, 8:20am

I tried to call this function to check if it’s worked but it has an error

def update(m_n):
    eps.x.array[:] = m_n
    eps_new = eps.vector.array-10**17*dJde.array
    mn.vector()[:] = eps_new
    return mn

update(0.1)
print(eps_new)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[25], line 1
----> 1 update(0.1)
      2 print(eps_new)

Cell In[23], line 4, in update(m_n)
      2 eps.x.array[:] = m_n
      3 eps_new = eps.vector.array-10**17*dJde.array
----> 4 mn.vector()[:] = eps_new
      5 return mn

AttributeError: 'float' object has no attribute 'vector'

Ekrem_Ekici · March 12, 2024, 11:57am

As long as you do not provide a MWE, I am afraid that nobody will be in a position to help you.

Pasuta · March 13, 2024, 10:54am

@dokken
I’m really sorry for that, here i tried again and something went wrong at def update_eps , I’ve been stuck here for many days and i can’t update the eps until now. Could you please help me.

and here is the research i try following https://www.sciencedirect.com/science/article/pii/S0017931022011930

import dolfinx
import ufl
from dolfinx.plot import vtk_mesh
import pyvista
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative
from dolfinx.mesh import CellType, create_rectangle
from dolfinx import fem, nls, io
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import time
from matplotlib import*
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner
import dolfinx.nls.petsc
import dolfinx.fem.petsc
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx import default_scalar_type
from dolfinx import *

import numpy as np
from mpi4py import MPI
import dolfinx
from ufl import grad, dot
import os
from dolfinx.cpp.mesh import CellType
from dolfinx.io import XDMFFile
import dolfinx.nls.petsc
import dolfinx.fem.petsc


#from dolfinx import has_patsc_complex
#Parameter
D0 = 1;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;


# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [20, 20], CellType.triangle)
V = fem.FunctionSpace(mesh, ("CG", 1))
c = fem.Function(V)
phi = TestFunction(V)

# The eps is the design variable
eps = fem.Function(V)
eps.x.array[:] = 0.8

dx = Measure("dx") ##integration over different parts of the domain
x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(0)) ##Neauman = 0

# Residual of the variational form of Poisson problem
R = inner(D0*eps**eta*grad(c), grad(phi))*dx - inner(k0*(1-eps)**gamma*c,phi)*dx - inner(gu,phi) *ds     

##Dirichlet 
def on_boundary(x):
    return np.isclose(x[0], 0)
dofs_L = fem.locate_dofs_geometrical(V , on_boundary)
bc_L = fem.dirichletbc(default_scalar_type(1), dofs_L, V) ##Dirichlet = 1
bcs = [bc_L]

#Forming and solving the linear system
problem_u = fem.petsc.NonlinearProblem(R, c, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem_u)
solver.rtol = 1e-16
solver.solve(c)

#how big is the solution
J = ufl.inner(k0*(1-eps)**gamma,c)*dx ##Objective function
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c) 


bc = []
#Adjective problem
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
#solve adj problem in lamda (lamda is adjoint variable)
lamda = problem_adj.solve()

# dJde is objactive function / eps
dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))
dJde = dolfinx.fem.petsc.assemble_vector(dJde_form)
dJde.assemble()

# dRde is the equation accounting for the transport of the species / eps
dRde_form = fem.form(adjoint(derivative(R,eps)))
dRde = dolfinx.fem.petsc.assemble_matrix(dRde_form)
dRde.assemble()


B = dolfinx.fem.petsc.assemble_vector(dJde_form)
B.assemble()

dRde.mult(lamda.vector, B)

# Update dJde
dJde += B
dJde.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Function to compute the objective
def J_eps(m_n):
    eps = m_n.x.array[:] 
    solver.solve(c)
    J_old = fem.assemble_scalar(J_form)
    return J_old

# Function to compute the gradient of the reaction rate w.r.t. eps
def grad_J(m_n,dJde=dJde):
    eps = m_n.x.array[:]
    
    # Direct solution
    solver.solve(c)

    # Adjoint solution
    lamda = problem_adj.solve()

    # Reassemble dRde
    dRde.zeroEntries()
    fem.petsc.assemble_matrix(dRde, dRde_form)
    dRde.assemble()

    # Compute gradient
    # zero gradient first
    with dJde.localForm() as dJde_loc:
        dJde_loc.set(0)

    fem.petsc.assemble_vector(dJde, dJde_form)
    dJde.assemble()
    
    B = dolfinx.fem.petsc.assemble_vector(dJde_form)
    B.assemble()
    
    dRde.mult(lamda.vector, B)
    dJde += B
    dJde.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    return dJde

def update_eps(m_n,result_grad_J):
    eps = m_n.x.array[:]
    eps_new = fem.Function(V)
    eps_new.x.array[:] = eps.x.array - 10**17*result_grad_J.x.array   
    m_n.x.array[:] = eps_new  
    return m_n

loop = 10
mn = fem.Function(V)
mn.x.array[:] = 0.5
result_grad_J = fem.Function(V)

for i in range(loop):

    #solve obj
    J_eps(mn)
    print("J_eps: ",J_eps(mn))
    
    #solve gradient
    grad_J(mn,dJde=dJde)
    result_grad_J.x.array[:] = grad_J(mn,dJde=dJde)
    print("grad_J: ",result_grad_J.x.array)

    #update eps
    update_eps(mn,result_grad_J)
    mn = update_eps(mn, result_grad_J)
    print("eps: ",mn.x.array.tolist())

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 18
     14 result_grad_J.x.array[:] = grad_J(mn,dJde=dJde)
     15 #print("grad_J: ",result_grad_J.x.array)
     16 
     17 #update eps
---> 18 update_eps(mn,result_grad_J)
     19 mn = update_eps(mn, result_grad_J)
     20 print("eps: ",mn.x.array.tolist())

Cell In[3], line 42, in update_eps(m_n, result_grad_J)
     40 eps = m_n.x.array[:]
     41 eps_new = fem.Function(V)
---> 42 eps_new.x.array[:] = eps.x.array - 10**17*result_grad_J.x.array   
     43 m_n.x.array[:] = eps_new  
     44 return m_n

AttributeError: 'numpy.ndarray' object has no attribute 'x'

dokken · March 13, 2024, 3:24pm

Dear @Pasuta,
I will make an exception and reply to this post, even if I believe it does not qualify as a minimal example.
Here are a couple of reasons:

Unused imports
Unused variables

However, as the error is quite clear, consider

Pasuta:

def update_eps(m_n,result_grad_J):
    eps = m_n.x.array[:]
    eps_new = fem.Function(V)
    eps_new.x.array[:] = eps.x.array - 10**17*result_grad_J.x.array   
    m_n.x.array[:] = eps_new  
    return m_n

Here, you first create a variable eps, which is an numpy.ndarray. You then try to access x.array of this variable in

Just remove x.array, i.e.

def update_eps(m_n,result_grad_J):
    eps = m_n.x.array[:]
    eps_new = fem.Function(V)
    eps_new.x.array[:] = eps - 10**17*result_grad_J.x.array   
    m_n.x.array[:] = eps_new  
    return m_n

and the code does contine.

Please note that I made an adjoint solver for DOLFINx at: Calculate sensitivity when the control is a fem function - #3 by dokken
which is a general wrapper for linear problems

Pasuta · March 14, 2024, 7:53am

Thank you for your answer, I will remove unused imports. I tried doing as you said but the results of eps didn’t come out , because i don’t understand much about looping or iteration so i don’t know why are the results not as expected and it seems like it can only be computed just 1 iteration or do you have any example about looping which similar to my problem.

def update_eps(m_n,result_grad_J):
    eps = m_n.x.array[:]
    eps_new = fem.Function(V)
    eps_new.x.array[:] = eps - 10**17*result_grad_J.x.array   
    m_n.x.array[:] = eps_new  
    return m_n

loop = 10
mn = fem.Function(V)
mn.x.array[:] = 0.1
result_grad_J = fem.Function(V)

for i in range(loop):

    #solve obj
    J_eps(mn)
    print("J_eps: ",J_eps(mn))
    
    #solve gradient
    grad_J(mn,dJde=dJde)
    result_grad_J.x.array[:] = grad_J(mn,dJde=dJde)
  
    #update eps
    update_eps(mn,result_grad_J)
    #mn = update_eps(mn, result_grad_J)
    print("eps: ",update_eps.x.array)

the result

J_eps:  (-0.03918748536752729+0j)

Topic		Replies	Views
Try to evaluate the sensitivity value but couldn't solve the correct value mathematics dolfinx	3	97	June 23, 2024
A question arose during the process of converting existing 55-line code of topology optimization to dolfinx mesh	5	984	June 19, 2023
Unstructured Mesh - Topological Optimization mesh mesh	3	89	August 5, 2024
Parallelisation of topology optimisation dolfinx	2	1068	June 13, 2022
Interpolation issue for topology optimisation dolfinx	2	415	June 8, 2022

Topology optimization loop cannot update the new design variable

Related topics