### Solving PDEs on mesh boundaries - unable to successfully call PETSc function 'MatSetValuesLocal'

197

views

2

In my other post https://www.allanswered.com/post/jrobe/solving-pdes-on-mesh-boundaries-better-way-to-calculate-gradient-along-mesh-boundary/ I was trying to solve a 1D PDE on the boundary of a square mesh. I have managed to get that example to work (either by setting "keepdiagonals=True" or by adding in a dummy term of zeros into my bilinear form) , but now I am having trouble with a slightly more complicated non-linear elasticity problem. Could anyone help me pinpoint the source of the error? Thanks

```
from dolfin import *
import numpy as np
import sympy as sp
ffc_options = {"optimize": True, "quadrature_degree":5}
if has_linear_algebra_backend("Epetra"):
parameters["linear_algebra_backend"] = "Epetra"
if __name__ == "__main__":
mesh = UnitSquareMesh(20,20)
# Define subdomains
# the bottom of the unit square, want to solve a pde along this boundary
class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] < DOLFIN_EPS
# The mesh boundary minus the bottom
class the_rest_bd(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[1]) >=0 and abs (x[0]) < DOLFIN_EPS or abs(x[0]-1) < DOLFIN_EPS) or x[1]>DOLFIN_EPS and on_boundary
# The inner mesh minus the bottom (NOT USED, assuming the and in the_rest_bd includes both boundaries and inner regions)
class the_rest_inner(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] -1) > DOLFIN_EPS and x[0] > DOLFIN_EPS and \
abs(x[1] -1) > DOLFIN_EPS and x[1] > DOLFIN_EPS
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domain = MeshFunction("size_t", mesh, mesh.topology().dim() )
boundaries.set_all(0)
sub_domain.set_all(0)
the_rest_label = 1
the_rest_bd = the_rest_bd()
the_rest_inner = the_rest_inner()
the_rest_bd.mark(boundaries, the_rest_label)
the_rest_inner.mark(sub_domain, the_rest_label)
bottom_label = 2
bottom_bd = bottom()
bottom_bd.mark(boundaries,bottom_label)
ds = ds(subdomain_data=boundaries)
gamma = Expression( ('x[0]','0'),degree = 1)
our_norm = FacetNormal(mesh)
our_tang = as_vector([our_norm[1], -our_norm[0]])
Gsub1 = our_tang
# Define function spaces
Ve = VectorElement("CG",mesh.ufl_cell(), degree=1 ,dim=2)
V = FunctionSpace(mesh, Ve)
delU = Function(V)
VScalar = FunctionSpace(mesh,'CG',1)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
pressureFunc = Function(VScalar)
gsub1 = Gsub1 + dot(grad(u),our_tang)
squG = sqrt(dot(Gsub1,Gsub1))
squg = sqrt(dot(gsub1,gsub1))
lam_nonLin = squg/squG
Ic_incompr = lam_nonLin + 1./lam_nonLin
mu = Constant(1.)
psi = mu/2.*( Ic_incompr - 3 )
# Pin, Pin
bc = DirichletBC(V, Constant(('0','0')), '(fabs(x[0]) < DOLFIN_EPS && fabs(x[1]) < DOLFIN_EPS) || (fabs(x[0]-1.) < DOLFIN_EPS && fabs(x[1]) < DOLFIN_EPS)', method = 'pointwise')
# the rest of the domain will be set to 0
bc_the_rest_bd = DirichletBC(V, Constant([0,0]), the_rest_bd)
bc_the_rest_inner = DirichletBC(V, Constant([0,0]), the_rest_inner)
bcs = [bc_the_rest_bd, bc_the_rest_inner, bc]
thickness = Constant(0.01)
f = Expression("0", degree=1)
dummy_vec = Expression(("1","1"), degree=1)
Pi_int = thickness*psi*sqrt(dot(Gsub1,Gsub1))*ds(2)
pressureFunc.vector()[:] = 0.01
D_Pi_int = derivative(Pi_int,u,v) #
F = D_Pi_int - pressureFunc*dot( v, our_norm )*ds(2) #residual DPI[v] ... \int t_0*dv*ds
K = derivative(F,u,du)
F = F + f*dot(dot(grad(u), dummy_vec),dot(grad(v),dummy_vec))*dx # add dummy zero terms so that matrix is properly preallocated (because no "keep diagonals" option for newton solver)
solve(F == 0, u, bcs, J = K, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}} ,form_compiler_parameters=ffc_options)
```

Currently I am getting the following error:

```
Traceback (most recent call last):
File "<ipython-input-2-02884d3593b2>", line 1, in <module>
runfile('/Users/alexanderniewiarowski/fenics/ale_disk/source/ale-stuff/skinny_inflation_2d.py', wdir='/Users/alexanderniewiarowski/fenics/ale_disk/source/ale-stuff')
File "/anaconda3/envs/py36/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
execfile(filename, namespace)
File "/anaconda3/envs/py36/lib/python3.6/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)
File "/Users/alexanderniewiarowski/fenics/ale_disk/source/ale-stuff/skinny_inflation_2d.py", line 97, in <module>
solve(F == 0, u, bcs, J = K, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}} ,form_compiler_parameters=ffc_options)
File "/anaconda3/envs/py36/lib/python3.6/site-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/anaconda3/envs/py36/lib/python3.6/site-packages/dolfin/fem/solving.py", line 349, 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 /Users/travis/miniconda3/conda-bld/fenics_1519648192279/work/dolfin-2017.2.0.post0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:
*** -------------------------------------------------------------------------
```

Community: FEniCS Project

2

Since the derivative of F is taken before adding the dummy term, the dummy term has no effect on the bilinear form used in each Newton iteration. If you compute K after adding the dummy term, it gets rid of this error, but then the nonlinear solver diverges, so there's more work to do...

written
4 months ago by
David Kamensky

Please login to add an answer/comment or follow this question.