How to set Neumann boundary conditions to solve adjoint variable

Hello, I want to find the sensitivity from the objective function but things went wrong when i solved the adjoint variable, the result was not as i expected. I want to set Neumann boundary conditions to be 0 on all sides of the design space, but I don’t know what’s wrong with the results. How can i fix this??

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
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 dolfinx import default_scalar_type
from petsc4py import PETSc

#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]], [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.5

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) 


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

#solve adj problem in lamda (lamda is adjoint variable)
lamda = problem_adj.solve()  #adjoint variable

print("The objective is",J_old)

lamda_ = lamda.x.array
min_ = min(lamda_)
max_ = max(lamda_)
print("min :",min_)
print("max :",max_)
The objective is (0.0788579733134757+0j)
min : (-1.0000000000000107+0j)
max : (-1.0000000000000013+0j)

but the actual result should be

The objective is: 0.0787966
 min -0.915644  max -1.31291e-34

@dokken Sorry for disturbing you at this time but I do not have any idea to solve this problem , Could you please suggest me or need more information, I really need to correct this result.

Dear @Pasuta

It is quite hard to give you any guidance when you do not specify mathematically what problems you are solving.

I have already sketched out how to set up the adjoint equations in DOLFINx at: Calculate sensitivity when the control is a fem function - #3 by dokken
which I have sent to you before.

I have explained the mathematical equation, the objective function and adjoint equation i want to solve in this link Calculate sensitivity equation of reaction-diffusion systems with dolfinx and I tried as your tutorial but had this error

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
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 dolfinx import default_scalar_type
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)
J = -ufl.inner((k0*(1-eps)**gamma),c)*dx ##Objective function
J_form = fem.form(J)

dRdc = adjoint(derivative(R, c))
lhs = fem.form(dRdc)
rhs = derivative(J, c)

lmbda = fem.Function(V)
lambda_0 = fem.Function(V)
lambda_0.x.array[:] = 0
bcs_adjoint = [fem.dirichletbc(lambda_0, bc._cpp_object.dof_indices()[0]) for bc in bcs]  
adjoint_problem = dolfinx.fem.petsc.LinearProblem(lhs, rhs, c=lmbda,bcs=bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Exception ignored in: <function LinearProblem.__del__ at 0x7f334b39d510>
Traceback (most recent call last):
  File "/usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 627, in __del__
    self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[7], line 12
     10 lambda_0.x.array[:] = 0
     11 bcs_adjoint = [fem.dirichletbc(lambda_0, bc._cpp_object.dof_indices()[0]) for bc in bcs]  
---> 12 adjoint_problem = dolfinx.fem.petsc.LinearProblem(lhs, rhs, c=lmbda,bcs=bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

TypeError: LinearProblem.__init__() got an unexpected keyword argument 'c'

Well, it is quite clear from thw error message. There is no c kwarg in linear solver. Did you mean u=lmbda to specify the unknown function that you want to populate with a solution?

I already set c as fem.function just like you did, so I thought it should to be c=lmbda

c = fem.Function(V)
R = inner(D0eps**etagrad(c), grad(phi))dx - inner(k0(1-eps)**gamma*c,phi)*dx - inner(gu,phi) *ds

self.u = fem.Function(V)
R = inner(grad(self.u), grad(v))*dx - inner(self.f, v)*dx

To me it is unclear what you are trying to achieve.
I think what you want is

adjoint_problem = dolfinx.fem.petsc.LinearProblem(lhs, rhs, u=lmbda,bcs=bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

see dolfinx.fem.petsc — DOLFINx 0.7.3 documentation