Calculate sensitivity equation of reaction-diffusion systems with dolfinx

First method, I try to solve my problem by following (1) Calculate sensitivity when the control is a fem function
The objective function is, which epsilon is the design variable


, the adjoint equation is

adj2
and the sensitivity
adj1

First, i solve the concentration in the function space V and the distribution seems good

import dolfinx
import ufl
from dolfinx.plot import vtk_mesh
import pyvista
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative
from dolfinx.mesh import CellType, create_rectangle
from dolfinx import fem, nls, io
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import time
from matplotlib import*
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner
import dolfinx.nls.petsc
import dolfinx.fem.petsc
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx import default_scalar_type
from dolfinx import *
#from dolfinx import has_patsc_complex
#Parameter
theta0 = 0.5;
theta = theta0;
D0 = 1.0;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;
Rfl = 0.0018;
URF = 0.1;

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

but in the second part of solving adjoint equation has an error ,when I apply this step that the link I following above is error

#Define the objective function
J = ufl.inner(k0*(1-eps)**gamma,c)*dx
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#Bilinear and linear form of the adjoint problem
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c) 

#Solve adjoint equation
bc = []
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
lambda = problem_adj.solve()

#compute sensitivity
dJde_form = fem.form(ufl.derivative(J, eps))

dJde = fem.petsc.assemble_vector(dJde_form)
dJde.assemble()

# partial derivative of R w.r.t. eps
dRde_form = fem.form(adjoint(derivative(R, eps)))
dRde = fem.petsc.assemble_matrix(dRde_form)
dRde.assemble()

# Calculate the gradient of J with respect to eps (the epsilon)
dJde += dRde*lambda.vector
dJde.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

error has occured at this line #compute sensitivity

dJde_form = fem.form(ufl.derivative(J, eps))

"TimeoutError: JIT compilation timed out, probably due to a failed previous compile.Try cleaning cache (e.g. remove /root/.cache/fenics/libffcx_forms_5eb899f22d9b05f26c47e06ab454b4561b19ac5a.c) or increase timeout option." 

then i try following from this link How to increase the timeout parameter in JIT compilation in FEniCSX? but still has error.

Second method, I try to substitute variable values instead because i already have the sensitivity equation as this pic

dj_0 = fem.Function(V) # The eps is the design variable
dj_0.x.array[:] = 0

dj_0 = dolfinx.fem.Expression((k0*gamma*(1-eps)**(gamma-1)*c)+((D0*eta*eps**(eta-1))*dot(grad(c),grad(lamda)))+(k0*gamma*c*lamda*(1-eps)**(gamma-1)),V.element.interpolation_points())
print(dj_0)

the result is, which i don’t know what does it mean

<dolfinx.fem.function.Expression object at 0x7f52f46480a0>

Now i try 2 methods

  1. Following the first link Calculate sensitivity when the control is a fem function
  2. Substituting variable values in the sensitivity equation directly

My questions are

  1. Which methods should i follow?
  2. If i choose the first method how could i figure out the error?
  3. The result when i substitute variable values in the sensitivity equation, what does it mean?

I tried copying and pasting method 1, and my script runs without issues. You may want to rename the variable lambda to lambda_, because the former is reserved by python as a keyword to define lambda functions.

Please try again deleting your cache (typically in $HOME/.cache/fenicsx)

how can i delete cache? :smiling_face_with_tear:

It is a directory in your home, I think the default is $HOME/.cache/fenicsx. Just locate it (in a terminal, or in whatever GUI you feel comfortable) and delete the directory.

Now I can fixed timeout error, but i have another issue


import dolfinx
import ufl
from dolfinx.plot import vtk_mesh
import pyvista
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative
from dolfinx.mesh import CellType, create_rectangle
from dolfinx import fem, nls, io
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import time
from matplotlib import*
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner
import dolfinx.nls.petsc
import dolfinx.fem.petsc
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx import default_scalar_type
from dolfinx import *

import numpy as np
from mpi4py import MPI
import dolfinx
from ufl import grad, dot

from dolfinx.cpp.mesh import CellType
from dolfinx.io import XDMFFile


#from dolfinx import has_patsc_complex
#Parameter
theta0 = 0.5;
theta = theta0;
D0 = 1.0;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;
Rfl = 0.0018;
URF = 0.1;

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

#Define the objective function
J = ufl.inner(k0*(1-eps)**gamma,c)*dx
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)

#Bilinear and linear form of the adjoint problem
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c) 

#Solve adjoint equation
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()

#compute sensitivity
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()

# partial derivative of R w.r.t. eps
dRde_form = fem.form(adjoint(derivative(R,eps)))
dRde = fem.petsc.create_matrix(dRde_form)
dRde.assemble()

# Calculate the gradient of J with respect to eps (the epsilon)
dJde += dRde*lamda.vector
dJde.ghostUpdate(addv=PETSc.InserMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print(dJde.array)

I have error in this line

dJde += dRde*lamda.vector
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 dJde += dRde*lamda.vector
      2 dJde.ghostUpdate(addv=PETSc.InserMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
      3 print(dJde.array)

File petsc4py/PETSc/Mat.pyx:395, in petsc4py.PETSc.Mat.__mul__()

File petsc4py/PETSc/petscmat.pxi:642, in petsc4py.PETSc.mat_mul()

TypeError: Cannot convert petsc4py.PETSc.Vec to petsc4py.PETSc.Mat

Could you help me :sob: :sob:

Use dRde.mult, see for instance; oasisx/src/oasisx/fracstep.py at main · ComputationalPhysiology/oasisx · GitHub

I tried doing as you mentioned but the error kept showing, i’m sorry but those examples are very difficult to understand and I don’t understand how can my equation give a solution that isn’t a real number??

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 dJde += dRde.mult*lamda.vector

File petsc4py/PETSc/Vec.pyx:93, in petsc4py.PETSc.Vec.__rmul__()

File petsc4py/PETSc/petscvec.pxi:325, in petsc4py.PETSc.vec_rmul()

File petsc4py/PETSc/petscvec.pxi:309, in petsc4py.PETSc.vec_mul()

File petsc4py/PETSc/petscvec.pxi:284, in petsc4py.PETSc.vec_imul()

File petsc4py/PETSc/PETSc.pyx:150, in petsc4py.PETSc.asScalar()

TypeError: must be real number, not builtin_function_or_method

Please look at the syntax in the referenced link:

self._A.mult(self._u1[i].vector, self._wrk_comp.vector)

which multiplies _A by _u1 and places the result in _wrk_comp

another question, do i have to do like this from your code?? Calculate sensitivity when the control is a fem function - #3 by dokken

class Problem():
    def __init__(self, V, bcs, u_obs):
        """Initialize all variational forms for the given problem, for a given set of boundary conditions and observed u"""
        # Residual of the variational form of Poisson problem
        self.f = fem.Function(V)  # The forcing f is the design variable
        self.u = fem.Function(V)
        R = inner(grad(self.u), grad(v))*dx - inner(self.f, v)*dx

because when i did this

self._A.mult(self._u1[i].vector, self._wrk_comp.vector)

it said

NameError: name 'self' is not defined

Without a minimal example it is hard to guide you. I use self as I am accessing a variable from my class.