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:
[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
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)