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
bcs_adjoint = [bc_L]
lamda = fem.Function(V)
phi_2 = TestFunction(V)
R_2 = inner(D0*eps**eta*grad(lamda), grad(phi_2))*dx - inner(k0*(1-eps)**gamma*lamda,phi_2)*dx-inner(k0*(1-eps)**gamma,phi_2)*dx
adjoint_u = fem.petsc.NonlinearProblem(R_2, lamda, bcs_adjoint)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,adjoint_u)
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_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)
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-----------------------------------------------