Hi,
I want to create a periodic micromodel similar to this tutorial but with a hyperelastic material. I’m using dolfinx_mpc, and I am trying to adapt this test case with a nonlinear MPC problem to my application.
This is giving me an error related to the NonlinearMPCProblem class. What error am I making in my problem definition?
Here is the MWE (the classes are unchanged from the example):
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, io, mesh
# from dolfinx.la.petsc import create_vector
from dolfinx.fem.petsc import create_vector
import dolfinx_mpc
import dolfinx.fem.petsc
import dolfinx.nls.petsc
import petsc4py.PETSc as PETSc
from ufl import ( Identity, grad, ln, tr, det, variable, derivative, TestFunction, TrialFunction )
"""
Solving functions copied from dolfinx_mpc/python/tests/test_nonlinear_assembly.py
"""
class NonlinearMPCProblem(dolfinx.fem.petsc.NonlinearProblem):
def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_options={}, jit_options={}):
self.mpc = mpc
super().__init__(F, u, bcs=bcs, J=J, form_compiler_options=form_compiler_options, jit_options=jit_options)
def F(self, x: PETSc.Vec, F: PETSc.Vec): # type: ignore
with F.localForm() as F_local:
F_local.set(0.0)
dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F)
# Apply boundary condition
dolfinx_mpc.apply_lifting(
F,
[self._a],
bcs=[self.bcs],
constraint=self.mpc,
x0=[x],
scale=dolfinx.default_scalar_type(-1.0),
)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore
dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0)
def J(self, x: PETSc.Vec, A: PETSc.Mat): # type: ignore
A.zeroEntries()
dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A)
A.assemble()
class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
def __init__(
self,
comm: MPI.Intracomm,
problem: NonlinearMPCProblem,
mpc: dolfinx_mpc.MultiPointConstraint,
):
"""A Newton solver for non-linear MPC problems."""
super().__init__(comm)
self.mpc = mpc
self.u_mpc = dolfinx.fem.Function(mpc.function_space)
# Create matrix and vector to be used for assembly of the non-linear
# MPC problem
self._A = dolfinx_mpc.cpp.mpc.create_matrix(problem.a._cpp_object, mpc._cpp_object)
self._b = create_vector(mpc.function_space.dofmap.index_map, mpc.function_space.dofmap.index_map_bs)
self.setF(problem.F, self._b)
self.setJ(problem.J, self._A)
self.set_form(problem.form)
self.set_update(self.update)
def update(self, solver: dolfinx.nls.petsc.NewtonSolver, dx: PETSc.Vec, x: PETSc.Vec): # type: ignore
# We need to use a vector created on the MPC's space to update ghosts
self.u_mpc.x.petsc_vec.array = x.array_r
self.u_mpc.x.petsc_vec.axpy(-1.0, dx)
self.u_mpc.x.petsc_vec.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD, # type: ignore
) # type: ignore
self.mpc.homogenize(self.u_mpc)
self.mpc.backsubstitution(self.u_mpc)
x.array = self.u_mpc.x.petsc_vec.array_r
x.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD, # type: ignore
) # type: ignore
def solve(self, u: dolfinx.fem.Function):
"""Solve non-linear problem into function u. Returns the number
of iterations and if the solver converged."""
n, converged = super().solve(u.x.petsc_vec)
u.x.scatter_forward()
return n, converged
@property
def A(self) -> PETSc.Mat: # type: ignore
"""Jacobian matrix"""
return self._A
@property
def b(self) -> PETSc.Vec: # type: ignore
"""Residual vector"""
return self._b
"""
Mesh
"""
length_plate = 1
# Create a simple 2D grid mesh (1x1) with 20 elements in each direction
domain = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=((0, 0), (length_plate, length_plate)),
n=(20, 20),
cell_type=mesh.CellType.triangle
)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
dim = domain.topology.dim
"""
Function spaces
"""
degree = 1
shape = (dim,)
V = fem.functionspace(domain, ("P", degree, shape))
"""
Boundary conditions
"""
def bot_left_corner(x):
return np.logical_and(np.isclose(x[1], 0.0), np.isclose(x[0], 0.0))
def right_boundary(x):
return np.isclose(x[0], length_plate)
def top_boundary(x):
return np.isclose(x[1], length_plate)
bot_left_dof = fem.locate_dofs_geometrical(
V, bot_left_corner
)
top_facets = mesh.locate_entities_boundary(domain, dim-1, top_boundary)
top_dofs = fem.locate_dofs_topological(V.sub(1), dim-1, top_facets)
u_top = fem.Constant(domain, PETSc.ScalarType(1.))
# Fixed bottom left corner, prescribed value on all top dofs
bcs = [fem.dirichletbc(np.zeros((dim,)), bot_left_dof, V), fem.dirichletbc(u_top, top_dofs, V.sub(1))]
a = np.array([1, 0])
def periodic_relation_left_right(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0] - 1.0
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
def periodic_relation_bottom_top(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0]
out_x[1] = x[1] - 1.0
out_x[2] = x[2]
return out_x
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, right_boundary, periodic_relation_left_right, bcs)
mpc.create_periodic_constraint_geometrical(V, top_boundary, periodic_relation_bottom_top, bcs)
mpc.finalize()
# Define variational problem
V_mpc = mpc.function_space
# unkown lives in MPC space
u = fem.Function(V_mpc, name="Displacement")
"""
Hyperelastic material model (compressible neo-Hookean model)
"""
# Identity tensor
Id = Identity(dim)
F = variable(Id + grad(u)) # Deformation gradient
C = F.T * F # Right Cauchy-Green tensor
# Invariants of deformation tensors
I1 = tr(C)
J = det(F)
E = 1e3
nu = 0.3
mu = fem.Constant(domain, E / 2 / (1 + nu))
lmbda = fem.Constant(domain, E * nu / (1 - 2 * nu) / (1 + nu))
# Stored strain energy density
psi = mu / 2 * (I1 - 3 - 2 * ln(J)) + lmbda / 2 * (J - 1) ** 2
"""
Problem setup
"""
dx = ufl.Measure("dx", domain=domain, metadata={"quadrature_degree": 2})
E_pot = psi * dx
v = TestFunction(V)
du = TrialFunction(V)
Residual = derivative(
E_pot, u, v
) # This is equivalent to Residual = inner(P, grad(v))*dx
Jacobian = derivative(Residual, u, du)
problem = NonlinearMPCProblem(Residual, u, mpc, bcs=bcs, J=J)
solver = NewtonSolverMPC(mesh.comm, problem, mpc)
solver.atol = 1e-4
solver.rtol = 1e-4
"""
Solving of problem
"""
out_file = "compression_hypel_MWE.xdmf"
with io.XDMFFile(domain.comm, out_file, "w") as xdmf:
xdmf.write_mesh(domain)
for n, disp in enumerate(np.linspace(0, -0.2, 11)[1:]):
u_top.value = disp
num_its, converged = solver.solve(u)
assert converged
u.x.scatter_forward() # updates ghost values for parallel computations
print(
f"Time step {n}, Number of iterations {num_its}, Disp {disp:.3f}"
)
with io.XDMFFile(domain.comm, out_file, "a") as xdmf:
xdmf.write_function(u, n + 1)
Snippets from the error:
in <module> problem = NonlinearMPCProblem(Residual, u, mpc, bcs=bcs, J=J)
super().__init__(F, u, bcs=bcs, J=J, form_compiler_options=form_compiler_options, jit_options=jit_options)
~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
self._a = _create_form(
~~~~~~~~~~~~^
J, form_compiler_options=form_compiler_options, jit_options=jit_options
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
NotImplementedError: Cannot take length of non-vector expression.
Perhaps related is the fact that I can’t import from dolfinx.la.petsc import create_vector, but instead use from dolfinx.fem.petsc import create_vector.
I appreciate any suggestions