# PETSc error code is: 63 (Argument out of range) while assembling system

**URL:** <https://fenicsproject.discourse.group/t/petsc-error-code-is-63-argument-out-of-range-while-assembling-system/6598>\
**Category:** Errors\
**Created:** [September 16, 2021, 2:08pm UTC](https://fenicsproject.discourse.group/t/petsc-error-code-is-63-argument-out-of-range-while-assembling-system/6598 "2021-09-16T14:08:53Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Emanuel\_Gebauer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/emanuel_gebauer/32/3004_2.png) [@Emanuel\_Gebauer](https://fenicsproject.discourse.group/u/Emanuel_Gebauer)\
**Post date:** [September 16, 2021, 2:08pm UTC](https://fenicsproject.discourse.group/t/petsc-error-code-is-63-argument-out-of-range-while-assembling-system/6598/1 "2021-09-16T14:08:53Z")

</div>

Hey,

I stumbled into the following error while trying to assemble a mixed variational system:  
\*\*\* -------------------------------------------------------------------------  
\*\*\* Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.  
\*\*\* Reason: PETSc error code is: 63 (Argument out of range).  
\*\*\* Where: This error was encountered inside /build/dolfin-GfTndI/dolfin-2019.2.0~git20201207.b495043/dolfin/la/PETScMatrix.cpp.  
\*\*\* Process: 0  
\*\*\* DOLFIN version: 2019.2.0.dev0  
\*\*\* Git changeset: unknown  
\*\*\* -------------------------------------------------------------------------

The system I’m trying to simulate follows the following equations:  
\begin{aligned} Eq0: &&-\nabla^2 u0 &= &1\\ Eq1: &&-\nabla^2 u0 -\nabla^2 u1 &= &1\\ \end{aligned}  
with a Dirichlet Boundary at the left with u0\_{left}=1, u1\_{left}=1. On the other boundaries there should be no flow going outside: \frac{\delta u\_0}{\delta \overrightarrow{n}}=0 and \frac{\delta u\_1}{\delta \overrightarrow{n}}=0.

Hereby u1 is only defined on the yellow subdomain, whereas u0 is defined on the whole domain (purple and yellow):

 ![SubDomain](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/33e874e86b0ddeb94c4e230156d7e7081f90a109.png)

Defining the TestFunctions v0 and v1 for equation 0 and 1 results in:  
\begin{aligned} Eq0: && \int\_{\Omega\_0} \nabla u0 \nabla v0 dx&=& &\int\_{\Omega\_0} v0 dx + \underbrace{\int\_{\delta\Omega\_0} \frac{\delta u0}{\delta \overrightarrow{n}} v0 ds}\_{=0} \\ Eq1: &&\int\_{\Omega\_1}\nabla u0\nabla v1 dx +\int\_{\Omega\_1}\nabla u1\nabla v1 dx &= &&\int\_{\Omega\_1} v1 dx + \underbrace{\int\_{\delta\Omega\_1} \frac{\delta u1}{\delta \overrightarrow{n}} v1 ds}\_{=0} + \underbrace{\int\_{\delta\Omega\_0} \frac{\delta u0}{\delta \overrightarrow{n}} v1 ds}\_{\neq0}\\ \end{aligned}  
Therefore my code is the following:

```auto
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(64, 64)

# Create Subdomain
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return x[1] <= 0.75 + tol and x[1] >= 0.25 - tol

materials = MeshFunction("size_t", mesh, 2)
materials.set_all(0)
subdomain = Omega_0()
subdomain.mark(materials, 1)

mesh_subdomain = MeshView.create(materials, 1)

# Introduce FunctionSpaces for different variables
v0 = FunctionSpace(mesh, 'P', 1)
v1 = FunctionSpace(mesh_subdomain, 'P', 1)

V = MixedFunctionSpace(v0, v1)

## Init the functions
u0, u1 = TrialFunctions(V) 
v0, v1 = TestFunctions(V)

# dx and ds for the integrals
dx_0 = Measure("dx", domain = V.sub_space(0).mesh())
dx_1 = Measure("dx", domain = V.sub_space(1).mesh())

ds_1 = Measure("ds", domain = V.sub_space(1).mesh())

# Boundary conditions
def boundary_Left(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0., tol)
bc_0 = DirichletBC(V.sub_space(0), Constant(1.), boundary_Left) 
bc_1 = DirichletBC(V.sub_space(1), Constant(1.), boundary_Left) 
bcs = [bc_0, bc_1]

n_1 = FacetNormal(V.sub_space(1).mesh())
### Set up the variational problem
a_0 = dot(grad(u0), grad(v0)) * dx_0     
a_1 = dot(grad(u1), grad(v1)) * dx_1 + \
       dot(grad(u0), grad(v1)) * dx_1 - \
       dot(n_1, grad(u0)) * v1 * ds_1    
a = a_0 + a_1

L_0 = Constant(1.) * v0 * dx_0
L_1 = Constant(1.) * v1 * dx_1
L = L_0 + L_1

du_ = Function(V)
system = assemble_mixed_system(a==L, du_, bcs)

```

At the last line (assemble\_mixed\_system) I get said error.  
Confusingly if I use a differnt Subdomain Omega\_1, then assembling the system works.

```auto
# Create Subdomain
class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return x[1] <= 0.75 + tol

```

Can anyone point me to what I am missing?  
Thanks in advance!

Emanuel

---

<div class="post-metadata">

**Author:** ![cdaversin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/cdaversin/32/44_2.png) [@cdaversin](https://fenicsproject.discourse.group/u/cdaversin)\
**Post date:** [September 28, 2021, 8:31am UTC](https://fenicsproject.discourse.group/t/petsc-error-code-is-63-argument-out-of-range-while-assembling-system/6598/2 "2021-09-28T08:31:01Z")

</div>

Hi Emanuel,  
It seems that the issue here is the surface integral `dot(n_1, grad(u0)) * v1 * ds_1` mixing functions `u0` and `v1` (not sure why it doesn’t crash with your other subdomain, but the result might not be right). Mixed-domains surface integrals are currently not supported. From the [paper introducing mixed-dimensional features](https://arxiv.org/pdf/1911.01166.pdf) : “the redefined measures should define integration over cells for a lower dimensional mesh rather than integration over facets for a higher dimensional mesh.”

You could possibly define an other lower-dimensional submesh corresponding to your boundary, with `dx`-type integral instead.

---

<div class="post-metadata">

**Author:** ![Emanuel\_Gebauer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/emanuel_gebauer/32/3004_2.png) [@Emanuel\_Gebauer](https://fenicsproject.discourse.group/u/Emanuel_Gebauer)\
**Post date:** [September 28, 2021, 11:26am UTC](https://fenicsproject.discourse.group/t/petsc-error-code-is-63-argument-out-of-range-while-assembling-system/6598/3 "2021-09-28T11:26:27Z")

</div>

Dear cdaversin,

thank you for the reply. In the meantime however I reimplemented my system using MixedElements instead of a MixedFunctionSpace. I think now it works as intended.

**Edit:**  
For reference: a quick and dirty implementation of the problem using MixedElements:

```auto
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(64, 64)

# Create Subdomain
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return x[1] <= 0.75 + tol and x[1] >= 0.25 - tol

materials = MeshFunction("size_t", mesh, 2)
materials.set_all(0)
subdomain = Omega_0()
subdomain.mark(materials, 1)

facet_materials = MeshFunction("size_t", mesh, 1)
facet_materials.set_all(0)
subdomain = Omega_0()
subdomain.mark(facet_materials, 1)

# Introduce FunctionSpaces for different variables
V0 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
V1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)

V = FunctionSpace(mesh, MixedElement([V0, V1]))

## Init the functions
u0, u1 = TrialFunctions(V) 
v0, v1 = TestFunctions(V)

# dx and ds for the integrals
dx = Measure("dx", domain = V.sub(0).mesh(), subdomain_data=materials)

ds_1 = Measure("ds", domain = V.sub(1).mesh(), subdomain_data=facet_materials, subdomain_id=1)

# Boundary conditions
def boundary_Left(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0., tol)

def boundary_Left_Omega0(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0., tol) and x[1] <= 0.75 + tol and x[1] >= 0.25 - tol
bc_0 = DirichletBC(V.sub(0), Constant(1.), boundary_Left) 
bc_1 = DirichletBC(V.sub(1), Constant(0.), boundary_Left_Omega0) 
bcs = [bc_0, bc_1]

n_1 = FacetNormal(V.sub(1).mesh())
### Set up the variational problem
a_0 = dot(grad(u0), grad(v0)) * dx     
a_1 = dot(grad(u1), grad(v1)) * dx(1) + \
       dot(grad(u0), grad(v1)) * dx(1) - \
       dot(n_1, grad(u0)) * v1 * ds_1    
a = a_0 + a_1

L_0 = Constant(1.) * v0 * dx
L_1 = Constant(1.) * v1 * dx(1)
L = L_0 + L_1

du_ = Function(V)
A = assemble(a, keep_diagonal=True)
b = assemble(L)
for bc in bcs:
    bc.apply(A,b)

A.ident_zeros()

solver = KrylovSolver('gmres', 'ilu')
prm = solver.parameters
prm['absolute_tolerance'] = 1E-14
prm['relative_tolerance'] = 1E-7
prm['maximum_iterations'] = 1000
prm['monitor_convergence'] = True
solver.solve(A, du_.vector(), b)
u0, u1 = du_.split()
plot(u1)
plt.show()

```