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!