Calculate sensitivity equation of reaction-diffusion systems with dolfinx

General
, ,
13

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)

lhs = ufl.adjoint(ufl.derivative(R,c))
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"})
lamda = problem_adj.solve()

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_form = fem.form(adjoint(derivative(R,eps)))
dRde = fem.petsc.create_matrix(dRde_form)
dRde.assemble()

#Gradient
# 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'