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