Missing diagonal entry when solving mixed-domain PDE

dolfinx
1

I'm trying to solve a coupled micromagnetic problem on two domains with fenicsx v0.9.0 that is installed with conda involving:

  • Magnetization (a 3D vector field),
  • Scalar magnetic potential, and
  • Lagrange multiplier to enforce unit norm of the magnetization.

These fields are defined in a mixed space, but only the scalar potential phi is defined over the entire domain. The magnetization and the Lagrange multiplier are defined only on a subdomain (inner region). The residual and jacobian are constructed with ufl.extract_blocks with entity_maps, and passed into the solver. However, I get the following error and quite unsure about what it means,

Any suggestions would be much appreciated.

Error

self._solver.solve(self.b, self.dx)
    ~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^
  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 2997

Main Code

import random
import ufl
import numpy as np
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem import Function, locate_dofs_geometrical, locate_dofs_topological
from mpi4py import MPI
from dolfinx import io
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

#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, inner_to_msh = mesh.create_submesh(msh, msh.topology.dim, cells_inner)[0:2]

# 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))

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

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)
dH = -ufl.grad(dphi)

m_var = ufl.variable(m)
l_var = ufl.variable(l)
H = -ufl.grad(phi)
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 = locate_dofs_geometrical(fs_phi, lambda x: np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0) & np.isclose(x[2], 0.0))
bc_phi1 = fem.dirichletbc(fem.Constant(inner_sm, default_scalar_type(1)), node, fs_phi)
bcs = [bc_phi1]

def initial_condition(x):
    values = np.zeros((3, x.shape[1]))
    for i in range(x.shape[1]):
        x_r, y_r, z_r = np.random.random(), np.random.random(), np.random.random()
        norm = np.sqrt(x_r**2 + y_r**2 + z_r**2) + 1e-12
        values[:, i] = np.array([x_r, y_r, z_r]) / norm
    return values
m.interpolate(initial_condition)
m.x.scatter_forward()

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

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)
F = F0 + F1 + F2
residual = fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)

m_, phi_, l_ = ufl.TestFunctions(V)
jac = ufl.derivative(F, m, m_) + ufl.derivative(F, phi, phi_) + ufl.derivative(F, l, l_)
J = fem.form(ufl.extract_blocks(jac), entity_maps=entity_maps)

solver = customNewtonSolver(residual, J, [m, phi, l], bcs)

tot_step = tmax/dt
counter, ii = 0, 0
m0.x.array[:] = m.x.array
while t < tmax:
    n, bConverged = solver.solve()

    m.x.scatter_forward()
    phi.x.scatter_forward()
    l.x.scatter_forward()

    if n==0:
        break

    # Update solution at previous time step (u_n)
    m0.x.array[:] = m.x.array
    t += dt

Solver


from dolfinx import  fem
from petsc4py import PETSc
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.cpp.la.petsc import scatter_local_vectors

