Trouble solving mixed-domain PDE with BlockedNewtonSolver

Caglar · April 16, 2025, 10:29am

I'm trying to use the BlockedNewtonSolver to solve a coupled micromagnetic problem in a mixed domain on fenicsx v0.9.0 involving,

Magnetization, m,
Magnetic potential, phi,
Lagrange multiplier, l,

Magnetic potential phi is defined over the entire domain. The magnetization and the Lagrange multiplier are defined only on a subdomain (inner region). I succesfully run this formulation in a single domain without MixedFunctionSpace and BlockedNewtonSolver.

The residual is constructed with ufl.extract_blocks using entity_maps, and passed into BlockedNewtonSolver. However, the solver fails with the error below.

Any suggestions would be much appreciated.

Error

Traceback (most recent call last):
  File "test.py", line 101, in 
    solver = BlockedNewtonSolver(residual,[m, phi, l], bcs_phi)
  File "/opt/homebrew/anaconda3/envs/fenicsx-env/lib/python3.13/site-packages/scifem/solvers.py", line 312, in __init__
    J = ufl.extract_blocks(sum(ufl.derivative(sum(F), u[i], du[i]) for i in range(len(u))))
                           ~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/anaconda3/envs/fenicsx-env/lib/python3.13/site-packages/scifem/solvers.py", line 312, in 
    J = ufl.extract_blocks(sum(ufl.derivative(sum(F), u[i], du[i]) for i in range(len(u))))
                                              ~~~^^^
TypeError: unsupported operand type(s) for +: 'int' and 'Form'
Exception ignored in: 
Traceback (most recent call last):
  File "/opt/homebrew/anaconda3/envs/fenicsx-env/lib/python3.13/site-packages/scifem/solvers.py", line 361, in __del__
    self._J.destroy()
AttributeError: 'BlockedNewtonSolver' object has no attribute '_J'

Main Solver Code

import ufl
import numpy as np
from dolfinx import mesh, fem, default_scalar_type, io, la
from mpi4py import MPI
from dolfinx import io
from scifem import BlockedNewtonSolver
from math import pi, sqrt
from customNewtonSolver import *


Amp, LE, GE, ZE= 1.e0, 1.e0, 1.e0, 1.e0

mu0D = 4 * pi * 1.e-7 * (GE*LE)/(Amp**2*ZE**2)
mu1D = 1.25 * 1.e-1 * (GE*LE)/(Amp**2*ZE**2)
K1D = 5 * 1.e2 * GE/(ZE*ZE*LE)
AexD = 1.3 * 1.e-11 * (GE*LE)/(ZE*ZE)
MsD = 8 * 1.e5 * Amp / LE
g0D = 2.21 * 1.e5 * (Amp*ZE)/GE

#Dimensions
x1D = 60*1.e-9*LE
x2D = 3*1.e-9*LE
x3D = 30*1.e-9*LE

#Non-Dimensional Parameters
mu0S , AexS, K1S = 1, 1, 1
MsS = sqrt(mu0D / K1D) * MsD

msh , cell_markers, facet_markers = io.gmshio.read_from_msh("trapezoidal_prisms.msh", MPI.COMM_WORLD,0, gdim=3)

msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
cells_outer = cell_markers.find(2)
cells_inner = cell_markers.find(1)

inner_sm, msh_to_inner = mesh.create_submesh(msh, msh.topology.dim, cells_inner)[0:2]

dx = ufl.Measure("dx", domain=msh)
dx_inner = ufl.Measure("dx", domain=inner_sm)

cell_imap = msh.topology.index_map(msh.topology.dim)
num_cells = cell_imap.size_local + cell_imap.num_ghosts
omega_to_gamma = np.full(num_cells, -1)
omega_to_gamma[msh_to_inner] = np.arange(len(msh_to_inner))
entity_maps = {inner_sm: omega_to_gamma}

# Material properties
mu0 = fem.Constant(msh, default_scalar_type(mu0S))
Ms = fem.Constant(msh, default_scalar_type(MsS))
Aexc = fem.Constant(inner_sm, default_scalar_type(AexS))
K1 = fem.Constant(inner_sm, default_scalar_type(K1S))
kl = fem.Constant(inner_sm, default_scalar_type(1.e5))
alpha = fem.Constant(inner_sm, default_scalar_type(1.e0))

