Derivative w.r.t. a quantity defined over a real_functionspace: Problem in parallel

dolfin-adjoint
,
1

Dear community/developers,

I am using FEniCSx 0.9 with Docker on Microsoft Windows. I am trying to code a constitutive parameter identification routine with the help of the adjoint method for gradient computation.

In the MWE posted below, you can see that I am interested in finding the derivative DresFE_Dθ (which will be used later to compute the adjoint sensitivity). resFE is simply the residual in the finite element problem and θ is the gathering of “nf” constitutive parameters (here in anisotropic elasticity, a vector with nf = 6 components). I could not find any better way to define θ than using scifem’s create_real_functionspace. I did not use dolfinx.fem.functionspace because I do not want to associate “nf” degrees of freedom to all nodes, but rather I want to have only a vector of size “nf” for the whole domain. Here is the MWE:

from ufl import TrialFunction, TestFunction, dot, inner, grad, derivative, adjoint, sym, as_tensor, as_matrix, as_vector, action, dx
from dolfinx import mesh, default_scalar_type
from dolfinx.fem import Function, functionspace, form
from dolfinx.fem.petsc import assemble_matrix
from scifem import create_real_functionspace
from mpi4py import MPI
import numpy as np

# PARALLEL MPI COMMUNICATOR
comm = MPI.COMM_WORLD
rank = comm.Get_rank()  # Get process rank (ID)
nProcess = comm.Get_size()  # Get total number of processes
print(f"Process {rank} out of {nProcess} process(es) is running.")

# MESH GENERATION
aspectRatio = 2.00
elementSize = 5
Lx = 100.0 # Domain dimension X
Ly = Lx*aspectRatio  # Domain dimension Y
Nx = np.floor(Lx/elementSize).astype(int)  # Number of elements X
Ny = np.floor(Ly/elementSize).astype(int) # Number of elements Y
domain = mesh.create_rectangle(comm, [(0.0, 0.0), (Lx, Ly)], [Nx, Ny], cell_type=mesh.CellType.quadrilateral)
gdim = domain.geometry.dim  # Geometrical dimension
element_order = 1
element_type = "Lagrange"
V = functionspace(domain, (element_type, element_order, (gdim, ))) # Displacement space

# ELASTICITY MATRIX IN VOIGT NOTATION
def elasticity_matrix_voigt(θ):
    """ Computes the 2D elasticity matrix in Voigt notation from material parameters. """
    return as_matrix([
        [θ[0], θ[1], θ[2]],
        [θ[1], θ[3], θ[4]],
        [θ[2], θ[4], θ[5]]])

# FINITE ELEMENT FORMULATION
u = TrialFunction(V)
v = TestFunction(V)

def strain(u_):
    """Computes the infinitesimal strain tensor"""
    return sym(grad(u_))

def stress(θ_, u_):
    """Computes the stress tensor in linear elasticity"""
    strain_tensor = strain(u_)
    strain_voigt = as_vector([strain_tensor[0,0], strain_tensor[1,1], 2*strain_tensor[0,1]])
    stress_voigt = dot(elasticity_matrix_voigt(θ_), strain_voigt)  # Voigt stress
    return as_tensor([[stress_voigt[0], stress_voigt[2]], [stress_voigt[2], stress_voigt[1]]])

# DEFINE THE RHS (I.E., LINEAR FORM): ZERO TRACTION & ZERO BODY FORCE
L = default_scalar_type(0.0)

# DEFINE THE FUNCTION SPACE FOR θ
nf = 6  # Number of material parameters (Voigt notation)
θ_V = create_real_functionspace(domain, (nf, ))
θ = Function(θ_V)

print(f"Process {rank}: θ = {θ.x.petsc_vec.getArray()}")

a = inner(stress(θ, u), strain(v)) * dx
U = Function(V)
resFE = action(a, U) - L

DresFE_Dθ = assemble_matrix(form(adjoint(derivative(resFE, θ))))
DresFE_Dθ.assemble()

local_nrows, local_ncols = DresFE_Dθ.getLocalSize()
print(f"Process {rank}: Local matrix size: {local_nrows} rows x {local_ncols} columns")

When I run the code using a single processor, I receive:

Process 0 out of 1 process(es) is running.
Process 0: θ = [0. 0. 0. 0. 0. 0.]
Process 0: Local matrix size: 6 rows x 1722 columns

and when I run the code in parallel with 2 processors, I receive:

Process 0 out of 2 process(es) is running.
Process 1 out of 2 process(es) is running.
Process 1: θ = []
Process 0: θ = [0. 0. 0. 0. 0. 0.]
Process 0: Local matrix size: 6 rows x 854 columns
Process 1: Local matrix size: 0 rows x 868 columns

As evident above, θ does not really exist on Process 1, and therefore, the derivative of resFE w.r.t. θ gets lost on Process 1, and so it causes wrong computations of the sensitivities (which are not shown here for brevity).

Any ideas and solutions to rectify this issue are very much appreciated.

2

getLocalSize does not include ghost entries, thus you won’t see the from the print.

I will illustrate this with a “simple” example:

from mpi4py import MPI
from dolfinx import mesh, fem
import dolfinx.fem.petsc
from scifem import create_real_functionspace
from ufl import TestFunction, TrialFunction, dx, inner
import numpy as np


domain = mesh.create_unit_square(
    MPI.COMM_WORLD,
    3,
    3,
    mesh.CellType.triangle,
    ghost_mode=dolfinx.mesh.GhostMode.shared_facet,
)
V = create_real_functionspace(domain, (3,))
Q = fem.functionspace(domain, ("Lagrange", 1, (3,)))
u = TrialFunction(V)
v = TestFunction(Q)
dx = dx(domain=domain)
a = u[0] * v[0] * dx + u[1] * v[1].dx(1) * dx + u[2] * v[0] * dx
A = fem.petsc.assemble_matrix(fem.form(a))
A.assemble()


def f(x):
    values = np.zeros((u.ufl_shape[0], x.shape[1]))
    values[0] = 1
    values[1] = 2
    values[2] = 4
    return values


b = dolfinx.fem.Function(V)
b.interpolate(f)



z = dolfinx.fem.Function(Q)
A.mult(b.x.petsc_vec, z.x.petsc_vec)
z.x.scatter_forward()

print(MPI.COMM_WORLD.rank, A.norm(0), dolfinx.la.norm(z.x))

here we have 3 components, and compute the matrix-vector product between the matrix and a corresponding vector.
You can inspect the number of owned entries, and number of ghosts with:


bs = V.dofmap.bs
print(
    MPI.COMM_WORLD.rank,
    b.x.array,
    bs * V.dofmap.index_map.size_local,
    bs * V.dofmap.index_map.num_ghosts,
)

which will show you all entries on the process ( split into local and ghosted entries)

2 Likes