Unable to successfully call PETSc function 'MatSetValuesLocal'

bbataszai · April 27, 2022, 8:40pm

Hi everyone

I try to implement the bidomain equations for cardiac dynamics on the heart (H) with a surrounding region (T) on the mixed dimensional branch with the great MeshView feature… I had a similar post with a stupid error here.

Now the example code gives a Unable to successfully call PETSc function 'MatSetValuesLocal' error and I don’t understand why.

The system has two equations with two function spaces, one defined over the heart only and one defined over both regions, the variational form for the two equations are:

(1) \int_\text{H} \text{d}x M_i \nabla v^{n+1} \nabla \psi + \int_\text{H} \text{d}x (M_i + M_e) \nabla u^{n+1} \nabla \psi) = \int_\text{H} \text{d}x (\frac{v^{n+1} - v^n}{\Delta t} - I_\text{ion}) \psi
and
(2) \int_\text{H} \text{d}x M_i \nabla v^{n+1} \nabla \phi + \int_\text{H} \text{d}x (M_i + M_e) \nabla u^{n+1} \nabla \phi) + \int_\text{T} \text{d}x M_t \nabla u \nabla \phi = \int_{H+T} 0 \cdot \phi \text{d}x.

The system has to be solved for u^{n+1} and v^{n+1} with test functions \psi define over H and \phi defined over H+T:

Translated to the dolfin variational form, I get (a minimal working example is below)

# heart space + # heart and Torso function space
H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
HaT = df.FunctionSpace(mesh, "Lagrange", 1)
W = df.MixedFunctionSpace(H, HaT)
psi_H, psi_HaT = df.TestFunctions(W)
v_mt, phi_et = df.TrialFunctions(W)

# eq1
a1 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
    + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)

L1 = (
    -df.dot(v_mt, psi_H) * dx_heart
    + v_mtm1 * psi_H * dx_heart
    - I_ion * psi_H * dx_heart
)

# eq2
a2 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
    + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
    + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso
)

L2 = df.Constant(0) * psi_HaT * dx_all

a = a1 + a2
L = L1 + L2

Running thie minimal working example

import dolfin as df
mesh = df.UnitSquareMesh(20, 20)
interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
# mark domain
for cell in df.cells(mesh):
    x = cell.midpoint().array()
    subdomains[cell] = int(interior.inside(x, False))

# create interior mesh
mesh_heart = df.MeshView.create(subdomains, 1)
mesh_torso = df.MeshView.create(subdomains, 0)

# Function spaces
H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
HaT = df.FunctionSpace(mesh, "Lagrange", 1)
W = df.MixedFunctionSpace(H, HaT)

# Variational form
dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh_torso)
dx_all = df.Measure("dx", mesh)

psi_H, psi_HaT = df.TestFunctions(W)
v_mt, phi_et = df.TrialFunctions(W)

# last timestep v and ionic current
# can this be done like this?
v_mtm1 = df.Function(H)
I_ion = df.Function(H)

# conductivities
M_i = df.Constant(1)
M_e = df.Constant(1)
M_t = df.Constant(1)

# eq1
a1 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
    + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)

L1 = (
    -df.dot(v_mt, psi_H) * dx_heart
    + v_mtm1 * psi_H * dx_heart
    - I_ion * psi_H * dx_heart
)

# eq2
a2 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
    + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
    + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso
)


L2 = df.Constant(0) * psi_HaT * dx_all

a = a1 + a2
L = L1 + L2


# some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)


def boundary(x, on_boundary):
    return on_boundary


bc = df.DirichletBC(HaT, u_D, boundary)

sol = df.Function(W)
df.solve(a == L, sol, bc)

results in

PetscDolfinErrorHandler: line '511', function 'MatSetValues_SeqAIJ', file '[...]/SOURCE.CODES/petsc-3.15.0/src/mat/impls/aij/seq/aij.c',
                       : error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------
New nonzero at (210,210) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
------------------------------------------------------------------------------
PetscDolfinErrorHandler: line '1401', function 'MatSetValues', file '[...]/SOURCE.CODES/petsc-3.15.0/src/mat/interface/matrix.c',
                       : error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------

------------------------------------------------------------------------------
PetscDolfinErrorHandler: line '2261', function 'MatSetValuesLocal', file '[...]/SOURCE.CODES/petsc-3.15.0/src/mat/interface/matrix.c',
                       : error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------