fs_m = fem.functionspace(inner_sm, ("Lagrange", 1, (msh.geometry.dim, )))
fs_l = fem.functionspace(inner_sm, ("Lagrange", 1))
fs_phi = fem.functionspace(msh, ("Lagrange", 1))
V = ufl.MixedFunctionSpace(fs_m, fs_phi, fs_l)

m = fem.Function(fs_m, name="Magnetization")
m0 = fem.Function(fs_m, name="Magnetization") #previous solution
phi = fem.Function(fs_phi, name="Magnetic Potential")
l = fem.Function(fs_l, name="Lagrange Multiplier")

dm, dphi, dl = ufl.TrialFunctions(V)
H = -ufl.grad(phi)
dH = -ufl.grad(dphi)

phi_var = ufl.variable(phi)
m_var = ufl.variable(m)
l_var = ufl.variable(l)
H_var = ufl.variable(H)

a=ufl.as_vector([1, 0, 0])

psi_mag = - 0.5 * mu0 * ufl.dot(H, H) - mu0 * Ms * ufl.dot(m, H)
psi_exc = Aexc * ufl.inner(ufl.grad(m), ufl.grad(m))
psi_ani = K1 * (1. - ufl.dot(m, a)**2)
psi_lag = l*(ufl.dot(m, m) - 1) - 0.5*l*l / kl
psi = (psi_mag + psi_exc + psi_ani + psi_lag)

B = -ufl.diff(psi, H_var)
T = ufl.diff(psi, m_var)

fdim = msh.topology.dim - 1
node = mesh.locate_entities(msh, 0, lambda x: np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0) & np.isclose(x[2], 0.0))

bc_phi1 = fem.dirichletbc(0.0, node, fs_m.sub(0))
bcs_phi = [bc_phi1]

t  = 0.
tmax = g0D * K1D * (1.e2 * ZE * 1.e-9) / (MsD * mu0D)
dt = tmax / 1.e7

F0 = ufl.dot(dH, B) * dx
F1 = alpha * ufl.dot(dm, T + (m-m0)/dt) * dx_inner  + alpha *  2 * Aexc * ufl.inner(ufl.grad(m), ufl.grad(dm)) * dx_inner
F2 = (ufl.dot(m,m) - 1 - l/kl) * dl * dx_inner + 2*l*ufl.dot(m, dm) * dx_inner

F = F0 + F1 + F2
residual = fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)

solver = BlockedNewtonSolver(residual,[m, phi, l], bcs_phi)

Gmsh Code

import gmsh
import sys
import pygmsh
from dolfinx.cpp.mesh import create_mesh
from mpi4py import MPI

gmsh.initialize()
gmsh.model.add("trapezoidal_prisms")

# Define the points
Lx = 0.37210420376762543  # Example length values
Ly = 0.01860521018838127
Lz = 0.18605210188381271

Ly_mid = Ly
Lz_mid = Lz
LMesh = 1.2*Ly
# Left prism points
p1 = gmsh.model.occ.addPoint(0, Ly-Ly, Lz-Lz, LMesh)
p2 = gmsh.model.occ.addPoint(0, Ly+Ly, Lz-Lz, LMesh)
p3 = gmsh.model.occ.addPoint(0, Ly+Ly, Lz+Lz, LMesh)
p4 = gmsh.model.occ.addPoint(0, Ly-Ly, Lz+Lz, LMesh)

# Right prism points
p9 = gmsh.model.occ.addPoint(Lx/2, Ly-Ly_mid, Lz-Lz_mid, LMesh)
p10 = gmsh.model.occ.addPoint(Lx/2,Ly+ Ly_mid,Lz -Lz_mid, LMesh)
p11 = gmsh.model.occ.addPoint(Lx/2,Ly+ Ly_mid,Lz+ Lz_mid, LMesh)
p12 = gmsh.model.occ.addPoint(Lx/2,Ly-Ly_mid, Lz+Lz_mid,LMesh)

