Unable to successfully call PETSc function 'MatSetValuesLocal'

variational formulation
5

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 :slight_smile:

Baltasar