Calculate sensitivity equation of reaction-diffusion systems with dolfinx

General
, ,
20

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-----------------------------------------------