Compute the sensitivity in parallel with DOLFINx

dolfinx
2

The problem with your code is that you start using scipy matrices to do the operations. If you continue using ‘pure’ dolfinx/ufl/petsc operations, everything works as expected.
Consider the minimal code below:

import dolfinx
import ufl
from dolfinx.generation import RectangleMesh
from dolfinx.mesh import CellType
from dolfinx.fem import VectorFunctionSpace, FunctionSpace, Function, DirichletBC
from dolfinx.nls import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

# mesh and function space
mesh = RectangleMesh(MPI.COMM_WORLD, [[0, 0, 0], [1, 1, 0]],
                     [3, 3], CellType.quadrilateral)
V = VectorFunctionSpace(mesh, ("CG", 1))
S = FunctionSpace(mesh, ("DG", 0))
u = Function(V)  # test function (displacment field)
v = ufl.TestFunction(V)  # trial function

# kinematics
lambda_ = 0.5769
mu = Function(S)  # in this case, we set the mu as a design variable
mu.vector.array = 0.3846
P = 2.0*mu*ufl.sym(ufl.grad(u)) + lambda_ * \
    ufl.tr(ufl.sym(ufl.grad(u)))*ufl.Identity(len(u))
dx = ufl.Measure("dx", metadata={"quadrature_degree": 4})
R = ufl.inner(ufl.grad(v), P)*dx

# boundary conditions


def left(x): return np.isclose(x[0], 0)
def right(x): return np.isclose(x[0], 1)


fdim = mesh.topology.dim - 1

u1 = Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
with u1.vector.localForm() as loc:
    loc.set(0)
bc1 = DirichletBC(u1, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
u2 = Function(V)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
with u2.vector.localForm() as loc:
    loc.set(1)
bc2 = DirichletBC(u2, dolfinx.fem.locate_dofs_topological(
    V, fdim, boundary_facets))
bcs = [bc1, bc2]

# solve the problem
problem = dolfinx.fem.NonlinearProblem(R, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)
# print(u.vector.array)

# set an arbitrary functinal
J = mu*ufl.dot(u, u)*dx
J_old = dolfinx.fem.assemble_scalar(J)

# partial derivative of J w.r.t. mu
dJdmu = dolfinx.fem.assemble_vector(ufl.derivative(J, mu))
dJdmu.assemble()

# partial derivative of R w.r.t. mu
dRdmu = dolfinx.fem.assemble_matrix(
    ufl.adjoint(ufl.derivative(R, mu)))  # partial derivative
dRdmu.assemble()

# reset the boundary condition
with u2.vector.localForm() as loc:
    loc.set(0.0)

# solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R, u))
rhs = -ufl.derivative(J, u)
problem = dolfinx.fem.LinearProblem(lhs, rhs, bcs=bcs)
lmbda = problem.solve()

dJdmu += dRdmu*lmbda.vector
dJdmu.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print(dJdmu.array)
4 Likes