# Mixed elements crashes with apply\_lifting (Bus Error)

**URL:** <https://fenicsproject.discourse.group/t/mixed-elements-crashes-with-apply-lifting-bus-error/18027>\
**Category:** General\
**Tags:** boundary-conditions, dolfinx\
**Created:** [July 25, 2025, 6:24pm UTC](https://fenicsproject.discourse.group/t/mixed-elements-crashes-with-apply-lifting-bus-error/18027 "2025-07-25T18:24:02Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![a-lyrae](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/a-lyrae/32/5878_2.png) [@a-lyrae](https://fenicsproject.discourse.group/u/a-lyrae)\
**Post date:** [July 25, 2025, 6:24pm UTC](https://fenicsproject.discourse.group/t/mixed-elements-crashes-with-apply-lifting-bus-error/18027/1 "2025-07-25T18:24:02Z")

</div>

Hello,

I’m solving a EHD problem where I need to solve two nonlinear equations simultaneously. For that, I implemented the problem by using mixed elements and needed to write a custom nonlinear solver. This crashes as soon as I use apply\_lifting with the following error:

```auto
[0]PETSC ERROR: Caught signal number 10 BUS: Bus Error, possibly illegal memory access
  Proc: [[56095,0],0]
  Errorcode: 59

```

To the best of my abilities, I couldn’t figure out the problem. Could you please help me? I’m using dolfinx 0.9.0.  
thanks, below a MWE

```auto

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh, fem
from dolfinx.fem import (dirichletbc, locate_dofs_topological, Constant)
from dolfinx.cpp import fem as cpp_fem
from dolfinx.nls.petsc import NewtonSolver
from basix.ufl import mixed_element

# Create unit square mesh
comm = MPI.COMM_WORLD
domain = mesh.create_unit_square(comm, 8, 8, mesh.CellType.triangle)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)

# Define continuous Lagrange space of degree 1
V = fem.functionspace(domain, ("Lagrange", 1))

# Create mixed space W = V x V
P1 = V.ufl_element()
mixed_el = mixed_element([P1, P1])
W = fem.functionspace(domain, mixed_el)

# Split mixed solution and test functions
u = fem.Function(W)
(p, w) = ufl.split(u)
(q, v) = ufl.TestFunctions(W)

# Build a simple nonlinear problem resembling your structure
# For demonstration we keep it simple:
# Solve: p + w^3 = f1, and p^2 + w = f2 on the domain
f1 = Constant(domain, 1.0)
f2 = Constant(domain, 2.0)

F = (p + w**3 - f1) * q * ufl.dx + (p**2 + w - f2) * v * ufl.dx

# Nonlinear problem class
class SimpleNonlinearProblem(fem.petsc.NonlinearProblem):
    def __init__ (self, F, u, bcs):
        super(). __init__ (F, u, bcs)
        self.u = u
        self.bcs = bcs
        self.F_form = F
        self.mesh = domain

    def F(self, b, x):
        # Copy current iterate vector into solution
        self.u.x.petsc_vec.copy(x)

        # Assemble residual
        formF = fem.form(self.F_form)
        with b.localForm() as bl:
            bl.set(0.0)
        fem.petsc.assemble_vector(b, formF)
        
        # Apply lifting
        fem.petsc.apply_lifting(
            b,
            [formF],
            bcs=[self.bcs],
            x0=[self.u.x.petsc_vec],
            alpha=-1.0
        )

        # Enforce Dirichlet BC explicitly
        fem.petsc.set_bc(b, self.bcs, x, alpha=-1.0)

        return 0

    def J(self, A, x):
        # Copy current iterate vector into solution
        self.u.x.petsc_vec.copy(x)

        # Compute Jacobian form
        J_form = ufl.derivative(self.F_form, self.u)
        formJ = fem.form(J_form)

        fem.petsc.assemble_matrix(A, formJ, self.bcs)
        A.assemble()

        return 0

# Define boundary facets (all four edges in this simple example)
def boundary_all(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

bdry_dofs_0 = locate_dofs_topological(W.sub(0), domain.topology.dim - 1,
                                      mesh.exterior_facet_indices(domain.topology))
bdry_dofs_1 = locate_dofs_topological(W.sub(1), domain.topology.dim - 1,
                                      mesh.exterior_facet_indices(domain.topology))

# Create Dirichlet BCs on each subspace
bc_zero_0 = dirichletbc(Constant(domain, 0.0), bdry_dofs_0, W.sub(0))
bc_zero_1 = dirichletbc(Constant(domain, 0.0), bdry_dofs_1, W.sub(1))

bcs = [bc_zero_0, bc_zero_1]

# Create the nonlinear problem and Newton solver
problem = SimpleNonlinearProblem(F, u, bcs)
solver = NewtonSolver(comm, problem)
solver.rtol = 1e-7
solver.atol = 1e-9
solver.max_it = 20
solver.monitor = True
solver.line_search = "basic"

# Initial guess zero
u.x.array[:] = 0.0
u.x.scatter_forward()

# Solve nonlinear system
num_iter, converged = solver.solve(u)
if comm.rank == 0:
    if converged:
        print(f"Newton solver converged in {num_iter} iterations.")
    else:
        print(f"Newton solver did NOT converge in {num_iter} iterations.")

# Access solution components
p_sol = u.sub(0).collapse()
w_sol = u.sub(1).collapse()

if comm.rank == 0:
    print("Solution p at dofs:", p_sol.x.array)
    print("Solution w at dofs:", w_sol.x.array)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 30, 2025, 10:04am UTC](https://fenicsproject.discourse.group/t/mixed-elements-crashes-with-apply-lifting-bus-error/18027/2 "2025-07-30T10:04:25Z")

</div>

The problem is that the lifting doesn’t take the linear-form (residual) into lifting, but the jacobian.

This is closer to what you’d like:

```auto
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import (dirichletbc, locate_dofs_topological, Constant)
from dolfinx.nls.petsc import NewtonSolver
from basix.ufl import mixed_element

# Create unit square mesh
comm = MPI.COMM_WORLD
domain = mesh.create_unit_square(comm, 8, 8, mesh.CellType.triangle)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)

# Define continuous Lagrange space of degree 1
V = fem.functionspace(domain, ("Lagrange", 1))

# Create mixed space W = V x V
P1 = V.ufl_element()
mixed_el = mixed_element([P1, P1])
W = fem.functionspace(domain, mixed_el)

# Split mixed solution and test functions
u = fem.Function(W)
(p, w) = ufl.split(u)
(q, v) = ufl.TestFunctions(W)

# Build a simple nonlinear problem resembling your structure
# For demonstration we keep it simple:
# Solve: p + w^3 = f1, and p^2 + w = f2 on the domain
f1 = Constant(domain, 1.0)
f2 = Constant(domain, 2.0)

F = (p + w**3 - f1) * q * ufl.dx + (p**2 + w - f2) * v * ufl.dx

# Nonlinear problem class
class SimpleNonlinearProblem(fem.petsc.NonlinearProblem):
    def __init__ (self, F, u, bcs):
        super(). __init__ (F, u, bcs)
        self.u = u
        self.bcs = bcs
        self.mesh = domain

    def F(self, b, x):
        # Copy current iterate vector into solution
        self.u.x.petsc_vec.copy(x)

        # Assemble residual
        with b.localForm() as bl:
            bl.set(0.0)
        fem.petsc.assemble_vector(b, self.L)
        
        # Apply lifting
        fem.petsc.apply_lifting(
            b,
            [self.a],
            bcs=[self.bcs],
            x0=[self.u.x.petsc_vec],
            alpha=-1.0
        )

        # Enforce Dirichlet BC explicitly
        fem.petsc.set_bc(b, self.bcs, x, alpha=-1.0)

        return 0

    def J(self, A, x):
        # Copy current iterate vector into solution
        self.u.x.petsc_vec.copy(x)

        # Compute Jacobian form
        fem.petsc.assemble_matrix(A, self.a, self.bcs)
        A.assemble()

        return 0

# Define boundary facets (all four edges in this simple example)
def boundary_all(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

bdry_dofs_0 = locate_dofs_topological(W.sub(0), domain.topology.dim - 1,
                                      mesh.exterior_facet_indices(domain.topology))
bdry_dofs_1 = locate_dofs_topological(W.sub(1), domain.topology.dim - 1,
                                      mesh.exterior_facet_indices(domain.topology))

# Create Dirichlet BCs on each subspace
bc_zero_0 = dirichletbc(Constant(domain, 0.0), bdry_dofs_0, W.sub(0))
bc_zero_1 = dirichletbc(Constant(domain, 0.0), bdry_dofs_1, W.sub(1))

bcs = [bc_zero_0, bc_zero_1]

# Create the nonlinear problem and Newton solver
problem = SimpleNonlinearProblem(F, u, bcs)
solver = NewtonSolver(comm, problem)
solver.rtol = 1e-7
solver.atol = 1e-9
solver.max_it = 20
solver.monitor = True
solver.line_search = "basic"

# Initial guess zero
u.x.array[:] = 0.0
u.x.scatter_forward()

# Solve nonlinear system
num_iter, converged = solver.solve(u)
if comm.rank == 0:
    if converged:
        print(f"Newton solver converged in {num_iter} iterations.")
    else:
        print(f"Newton solver did NOT converge in {num_iter} iterations.")

# Access solution components
p_sol = u.sub(0).collapse()
w_sol = u.sub(1).collapse()

if comm.rank == 0:
    print("Solution p at dofs:", p_sol.x.array)
    print("Solution w at dofs:", w_sol.x.array)

```

All-though I don’t understand what you want a custom newton solver for