p5 = gmsh.model.occ.addPoint(Lx,Ly -Ly,Lz -Lz, LMesh)
p6 = gmsh.model.occ.addPoint(Lx,Ly+ Ly,Lz -Lz, LMesh)
p7 = gmsh.model.occ.addPoint(Lx,Ly+ Ly,Lz+ Lz, LMesh)
p8 = gmsh.model.occ.addPoint(Lx,Ly -Ly,Lz+ Lz, LMesh)

lines = []
lines.append(gmsh.model.occ.add_line(p1, p2))
lines.append(gmsh.model.occ.add_line(p2, p3))
lines.append(gmsh.model.occ.add_line(p3, p4))
lines.append(gmsh.model.occ.add_line(p4, p1))
lines.append(gmsh.model.occ.add_line(p5, p6))
lines.append(gmsh.model.occ.add_line(p6, p7))
lines.append(gmsh.model.occ.add_line(p7, p8))
lines.append(gmsh.model.occ.add_line(p8, p5))
lines.append(gmsh.model.occ.add_line(p2, p10))
lines.append(gmsh.model.occ.add_line(p10, p11))
lines.append(gmsh.model.occ.add_line(p11, p3))
lines.append(gmsh.model.occ.add_line(p10, p6))
lines.append(gmsh.model.occ.add_line(p7, p11))
lines.append(gmsh.model.occ.add_line(p1, p9))
lines.append(gmsh.model.occ.add_line(p9, p12))
lines.append(gmsh.model.occ.add_line(p12, p4))
lines.append(gmsh.model.occ.add_line(p9, p5))
lines.append(gmsh.model.occ.add_line(p8, p12))
lines.append(gmsh.model.occ.add_line(p9, p10))
lines.append(gmsh.model.occ.add_line(p11, p12))

surfs = []
# Create surfaces and volumes
SN = 1
loop1 = gmsh.model.occ.addCurveLoop([lines[0], lines[1], lines[2], lines[3]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop1]))

loop2 = gmsh.model.occ.addCurveLoop([lines[4], lines[5], lines[6], lines[7]])
surfs.append( gmsh.model.occ.addPlaneSurface([loop2]))

loop3 = gmsh.model.occ.addCurveLoop([lines[8], lines[9], lines[10], -lines[1]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop3]))

loop4 = gmsh.model.occ.addCurveLoop([lines[11], lines[5], lines[12], -lines[9]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop4]))

loop5 = gmsh.model.occ.addCurveLoop([lines[13], lines[14], lines[15], lines[3]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop5]))

loop6 = gmsh.model.occ.addCurveLoop([lines[16], -lines[7], lines[17], -lines[14]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop6]))

loop7 = gmsh.model.occ.addCurveLoop([lines[8], -lines[18], -lines[13], lines[0]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop7]))

loop8 = gmsh.model.occ.addCurveLoop([lines[11], -lines[4], -lines[16], lines[18]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop8]))

loop9 = gmsh.model.occ.addCurveLoop([-lines[10], lines[19], lines[15], -lines[2]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop9]))

loop10 = gmsh.model.occ.addCurveLoop([-lines[12], lines[6], lines[17], -lines[19]])
surfs.append(gmsh.model.occ.addPlaneSurface([loop10]))
gmsh.model.occ.synchronize()

gmsh.model.occ.add_surface_loop([surfs[0], surfs[1], surfs[2], surfs[3], surfs[4], surfs[5], surfs[6], surfs[7], surfs[8], surfs[9]], 1)
gmsh.model.occ.synchronize()
vol = gmsh.model.occ.addVolume([1],1)
gmsh.model.occ.synchronize()

outer_box = gmsh.model.occ.addSphere(Lx/2, Ly, Lz,1*Lx,2)
gmsh.model.occ.synchronize()

entities, map_to_input = gmsh.model.occ.fragment([(3, 2)],[(3, 1)])
gmsh.model.occ.synchronize()

gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 8*LMesh)  # Finer mesh
gmsh.model.occ.synchronize()

left_surf = [idx for (dim, idx) in map_to_input[1] if dim == 3]
right_surf = [idx for (dim, idx) in map_to_input[0] if dim == 3 and idx not in left_surf]
gmsh.model.addPhysicalGroup(3, left_surf, tag=1)
gmsh.model.addPhysicalGroup(3, right_surf, tag=2)


