# Try to evaluate the sensitivity value but couldn't solve the correct value
**URL:** https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223
**Category:** mathematics
**Tags:** dolfinx
**Created:** [April 11, 2024, 4:31pm UTC](https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223 "2024-04-11T16:31:50Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@Pasuta](https://fenicsproject.discourse.group/u/Pasuta)
#### Post date: [April 11, 2024, 4:31pm UTC](https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223/1 "2024-04-11T16:31:50Z")
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
.
```auto
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
```auto
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)
#solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c)
def on_boundary_1(x):
return np.isclose(x[0], 0)
dofs_L_adj = fem.locate_dofs_geometrical(V , on_boundary_1)
bc_L_adj = fem.dirichletbc(default_scalar_type(0), dofs_L_adj, V) ##Dirichlet = 1
bcs_adjoint = [bc_L_adj]
#Adjoint problem
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"})
#solve adj problem in lamda (lamda is adjoint variable)
lamda = problem_adj.solve()
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)
grad_c = grad(c)
grad_lamda = grad(lamda)
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)
print("dJde :",dJde)
```
The results are
```auto
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](https://fenicsproject.discourse.group/t/calculate-sensitivity-equation-of-reaction-diffusion-systems-with-dolfinx/13545/20)) 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))))
```auto
# 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 = dJde.copy()
dRde.mult(lamda.vector, B)
# Update dJde
dJde += B
dJde.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
```
I will be very pleased. If someone can help solve this problem.
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [April 12, 2024, 11:45am UTC](https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223/2 "2024-04-12T11:45:34Z")
> [@Pasuta](#):
>
> `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)`
Here you are mixing numpy arrays (`.x.array` with ufl objects.
> [@Pasuta](#):
>
> 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](https://fenicsproject.discourse.group/t/calculate-sensitivity-equation-of-reaction-diffusion-systems-with-dolfinx/13545/20)) 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))))
>
> ```auto
>
> ```
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
---
### Author:  [@Pasuta](https://fenicsproject.discourse.group/u/Pasuta)
#### Post date: [April 16, 2024, 3:09pm UTC](https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223/3 "2024-04-16T15:09:58Z")
As you can see the equation shown in the picture ,the first term equation is the same term as dJde in this line
```auto
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)
```
```auto
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 = (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])+(k0gamma_((1-eps.x.array)\*\*(eta-1))_c.x.array_lamda.x.array)
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [June 23, 2024, 7:13pm UTC](https://fenicsproject.discourse.group/t/try-to-evaluate-the-sensitivity-value-but-couldnt-solve-the-correct-value/14223/4 "2024-06-23T19:13:05Z")
> [@Pasuta](#):
>
> 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.
I don’t understand why the maximum value should be -1. Could you please elaborate on this?
> [@Pasuta](#):
>
> `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)`
`dJde = -1 * 2 * (0.5)*c`
with your current settings, so the first term should be `dJde=-c`?
If you are talking about:
> [@Pasuta](#):
>
> ```auto
> 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()
>
> ```
then this is a very different derivative than the first one, as it is an integral over the domain, which is not the same as the values at the dofs.