class customNewtonSolver():
    max_iterations: int
    bcs: list[fem.DirichletBC]
    A: PETSc.Mat
    b: PETSc.Vec
    J: fem.Form
    b: fem.Form
    dx: PETSc.Vec

    def __init__(self, F: list[fem.form], J: list[list[fem.form]], w: list[fem.Function],
                 bcs: list[fem.DirichletBC] | None = None, max_iterations: int = 5,
                 petsc_options: dict[str, str | float | int | None] = None,
                 problem_prefix="newton"):
        self.max_iterations = max_iterations
        self.bcs = [] if bcs is None else bcs
        self.b = fem.petsc.create_vector_block(F)
        self.F = F
        self.J = J
        self.A = fem.petsc.create_matrix_block(J)
        self.dx = self.A.createVecLeft()
        self.w = w
        self.x = fem.petsc.create_vector_block(F)

        # Set PETSc options
        opts = PETSc.Options()
        opts.prefixPush(problem_prefix)
        if petsc_options is not None:
            for k, v in petsc_options.items():
                opts[k] = v
        opts.prefixPop()

        # Define KSP solver
        self._solver = PETSc.KSP().create(self.b.getComm().tompi4py())
        self._solver.setOperators(self.A)
        self._solver.setFromOptions()

        # Set matrix and vector PETSc options
        self.A.setOptionsPrefix(problem_prefix)
        self.A.setFromOptions()
        self.b.setOptionsPrefix(problem_prefix)
        self.b.setFromOptions()

    def solve(self, tol=1e-6, beta=1.0):
        i = 0

        while i < self.max_iterations:
            scatter_local_vectors(
                self.x,
                [si.x.petsc_vec.array_r for si in self.w],
                [
                    (si.function_space.dofmap.index_map, si.function_space.dofmap.index_map_bs)
                    for si in self.w
                ])
            self.x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

            # Assemble F(u_{i-1}) - J(u_D - u_{i-1}) and set du|_bc= u_D - u_{i-1}
            self.b.zeroEntries()
            fem.petsc.assemble_vector_block(self.b, self.F, self.J, bcs=self.bcs, x0=self.x)
            self.b.ghostUpdate(PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD)
            # Assemble Jacobian
            self.A.zeroEntries()
            fem.petsc.assemble_matrix_block(self.A, self.J, bcs=self.bcs)
            self.A.assemble()

            self._solver.solve(self.b, self.dx)
            assert self._solver.getConvergedReason() > 0, "Linear solver did not converge"
            offset_start = 0
            for s in self.w:
                num_sub_dofs = s.function_space.dofmap.index_map.size_local * s.function_space.dofmap.index_map_bs
                s.x.array[:num_sub_dofs] -= beta * self.dx.array_r[offset_start:offset_start + num_sub_dofs]
                s.x.scatter_forward()
                offset_start += num_sub_dofs
            # Compute norm of update

            correction_norm = self.dx.norm(0)
            print(f"Iteration {i}: Correction norm {correction_norm}")
            if correction_norm < tol:
                break
            i += 1

Mesh Generation 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
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, surfaces, loops and volume creation code] ...

gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 8*LMesh)
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")
2

I do not have pygmsh and hence cannot execute the code and comment with certainty. But I feel you might be missing a term like this in your weak form:

zero = fem.Constant(msh, default_scalar_type(0))
F3 = ufl.dot(zero*phi, phi)*dx
F = F0 + F1 + F2 + F3

Additionally - and without understanding your model in detail - use of lagrange multiplier generally requires homogeneous boundary condition on fs_l. If I am right, you will additionally need to add another term like F3 for fs_l as well to your weak form.

1 Like
3

Thank you for your answer. Unfortunately, this didn't help.

I don't think it is related to the Lagrange term because problem persists even when I use only "m" and "phi".

I've created another MWE that includes two domains: in one, linear elasticity is coupled with magnetostatics, and in the other, it's just linear elasticity. I created the function space for displacement on the full mesh, and a separate function space for the magnetic potential on the submesh. When I use the coupling across the whole domain without a mixed function space, everything works. But when I try to couple them across multiple domains, I get the same error.

MWE

import ufl
import numpy as np
from dolfinx import mesh, fem, default_scalar_type
from mpi4py import MPI
from dolfinx import io
from scifem import BlockedNewtonSolver
from customNewtonSolver import *
import numpy.typing as npt

nx, ny, nz = 20, 20, 20  # Number of elements in x and y directions
points = [[0., 0., 0.], [1., 1., 1.]]  # Domain corners
nElement = [nx, ny, nz]
msh = mesh.create_box(MPI.COMM_WORLD, points, nElement, mesh.CellType.tetrahedron)

def up(x, tol=1e-14):
    return x[2] >= 0.5 - tol

tdim = msh.topology.dim
cell_map = msh.topology.index_map(tdim)
num_cells_local = cell_map.size_local + cell_map.num_ghosts
marker = np.ones(num_cells_local, dtype=np.int32)

marker[mesh.locate_entities(msh, tdim, up)] = 2
cell_tags = mesh.meshtags(msh, msh.topology.dim, np.arange(num_cells_local, dtype=np.int32), marker)