gmsh.model.mesh.generate(3)
gmsh.write("trapezoidal_prisms.msh")

Caglar · April 18, 2025, 12:40pm

I was able to initialize the solver by using this dx definition,


marker = np.ones(num_cells, dtype=np.int32)
marker[cells_inner] = 2
cell_tag = mesh.meshtags(msh, 3, np.arange(num_cells, dtype=np.int32), marker)
dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tag)

F0 = ufl.dot(dH, B) * dx
F1 = alpha * ufl.dot(dm, T + (m-m0)/dt) * dx(2)  + alpha *  2 * Aexc * ufl.inner(ufl.grad(m), ufl.grad(dm)) * dx(2)
F2 = (ufl.dot(m,m) - 1 - l/kl) * dl * dx(2) + 2*l*ufl.dot(m, dm) * dx(2)

dokken · April 24, 2025, 8:42pm

Could you try not compiling the residual before passing it into the BlockedNewtonSolver, i.e.

solver = BlockedNewtonSolver(ufl.extract_blocks(F), [m, phi, l], bcs_phi)

Caglar · April 25, 2025, 9:00am

This didn’t help but I was able to solve the problem (at least initialized the solver without errors) using the “dx” definition that I shared in my last comment.

dokken · April 25, 2025, 11:44am

There is quite alot going on in your problem that concerns me.

One is your usage of ufl.variable (they are not used in your variational forms).

Secondly, you should use:
dm, dphi, dl = ufl.TestFunctions(V)
rather than

to define your residuals F0, F1, F2, as otherwise ufl.extract_blocks will not behave as expected.

Caglar · April 25, 2025, 1:36pm

Thank you for pointing those out. I made a simpler MWE as in posted https://fenicsproject.discourse.group/t/missing-diagonal-entry-when-solving-mixed-domain-pde/17469/3?u=caglar and made the following changes as you suggested,

du, dphi = ufl.TestFunctions(V)

eps = ufl.variable(ufl.sym(ufl.grad(u)))
H = ufl.variable(-ufl.grad(phi))

and got rid of “eps_var” and “H_var”,

u_, phi_ = ufl.TrialFunctions(V)
jac = ufl.derivative(F, u, u_) + ufl.derivative(F, phi, phi_)

However, I still get the following error,

File "petsc4py/PETSc/KSP.pyx", line 1782, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 73
[0] KSPSolve() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/ksp/ksp/interface/itfunc.c:1094
[0] KSPSolve_Private() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/ksp/ksp/interface/itfunc.c:427
[0] PCSetUp() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/ksp/pc/interface/precon.c:1086
[0] PCSetUp_ILU() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/ksp/pc/impls/factor/ilu/ilu.c:135
[0] MatILUFactorSymbolic() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/mat/interface/matrix.c:7094
[0] MatILUFactorSymbolic_SeqAIJ() at /Users/runner/miniforge3/conda-bld/bld/rattler-build_petsc_1738766455/work/src/mat/impls/aij/seq/aijfact.c:1738
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 0
Exception ignored in: <function _DeleteDummyThreadOnDel.__del__ at 0x104c06b60>
Traceback (most recent call last):
  File "/Applications/PyCharm.app/Contents/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_pep_669_tracing.py", line 561, in py_raise_callback
  File "/opt/homebrew/anaconda3/envs/fenicsx-env/lib/python3.13/threading.py", line 1435, in current_thread
TypeError: 'NoneType' object is not subscriptable

Caglar · April 28, 2025, 2:02pm

I apologize, this solves my problem.

Topic		Replies	Views
Missing diagonal entry when solving mixed-domain PDE dolfinx	4	70	April 28, 2025
BlockedNewtonSolver() and memory release problem General dolfinx	3	81	January 1, 2025
Variational foams in mixedfunctionspace General	1	46	December 11, 2024
Lagrange multipliers for Fictitious Domain Method General	6	185	June 19, 2024
Integrating on subdomain	30	2369	November 11, 2022

Trouble solving mixed-domain PDE with BlockedNewtonSolver

Error

Main Solver Code

Gmsh Code

Related topics