------------------------------------------------------------------------------
PETSc error in '[...]/SOURCE.CODES/dolfin-2021-dev/dolfin/la/PETScMatrix.cpp', 'MatSetValuesLocal'
PETSc error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------
New nonzero at (210,210) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
------------------------------------------------------------------------------
Elapsed wall, usr, sys time: 0.000254798, 0, 0 ([MixedAssembler] Assemble cells)
Traceback (most recent call last):
  File "./temp.py", line 87, in <module>
    df.solve(a == L, sol, bc)
  File "[...]/.local/lib/python3.7/site-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "[...]/.local/lib/python3.7/site-packages/dolfin/fem/solving.py", line 264, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside [...]/SOURCE.CODES/dolfin-2021-dev/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  b55804ecca7d010ff976967af869571b56364975
*** -------------------------------------------------------------------------

Does anyone have an idea how I can debug this? I haven’t found a solution in the forum entries for mixed domain PDEs or for that specific error or in the undocumented examples here.

Any help would be very much appreciated!

tomfre · May 2, 2022, 8:15am

Hello,

this appears to be the same problem I found recently too, where the coupling of the subdomain function with the test functions on the whole domain gives the PETCs error 63. See my post here.

Did you find any solution for this yet?

bbataszai · May 2, 2022, 10:01am

Hi Tom,
unfortunately I don’t know how to continue here. I am waiting and I hope, that perhaps one of the experts finds this topic. Apparently coupling is difficult according to this thread. However in the power point here, slide 5, something very similar is shown (without full code unfortunately).

Probably you have already found this post but I wasn’t able to solve the coupled domain problem with this.

I have the feeling that I miss something when coupling a whole domain problem with a sub domain problem and I haven’t found an answer in the forum or somewhere else. I guess splitting the whole domain problem to each subdomain and imposing interior boundary conditions with Lagrange multiplier like here might be the solution.

I think it would be great to have an example on this from people much more experienced (or even developing) dolfin / dolfinx, perhaps someone finds the time However I really very much appreciate their work and I think its a great tool… Let’s see!

cdaversin · May 2, 2022, 11:59am

Hi,

Thanks for sending a minimal example, I was able to reproduce the error you mention.

There are two “issues” I noticed in your code :
1- You define (and use in your variational form) the Measure dx_torso = df.Measure("dx", mesh_torso) from mesh_torso, which is not used in the definition of your mixed function space. Your function space W = df.MixedFunctionSpace(H, HaT) only involve H built from mesh_heart and HaT from the whole mesh.

Since u is defined over H+T, I have replaced it with :

dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh)

(note that dx_all can then be replaced by dx_torso)

2- The term -df.dot(v_mt, psi_H) * dx_heart in L1 should be in the bilinear form a1 instead, since it involves a trial function :

a1 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
    + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)
a1 += df.dot(v_mt, psi_H) * dx_heart

L1 = - I_ion * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart

The code doesn’t crash anymore with these changes. Hope that helps,

Cécile

bbataszai · May 3, 2022, 9:10am

Hi Cécile,

Thank you a lot. Oh my god, I checked the code so often to not have a stupid error like this in it. Sorry for bothering you with this.

Something I don’t understand happens with the measures: Depending on how they are specified to integrate the last term of a2 only over the domain T the code either results in a Error: Unable to successfully call PETSc function 'MatSetValuesLocal' or it solves. Can you help me understand why this is the case? (perhaps related to this post from 2019)

Full, running examples are below. To result in the error I specify restrict the \int_T \text{d}x M_t \phi_{et} \psi_{HaT} term in a2 as follows:

# define dx_torso with subdomain data
dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh, subdomain_data=subdomains)
# Code as before [...]
a2 =  (
  M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart 
  + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
)
# Integrate only over $T$ -- next line: dx_torso -> dx_torso(0)
a2+= M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso(0)

To obtain running code, I use the measure that goes with the trial function:

dx_heart = df.Measure("dx", mesh_heart)
dx_torso = df.Measure("dx", mesh, subdomain_data=subdomains)
# Code as before [...]

# Measures in line 
a2 = (
   M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
   + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) *dx_torso(1)
   + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_torso(0)
)

It is kind of difficult to ask a clear question… Does this mean I have to be extra careful, which measure I use? Is the rule of thumb, to use the measure that goes with the trial function a good one?

The example codes are here:

MatSetValuesLocal - error example:

