Optimal control - inverse problem in hyperelasticity

dolfin-adjoint
1

Hello everyone,

I am trying to adapt the example “Optimal control in DOLFINx interfacing with scipy” (Optimal control in DOLFINx interfacing with scipy — FEniCS Workshop) to solve an inverse problem in hyperelasticity.

I am having the following issue: ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Exp to expression depending on form argument v_0.

I know where it does come from, but have no clue how to fix it, since it is the expression of the nonlinear forward problem.

The error appears when I am trying to create the adjoint problem:
adj_problem = LinearSolver(ufl.replace(dFdu_adj, {uh: v}), -dJdu, lmbda, bcs)

Below is a piece of code (I cleaned a little bit):

domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 8, 8, mesh.CellType.hexahedron)

# -----------------------------------------------------------------------------
ud = prob_hiperelastico(domain)
# -----------------------------------------------------------------------------

V = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
W = dolfinx.fem.functionspace(domain, ("Lagrange", 1,))

du = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
uh = dolfinx.fem.Function(V)

# -----------------------------------------------------------------------------
# Equacao de Estado (Fun)
# -----------------------------------------------------------------------------

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

def right(x):
    return np.isclose(x[0], 1)
    
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

# Kinematics
# Spatial dimension
d = len(uh)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(uh))
J = ufl.variable(ufl.det(F))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
E = 0.5*(C - I)           

f0 = as_vector([ 1.0, 0.0, 0.0 ])
s0 = as_vector([ 0.0, 1.0, 0.0 ])
n0 = as_vector([ 0.0, 0.0, 1.0 ])

# -----------------------------------------------------------------------------
# Constitutive model
# ----------------------------------------------------------------------------- 

def ctrue_expr(x):
    return 4.0 - 3.8 * ((x[0] - 0.7)**2 + (x[1] - 0.5)**2 + (x[2] - 0.5)**2 < 0.3**2)

Va = fem.functionspace(domain, ("Lagrange", 1))  
cj_true = fem.Function(Va)
cj_true.interpolate(ctrue_expr)

# C = default_scalar_type(2.0)
CC = cj_true
bf = default_scalar_type(8.0)
bt = default_scalar_type(2.0)
bfs = default_scalar_type(4.0)
e1 = f0
e2 = s0
e3 = n0
kappa = 1e2
kappa = fem.Constant(domain,kappa)

Cs = J**(-2/3) * F.T * F  
Es = 0.5 * (Cs - I)      

E11 = ufl.inner(Es*e1, e1)

Q = bf*E11**2 

# energia
Wpassive = CC/2.0 * (ufl.exp(Q) - 1)
Wactive  = 0.0
Wincompr = kappa * (J*ufl.ln(J) - J + 1)

strain_energy = Wpassive + Wactive + Wincompr
P = ufl.diff(strain_energy, F)


p_right = fem.Constant(domain,0.0) 

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

# Definition of the weak form:
N = ufl.FacetNormal(domain)
Gext = -p_right * ufl.inner(v, det(F) * ufl.inv(F) * N) * ds(2) #ds(2) = right boundary

# forward prob
Fun = ufl.inner(P, ufl.grad(v)) * dx + Gext 

# -------------------------------------------------------------------------------

# explorar
alpha = dolfinx.fem.Constant(domain, 1.0e-10)

x = ufl.SpatialCoordinate(domain)
J =  (1/2) * ufl.inner(uh - ud, uh - ud) * ufl.dx + (alpha/2) * ufl.inner(ufl.grad(CC), ufl.grad(CC)) * ufl.dx


forward_problem = NonLinearSolver(Fun, uh=uh, bcs=bcs)

lmbda = dolfinx.fem.Function(V)
dFdu = ufl.derivative(Fun, uh, du)
dFdu_adj = ufl.adjoint(dFdu)
dJdu = ufl.derivative(J, uh, v)


adj_problem = LinearSolver(
    ufl.replace(dFdu_adj, {uh: v}),
    -dJdu,
    lmbda,
    bcs
)

I hope this can make the question clear.
Any help is appreciated.

Best regards,

Bernardo M. Rocha

2

Your code is not reproducible,as you have not:

  1. Included all import statements
  2. have undefined variables, such as

If you change the line

to

adj_problem = LinearSolver(dFdu_adj, -dJdu, lmbda, bcs)

the code should run with no problem.

I cannot recall why I added the replace function there in the first place. Will have to revisit that..

3

Thanks a lot. Indeed it worked nicely now.

Sorry for not posting the entire working code, it was too messy.

Anyway, thanks a lot.