NonlinearMPCProblem: error with periodic structure and nonlinear material

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

First of all, you should change your constraint to:

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) & np.invert(right_boundary(x))

to avoid constraining the same dof twice.

Secondly, you are passing in the wrong value for the jacobian:

This should be

problem = NonlinearMPCProblem(Residual, u, mpc, bcs=bcs, J=Jacobian)
solver = NewtonSolverMPC(domain.comm, problem, mpc)
solver.atol = 1e-4
solver.rtol = 1e-4

Additionally, your change of create_vector is not correct. What version of DOLFINx are you using?

1 Like

Thank you dokken! That fixed the error I had before - now I indeed get a mistake for the create_vector function.
I’m using dolfinx 0.9.0 and dolfinx_mpc 0.9.1.

If you are using v0.9.1, you need to look at the demo for the 0.9.1 version:dolfinx_mpc/python/tests/test_nonlinear_assembly.py at v0.9.1 · jorgensd/dolfinx_mpc · GitHub

1 Like

I see, thank you for that, the create_vector function changed. I’ve updated that part of the code.

Now I don’t get any direct errors, but the Newton solver still can not converge (for any step, even with lenient tolerance settings). If I set the loading to 0 then I do “converge” and the code runs.
Am I still making a mistake in my boundary conditions?

Here is an updated code with the previously discussed changes:

import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, io, mesh
import dolfinx_mpc
import dolfinx.la as _la
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 = _la.create_petsc_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) & np.invert(right_boundary(x))

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=Jacobian)
solver = NewtonSolverMPC(domain.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)