import dolfin as df

mesh = df.UnitSquareMesh(20, 20)
interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
# mark domain
for cell in df.cells(mesh):
    x = cell.midpoint().array()
    subdomains[cell] = int(interior.inside(x, False))

# create interior mesh
mesh_heart = df.MeshView.create(subdomains, 1)

# Function spaces
H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
HaT = df.FunctionSpace(mesh, "Lagrange", 1)
W = df.MixedFunctionSpace(H, HaT)

# Variational form
dx_heart = df.Measure("dx", mesh_heart)
dx_all = df.Measure("dx", mesh, subdomain_data=subdomains)

psi_H, psi_HaT = df.TestFunctions(W)
v_mt, phi_et = df.TrialFunctions(W)

# last timestep v and ionic current
# can this be done like this?
v_mtm1 = df.Function(H)
I_ion = df.Function(H)

# conductivities
M_i = df.Constant(1)
M_e = df.Constant(1)
M_t = df.Constant(1)

# eq1
a1 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
    + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)
a1 += df.dot(v_mt, psi_H) * dx_heart

L1 = -I_ion * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart

# eq2
a2 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
    + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_heart
    + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(0)
)


L2 = df.Constant(0) * psi_HaT * dx_all

a = a1 + a2
L = L1 + L2


# some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)


def boundary(x, on_boundary):
    return on_boundary


bc = df.DirichletBC(HaT, u_D, boundary)

sol = df.Function(W)
df.solve(a == L, sol, bc)

Working example

import dolfin as df

mesh = df.UnitSquareMesh(20, 20)
interior = df.CompiledSubDomain("std::max(fabs(x[0] - 0.5), fabs(x[1] - 0.5)) < 0.25")
subdomains = df.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
# mark domain
for cell in df.cells(mesh):
    x = cell.midpoint().array()
    subdomains[cell] = int(interior.inside(x, False))

# create interior mesh
mesh_heart = df.MeshView.create(subdomains, 1)

# Function spaces
H = df.FunctionSpace(mesh_heart, "Lagrange", 1)
HaT = df.FunctionSpace(mesh, "Lagrange", 1)
W = df.MixedFunctionSpace(H, HaT)

# Variational form
dx_heart = df.Measure("dx", mesh_heart)
dx_all = df.Measure("dx", mesh, subdomain_data=subdomains)

psi_H, psi_HaT = df.TestFunctions(W)
v_mt, phi_et = df.TrialFunctions(W)

# last timestep v and ionic current
# can this be done like this?
v_mtm1 = df.Function(H)
I_ion = df.Function(H)

# conductivities
M_i = df.Constant(1)
M_e = df.Constant(1)
M_t = df.Constant(1)

# eq1
a1 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_H)) * dx_heart
    + M_i * df.inner(df.grad(phi_et), df.grad(psi_H)) * dx_heart
)
a1 += df.dot(v_mt, psi_H) * dx_heart

L1 = -I_ion * psi_H * dx_heart + v_mtm1 * psi_H * dx_heart

# eq2
a2 = (
    M_i * df.inner(df.grad(v_mt), df.grad(psi_HaT)) * dx_heart
    + (M_i + M_e) * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(1)
    + M_t * df.inner(df.grad(phi_et), df.grad(psi_HaT)) * dx_all(0)
)


L2 = df.Constant(0) * psi_HaT * dx_all

a = a1 + a2
L = L1 + L2


# some Boundary condition to test integration (phi otherwise wit non-zero nullspace)
u_D = df.Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)


def boundary(x, on_boundary):
    return on_boundary


bc = df.DirichletBC(HaT, u_D, boundary)

sol = df.Function(W)
df.solve(a == L, sol, bc)

Again, Cécile, thank you for answering and for your efforts

Baltasar

Topic		Replies	Views
Error while trying to solve coupled PDEs variational formulation	4	713	March 22, 2021
Error: Unable to successfully call PETSc function 'MatSetValuesLocal'. Reason: PETSc error code is: 63 (Argument out of range) Errors	2	2104	July 6, 2021
PETSc function error Errors	5	463	April 29, 2022
PETSs Error Code 63 when trying to do Nonlinear Problem	2	677	February 15, 2022
PETSc error 63 for a simple nonlinear 1D simulation General	2	218	January 30, 2023

Unable to successfully call PETSc function 'MatSetValuesLocal'

Related Topics