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

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.

Similar posts:
Search »