sm, sm_cell_map, sm_vertex_map, _ = mesh.create_submesh(msh, cell_tags.dim, cell_tags.find(2))

sm.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

mesh_to_heat_entity = np.full(num_cells_local, -1, dtype=np.int32)
mesh_to_heat_entity[sm_cell_map] = np.arange(len(sm_cell_map), dtype=np.int32)
entity_maps = {sm: mesh_to_heat_entity}

# Material properties
E  = fem.Constant(msh, default_scalar_type(21000.0))
nu = fem.Constant(msh, default_scalar_type(3.0e-1))
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1 - nu) * (1 + nu))

fs_u = fem.functionspace(msh, ("Lagrange", 1, (msh.geometry.dim, )))
fs_phi = fem.functionspace(sm, ("Lagrange", 1))
V = ufl.MixedFunctionSpace(fs_u, fs_phi)

u = fem.Function(fs_u, name="Magnetization")
phi = fem.Function(fs_phi, name="Magnetic Potential")

du, dphi = ufl.TrialFunctions(V)
eps = ufl.sym(ufl.grad(u))
eps_var = ufl.variable(ufl.sym(ufl.grad(u)))
deps = ufl.sym(ufl.grad(du))
H = -ufl.grad(phi)
H_var = ufl.variable(-ufl.grad(phi))
dH = -ufl.grad(dphi)

a=ufl.as_vector([1, 0, 0])
psi_phi = - 0.5 * ufl.dot(H, H)
psi_u = (0.5 * lmbda * ufl.tr(eps) * ufl.tr(eps) + mu * ufl.inner(eps, eps))
psi_mix =  ufl.tr(eps)*ufl.dot(H, H) + ufl.tr(ufl.dot(eps, ufl.outer(a,a)))*ufl.dot(H, H) + ufl.tr(ufl.dot(eps, ufl.outer(H, a)))
psi = psi_phi + psi_u + psi_mix

sig = ufl.diff(psi, eps_var)
B = -ufl.diff(psi, H_var)

fdim = msh.topology.dim - 1
facets_left = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 0.0))
facets_right = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[0], 1.0))

dofs_left_u = fem.locate_dofs_topological(fs_u.sub(0), fdim, facets_left)
dofs_right_u = fem.locate_dofs_topological(fs_u.sub(0), fdim, facets_right)

bc_u_left = fem.dirichletbc(0.0, dofs_left_u, fs_u.sub(0))
bc_u_right = fem.dirichletbc(1.0, dofs_right_u, fs_u.sub(0))

bcs_u = [bc_u_left, bc_u_right]

dx = ufl.Measure("dx", domain=msh, subdomain_data=cell_tags, subdomain_id=2)

F =  ufl.dot(dH, B) * dx(2) + ufl.inner(deps, sig) * dx
residual = fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)

a, b = ufl.TestFunctions(V)
jac = ufl.derivative(F, u, a) + ufl.derivative(F, phi, b)
J = fem.form(ufl.extract_blocks(jac), entity_maps=entity_maps)

solver = customNewtonSolver(residual, J, [u, phi], bcs_u)
solver.solve()
4

I managed to run your script. I do not know if it is the main problem or just a symptom but something is certainly off with your matrices. Doing the following:

print("Size of A matrix: {}".format(self.A.getSize()))        
print("Size of x vector: {}".format(self.b.getSize()))

gives

Size of A matrix: (32634, 32634)
Size of b vector: 55566

With a quick look, I am unable to tell where exactly it went wrong. Perhaps it is somewhere when you extract blocks? I hope someone more knowledgeable with the internals is able to notice and respond. Could you possibly reconstruct the residue and jacobian as shown in this example? I just checked a nonlinear model involving mixed spaces I have which is made like that and it works.

5

Thank you very much for pointing this out and helping, following solved my problem, Trouble solving mixed-domain PDE with BlockedNewtonSolver - #5 by dokken

1 Like