Non-linear Multiphysics on subdomain error

dolfinx
,
1

Hello.

I am trying to solve a coupled system of non-linear equations that define on part of mesh, following examples Multiphysics: Solving PDEs on subdomains and Coupling PDEs of multiple dimensions from FEniCS-workshop. In following minimal example i get error:

 Cannot insert rows that do not exist in the IndexMap.

 File "/home/admin_ubuntu/Py/FenicsxSem/test.py", line 37, in <module>
  resid,

...<15 lines>...
    "pc_type": "lu",
    ^
 RuntimeError: Cannot insert rows that do not exist in the IndexMap.

But when i change function space domain from submesh to whole mesh error disappears
( QWspaceN = fem.functionspace(msh_S, (“CG”, 1) → fem.functionspace(msh, (“CG”, 1) )
QWspaceP = fem.functionspace(msh_S, (“CG”, 1)) → fem.functionspace(msh, (“CG”, 1))
But from Multiphysics: Solving PDEs on subdomains seem that i can define function space on sub mesh.

How to correct use multyphisics style on submesh with non-linear problem?
Please help. (dolfinx version is 0.10.0)
Min example:

import dolfinx
from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
import numpy as np
from dolfinx.fem.petsc import NonlinearProblem
import matplotlib.pyplot as plt
domain = [[-1, 0], [1, 1]]
n_cells = [100, 1000]
msh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, domain, n_cells, cell_type=mesh.CellType.quadrilateral )

mask = lambda x: (x[1] <= 0.6) & (x[1] >= 0.4)
cell_map = msh.topology.index_map(msh.topology.dim)
num_cells_local = cell_map.size_local + cell_map.num_ghosts
marker = np.ones(num_cells_local, dtype=np.int32)
marker[dolfinx.mesh.locate_entities(msh, msh.topology.dim, mask)] = 2

subdomains = dolfinx.mesh.meshtags(msh, msh.topology.dim, np.arange(num_cells_local, dtype=np.int32), marker)
msh_S, msh_S_cell_map, msh_S_vertex_map, _ = dolfinx.mesh.create_submesh(msh, msh.topology.dim, subdomains.values == 2 )

dx = ufl.dx( domain = msh, subdomain_data=subdomains )
dx2 = dx(2)
QWspaceN = fem.functionspace(msh_S, ("CG", 1)) 
QWspaceP = fem.functionspace(msh_S, ("CG", 1)) 
W = ufl.MixedFunctionSpace(QWspaceN, QWspaceP)
dqphiN2D, dqphiP2D = ufl.TrialFunctions(W) 
v_2d_n, v_2d_p = ufl.TestFunctions(W)
qphiN2D = fem.Function(QWspaceN)
qphiP2D = fem.Function(QWspaceP)
FN2d =   ufl.exp(qphiN2D)*ufl.dot( ufl.grad(qphiN2D) , ufl.grad(v_2d_n) ) * dx2 
FP2d = - ufl.exp(-qphiP2D)*ufl.dot( ufl.grad(qphiP2D) , ufl.grad(v_2d_p) ) * dx2
F =  FN2d + FP2d 
Jac = ufl.derivative(F, qphiN2D, dqphiN2D) + ufl.derivative(F, qphiP2D, dqphiP2D)

resid = ufl.extract_blocks(F)
JC = ufl.extract_blocks(Jac)
solverCurr = NonlinearProblem(
    resid,
    u = [qphiN2D, qphiP2D ],
    J = JC ,
    bcs  = [],
    entity_maps = [msh_S_cell_map],
    petsc_options={
            "snes_monitor": None,
            "ksp_type": "preonly",
            "pc_type": "lu",
            "pc_factor_mat_solver_type": "mumps",
            "mat_mumps_icntl_14": 120,
            "ksp_error_if_not_converged": True,
            "snes_error_if_not_converged": True,
                },
            petsc_options_prefix="signorini_",
        )
2

This line is the issue, as subdomains.values ==2 returns a boolean array of [True, False, …] rathre than the indicies of the cells.
To solve this, replace the offending line by:

msh_S, msh_S_cell_map, msh_S_vertex_map, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, subdomains.find(2)
)

or

msh_S, msh_S_cell_map, msh_S_vertex_map, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, subdomains.indices[subdomains.values == 2]
)
3

Thank you for correction. Now it is work.