# Calculate sensitivity equation of reaction-diffusion systems with dolfinx

First method, I try to solve my problem by following (1) Calculate sensitivity when the control is a fem function
The objective function is, which epsilon is the design variable

, the adjoint equation is

and the sensitivity

First, i solve the concentration in the function space V and the distribution seems good

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 *
#from dolfinx import has_patsc_complex
#Parameter
theta0 = 0.5;
theta = theta0;
D0 = 1.0;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;
Rfl = 0.0018;
URF = 0.1;

# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [20, 20], CellType.triangle)
#domainU = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(mesh, ("CG", 1))

c = fem.Function(V)
phi = TestFunction(V)
eps = fem.Function(V) # The eps is the design variable
eps.x.array[:] = 0.5

dx = Measure("dx")

x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(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)
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
ch = solver.solve(c)

but in the second part of solving adjoint equation has an error ,when I apply this step that the link I following above is error

#Define the objective function
J = ufl.inner(k0*(1-eps)**gamma,c)*dx
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#Bilinear and linear form of the adjoint problem
rhs = -ufl.derivative(J,c)

bc = []
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

#compute sensitivity
dJde_form = fem.form(ufl.derivative(J, eps))

dJde = fem.petsc.assemble_vector(dJde_form)
dJde.assemble()

# partial derivative of R w.r.t. eps
dRde_form = fem.form(adjoint(derivative(R, eps)))
dRde = fem.petsc.assemble_matrix(dRde_form)
dRde.assemble()

# Calculate the gradient of J with respect to eps (the epsilon)
dJde += dRde*lambda.vector

error has occured at this line #compute sensitivity

dJde_form = fem.form(ufl.derivative(J, eps))

"TimeoutError: JIT compilation timed out, probably due to a failed previous compile.Try cleaning cache (e.g. remove /root/.cache/fenics/libffcx_forms_5eb899f22d9b05f26c47e06ab454b4561b19ac5a.c) or increase timeout option."

then i try following from this link How to increase the timeout parameter in JIT compilation in FEniCSX? but still has error.

Second method, I try to substitute variable values instead because i already have the sensitivity equation as this pic

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

print(dj_0)

the result is, which i don’t know what does it mean

<dolfinx.fem.function.Expression object at 0x7f52f46480a0>

Now i try 2 methods

1. Following the first link Calculate sensitivity when the control is a fem function
2. Substituting variable values in the sensitivity equation directly

My questions are

1. Which methods should i follow?
2. If i choose the first method how could i figure out the error?
3. The result when i substitute variable values in the sensitivity equation, what does it mean?

I tried copying and pasting method 1, and my script runs without issues. You may want to rename the variable lambda to lambda_, because the former is reserved by python as a keyword to define lambda functions.

Please try again deleting your cache (typically in \$HOME/.cache/fenicsx)

how can i delete cache?

It is a directory in your home, I think the default is \$HOME/.cache/fenicsx. Just locate it (in a terminal, or in whatever GUI you feel comfortable) and delete the directory.

Now I can fixed timeout error, but i have another issue

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

from dolfinx.cpp.mesh import CellType
from dolfinx.io import XDMFFile

#from dolfinx import has_patsc_complex
#Parameter
theta0 = 0.5;
theta = theta0;
D0 = 1.0;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;
Rfl = 0.0018;
URF = 0.1;

# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [20, 20], CellType.triangle)
#domainU = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(mesh, ("CG", 1))

c = fem.Function(V)
phi = TestFunction(V)
eps = fem.Function(V) # The eps is the design variable
eps.x.array[:] = 0.5

dx = Measure("dx")

x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(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)
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
ch = solver.solve(c)

#Define the objective function
J = ufl.inner(k0*(1-eps)**gamma,c)*dx
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#Bilinear and linear form of the adjoint problem
rhs = -ufl.derivative(J,c)

lamda = fem.Function(V)
bc = []
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

