# Try to evaluate the sensitivity value but couldn't solve the correct value

I am now solving a topology optimization problem that requires solving for sensitivity by substituting concentration, eps as the design variable and lamda into this equation
.

``````dJde = (k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array) + (D0*eta*(eps.x.array**(eta-1)))*(grad_c[0]*grad_lamda[0]+grad_c[1]*grad_lamda[1])+(k0*gamma*((1-eps.x.array)**(eta-1))*c.x.array*lamda.x.array)
``````

But the answer results are as below. which seems impossible to calculate

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

#Parameter
D0 = 0.1;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;

# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [100, 100], 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.5

#boundary_markers = FacetFunction('size_t',mesh)
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     # " - inner(gu,phi) *ds " is boundary neauman

##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)

rhs = -ufl.derivative(J,c)

def on_boundary_1(x):
return np.isclose(x[0], 0)

print("The objective is",J_old)

lamda_ = lamda.x.array
min_ = min(lamda_)
max_ = max(lamda_)
print("lamda min :",min_)
print("lamda max :",max_)

# Update dJde (sensitivity)
dJde = fem.Function(V)

print("dJde :",dJde)

``````

The results are

``````The objective is (0.07877750212780599+0j)
lamda min : (-0.9155314560137311+0j)
lamda max : 0j
dJde : [Sum(FloatValue(0.07740840145414743), Sum(FloatValue(-0.08455854979804014), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
Sum(FloatValue(0.07737913038121018), Sum(FloatValue(-0.08452332241286431), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
Sum(FloatValue(0.07737327161484303), Sum(FloatValue(-0.08451627181655852), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
...
Sum(FloatValue(-1.0), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),)))))))
Sum(FloatValue(0.030059484951120978), Sum(FloatValue(-0.9689781605244571), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),))))))))
Sum(FloatValue(-1.0), Product(FloatValue(0.1), Sum(Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(0),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(0),)))), Product(Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 80)), MultiIndex((FixedIndex(1),))), Indexed(Grad(Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False), (2,)), 20), Basix element (P, triangle, 1, gll_warped, unset, False)), 82)), MultiIndex((FixedIndex(1),)))))))]
``````

I had also tried this method(You can see more explanations at this link Calculate sensitivity equation of reaction-diffusion systems with dolfinx - #20 by Pasuta) but the answer was also wrong and I think there is something wrong with this line.

dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))

``````# 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 = dolfinx.fem.petsc.assemble_matrix(dRde_form)
dRde.assemble()

B = dJde.copy()
dRde.mult(lamda.vector, B)

# Update dJde
dJde += B
``````

I will be very pleased. If someone can help solve this problem.

Here you are mixing numpy arrays (`.x.array` with ufl objects.

Please elaborate on what you imply by «wrong». This gives you the symbolic expression for dJ/d(eps), which you in turn can assemble into a vector

As you can see the equation shown in the picture ,the first term equation is the same term as dJde in this line

``````dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))
dJde1 = dolfinx.fem.petsc.assemble_vector(dJde_form)
dJde1.assemble()
print("dJde1:",dJde1.array)
``````
``````dJde1: [-4.22792750e-06+0.j -1.40883944e-06+0.j -4.22651812e-06+0.j ...
-4.94829826e-05+0.j -4.87114826e-05+0.j -1.65374090e-05+0.j]
``````

If I try substituting the values ​​into the equation directly (in first term only), the variable is defined as follows and the maximum value of c is 1, so the maximum value of djde1 should be -1, But the maximum value printed is -1.65374090e-05 So I think that line of code is incorrect.

#Parameter
D0 = 0.1
k0 = -1.0
eta = 2.0
gamma = 2.0
eps = 0.5

And for this, I tried so many ways but still couldn’t find the correct way, can you correct it?? and I don’t know what is different with numpy array and ufl objects.

`dJde = -1 * 2 * (0.5)*c`
with your current settings, so the first term should be `dJde=-c`?