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

197
views
2
4 months ago by

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

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)

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

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
***
***
*** 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