#compute sensitivity
dJde_form = dolfinx.fem.form(ufl.derivative(J, eps, ufl.conj(ufl.TestFunction(eps.function_space))))
dJde = fem.petsc.create_vector(dJde_form)
dJde.assemble()

# partial derivative of R w.r.t. eps
dRde = fem.petsc.create_matrix(dRde_form)
dRde.assemble()

# Calculate the gradient of J with respect to eps (the epsilon)
dJde += dRde*lamda.vector
print(dJde.array)

I have error in this line

dJde += dRde*lamda.vector
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 dJde += dRde*lamda.vector
3 print(dJde.array)

File petsc4py/PETSc/Mat.pyx:395, in petsc4py.PETSc.Mat.__mul__()

File petsc4py/PETSc/petscmat.pxi:642, in petsc4py.PETSc.mat_mul()

TypeError: Cannot convert petsc4py.PETSc.Vec to petsc4py.PETSc.Mat

Could you help me

Use dRde.mult, see for instance; oasisx/src/oasisx/fracstep.py at main · ComputationalPhysiology/oasisx · GitHub

I tried doing as you mentioned but the error kept showing, i’m sorry but those examples are very difficult to understand and I don’t understand how can my equation give a solution that isn’t a real number??

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 dJde += dRde.mult*lamda.vector

File petsc4py/PETSc/Vec.pyx:93, in petsc4py.PETSc.Vec.__rmul__()

File petsc4py/PETSc/petscvec.pxi:325, in petsc4py.PETSc.vec_rmul()

File petsc4py/PETSc/petscvec.pxi:309, in petsc4py.PETSc.vec_mul()

File petsc4py/PETSc/petscvec.pxi:284, in petsc4py.PETSc.vec_imul()

File petsc4py/PETSc/PETSc.pyx:150, in petsc4py.PETSc.asScalar()

TypeError: must be real number, not builtin_function_or_method

Please look at the syntax in the referenced link:

self._A.mult(self._u1[i].vector, self._wrk_comp.vector)

which multiplies _A by _u1 and places the result in _wrk_comp

another question, do i have to do like this from your code?? Calculate sensitivity when the control is a fem function - #3 by dokken

class Problem():
def __init__(self, V, bcs, u_obs):
"""Initialize all variational forms for the given problem, for a given set of boundary conditions and observed u"""
# Residual of the variational form of Poisson problem
self.f = fem.Function(V)  # The forcing f is the design variable
self.u = fem.Function(V)
R = inner(grad(self.u), grad(v))*dx - inner(self.f, v)*dx

because when i did this

self._A.mult(self._u1[i].vector, self._wrk_comp.vector)

it said

NameError: name 'self' is not defined

Without a minimal example it is hard to guide you. I use self as I am accessing a variable from my class.

That’s my problem because there is no minimal example so i may ask for your help if you could do anything to my code . this part is the hardest for me. i’ve been stuck here for weeks as you can see. I know you’re very busy so I appreciate any time you can give me

The problem is that you provide a code, I suggest a fix, and then you only provide a minimal snippet that it is unclear where fits into your code, thus I cannot determine where things go wrong for you, as I cannot reproduce your error message. Like, you do not have a class Problem in your previous code, so why would you have that in the new code?
The point of my previous reply is that if you want to compute b=Ax, you do A.mult(x.vector, b.vector)

I’m really sorry for my questions, actually my confusing is what are _A and _u1, i supposed to ask you that , I want to compute dJde += dRde*lamda.vector which dJde is B, dRde is A and lamda is x?? as you said b = Ax so i add this dRde.mult(lamda.vector, dJde.vector)

dRde.mult(lamda.vector, dJde.vector)

# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [20, 20], CellType.triangle)
#domainU = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(mesh, ("CG", 1))

c = fem.Function(V)
phi = TestFunction(V)
eps = fem.Function(V) # The eps is the design variable
eps.x.array[:] = 0.5

dx = Measure("dx")

x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(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)
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
ch = solver.solve(c)

J = ufl.inner(k0*(1-eps)**gamma,c)*dx
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

rhs = -ufl.derivative(J,c)
lamda = fem.Function(V)
bc = []
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

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

dRde = fem.petsc.create_matrix(dRde_form)
dRde.assemble()

# b=Ax, dJde = b, dRde = A, lamda = x
#dJde += dRde*lamda.vector
#A.mult(x.vector, b.vector)
dRde.mult(lamda.vector, dJde.vector)

and i got error

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[84], line 4
1 # b=Ax, dJde = b, dRde = A, lamda = x
2 #dJde += dRde*lamda.vector
3 #A.mult(x.vector, b.vector)
----> 4 dRde.mult(lamda.vector, dJde.vector)

AttributeError: 'petsc4py.PETSc.Vec' object has no attribute 'vector'

I’m so sorry to bother you many times , but can you help me fix this error again. Thank you.

dJde is already a vector, so just use dJde instead of dJde.vector

Thank you so much for your help, but i don’t know why my solution is like this

dRde.mult(lamda.vector, dJde)
print(dJde.array)

Did i do something wrong?? Can you suggest me

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
0.+0.j]

dJde is the zero vector, so that result is not surprising. If you want it to be nonzero you should use dolfinx.fem.petsc.assemble_vector to assemble any form you’d like.

from this equation dJde += dRde*lambda.vector

I apply b=Ax → A.mult(x.vector, b.vector)

dRde.mult(lamda.vector, B)

# Update dJde
dJde += B

but it has error ‘B’ is not defined, How can i define ‘B’

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[17], line 112
92 dRde.assemble()
95 # b=Ax, dJde = b, dRde = A, lamda = x
96 #dJde += dRde*lamda.vector
97 #A.mult(x.vector, b.vector)
(...)
107 #B = fem.Function(V)
108 #B.x.array[:] = 0
--> 112 dRde.mult(lamda.vector, B)
114 # Update dJde
115 dJde += B

NameError: name 'B' is not defined

You need to define it before use it, simply:

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

I tried as you said but the result is still wrong ,why it has to define B as dJde

B = dJde.copy()

This is my lastest code , I have check the result (The objective, concentration(c), lamda) all of these are correct except the dJde(sensitivity) and I think something is not correct at #update dJde, I don’t know how to fix this.

import dolfinx
import numpy as np
import dolfinx.nls.petsc
import ufl
import dolfinx.fem.petsc
import pyvista
import scipy
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)
ds = Measure('ds')
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]

problem_u = fem.petsc.NonlinearProblem(R, c, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem_u)
solver.rtol = 1e-16
solver.solve(c)

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

lamda = fem.Function(V)
phi_2 = TestFunction(V)

solver.rtol = 1e-16
solver.solve(lamda)

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

# 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

print("The objective is :",J_old)
c_ = c.x.array
min_ = min(c_)
max_ = max(c_)
print(" ")
print("c min :",min_, "c max :",max_)
lamda_ = lamda.x.array
min_ = min(lamda_)
max_ = max(lamda_)
print(" ")
print("lamda min :",min_,"lamda max :",max_)

dJde_ = dJde.array
min_ = min(dJde_)
max_ = max(dJde_)
print(" ")
print("dJde min :",min_,"dJde max :",max_)

The result is

The objective is : (0.07877750212780599+0j)

c min : (0.08446854398718186+0j) c max : (1+0j)

lamda min : (-0.9155314560137318+0j) lamda max : 0j

dJde min : (-8.858601822635755e-07+0j) dJde max : (-1.1906820423691968e-07+0j)

but the actual result from original code(C++) is

The objective is: 0.0787966
-- C :
min 0.0843556  max 1
-- lamda :
min -0.915644  max -1.31291e-34

sensitivity(dJde) min : -0.00712961 sensitivity max : -0.00676388-----------------------------------------------