Error: Unable to successfully call PETSc function 'MatSetValuesLocal'. Reason: PETSc error code is: 63 (Argument out of range). Where: This error was encountered inside ./dolfin/la/PETScMatrix.cpp

General
1

Hi, I am not able to solve the following error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_2676/3462195073.py in <module>
      5     dw = Constant(0.0005) #5e-4
      6     w += dw
----> 7     bcs = [DirichletBC(V.sub(1), Constant(0.), bottom), DirichletBC(V.sub(1), w, top)]
      8     u = Function(V)
      9     solve(a == l, u, bcs)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/dirichletbc.py in __init__(self, *args, **kwargs)
     74             if not hasattr(args[1], "_cpp_object"):
     75                 if isinstance(args[1], ufl.classes.Expr):
---> 76                     expr = project(args[1], args[0])  # FIXME: This should really be interpolaton (project is expensive)
     77                 else:
     78                     expr = Constant(args[1])

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
    130 
    131     # Assemble linear system
--> 132     A, b = assemble_system(a, L, bcs=bcs,
    133                            form_compiler_parameters=form_compiler_parameters)
    134 

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble_system(A_form, b_form, bcs, x0, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend)
    428         assembler.assemble(A_tensor, b_tensor, x0)
    429     else:
--> 430         assembler.assemble(A_tensor, b_tensor)
    431 
    432     return A_tensor, b_tensor

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 ./dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Here is the Minimum Working Example -

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

domain= Rectangle(Point(0.,0.), Point(1., 1.))
mesh = generate_mesh(domain, 30)

E = Constant(200e9) #E = 200 GPa 
nu = Constant(0.3)

model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

def eps(v):
    return sym(grad(v))

def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

f = 0.
f = Constant((0, f))

V = VectorFunctionSpace(mesh, 'CG', 1)
W = FunctionSpace(mesh, 'DG', 0)

du = TrialFunction(V)
d = du.geometric_dimension()
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx

def bottom(x, on_boundary):
    return near(x[1], 0.)

def top(x, on_boundary):
    return near(x[1], 1.)

w = Constant(0.001) #1e-3
bcs = [DirichletBC(V.sub(1), Constant(0.), bottom), DirichletBC(V.sub(1), w, top)]
u = Function(V)
solve(a == l, u, bcs)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, W)
print('max u:', u_magnitude.vector().max())
    
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
von_Mises = project(von_Mises, W)
von_Mises.vector().max()
print('max vonmises:', von_Mises.vector().max())

#code
sig0 = 450e6 # yield stress = 450 MPa

while von_Mises.vector().max() < sig0:
    dw = Constant(0.0005) #5e-4
    w += dw
    bcs = [DirichletBC(V.sub(1), Constant(0.), bottom), DirichletBC(V.sub(1), w, top)]
    u = Function(V)
    solve(a == l, u, bcs)
    u_magnitude = sqrt(dot(u, u))
    u_magnitude = project(u_magnitude, W)
    print('max u:', u_magnitude.vector().max())
    
    s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
    von_Mises = sqrt(3./2*inner(s, s))
    von_Mises = project(von_Mises, W)
    print('max vonmises:',von_Mises.vector().max())

Kindly guide me regarding the above.

Thank You.

2

Replace this with:

    w.assign(float(w) + float(dw))

Also note that there is room for major optimizations in this code, as you are using the projection function inside the loop, and the matrix A does not change over time (only the bc on the rhs does).

3

Thanks @dokken, for the help. Code is working fine now.

Based on your above suggestion, I tried to do a optimization, here is my code:

from dolfin import *
from mshr import *

domain= Rectangle(Point(0.,0.), Point(1., 1.))
mesh = generate_mesh(domain, 30)

# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
u, v = TrialFunction(W), TestFunction(W)
d = u.geometric_dimension()

# Introduce manually the material parameters
E = Constant(200e3) #E = 200 GPa 
nu = Constant(0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

# Constituive functions
def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

# Boundary conditions
top = CompiledSubDomain("near(x[1], 1.) && on_boundary")
bot = CompiledSubDomain("near(x[1], 0.) && on_boundary")

load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)

ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

# Variational form
unew, uold = Function(W), Function(W)

E_du = inner(sigma(u), epsilon(v))*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
solver_disp = LinearVariationalSolver(p_disp)

# Initialization of the iterative procedure and output requests
t = 0
u_r = 0.01
deltaT  = 0.005
iter = 0
conc_f = File ("./Displacement/displacement.pvd")
fname1 = open('ForcevsDisp.txt', 'w')
fname2 = open('stress[1, 1]vsstrain[1, 1].txt', 'w')

s = sigma(unew) - (1./3)*tr(sigma(unew))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
von_Mises = project(von_Mises, WW)
vm = von_Mises.vector().max()
print(vm)

 # yield stress = 600 MPa
while vm < 600:
    t += deltaT    
    load.t=t*u_r
    iter += 1
    
    solver_disp.solve()
    err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    print(err_u)
                    
    uold.assign(unew)

    s = sigma(unew) - (1./3)*tr(sigma(unew))*Identity(d)  # deviatoric stress
    von_Mises = sqrt(3./2*inner(s, s))
    von_Mises = project(von_Mises, WW)
    vm = von_Mises.vector().max()
    h = plot(von_Mises, title='Stress intensity');
    
    if vm > 600:
        break
        
    print(vm)
    print ('Iterations:', iter, ', Total time', t)

    if round(t*1e4) % 10 == 0:
        conc_f << unew

        Traction = dot(sigma(unew),n)
        fy = Traction[1]*ds(1)
                
        fname1.write(str(t*u_r) + "\t")
        fname1.write(str(assemble(fy)) + "\n")
                
        strain = epsilon(unew)
        strain_function = project(strain[1, 1], V)
        strain_at_point = strain_function((0.5, 1))
     
        stress = sigma(unew)
        stress_function = project(stress[1, 1], V)
        stress_at_point = stress_function((0.5, 1))
        
        fname2.write(str(strain_at_point) + "\t")
        fname2.write(str(stress_at_point) + "\n")
fname1.close()
fname2.close()
print ('Simulation completed')

What I expected here is err_u to decrease with each iteration and compare it with a tolerance value say tol = 1e-6 but in the first iteration itself, err_u = 2.895049351827927e-05 and it is been constant throughout all iterations.
How can I compare the l2 norm of the residual vector with a tolerance value and say that when that l2 norm value is less than tolerance I have given, the iterative solver converges ?
Thanks.

PS: Some part of the code is adapted from this paper and I basically wanted to use the method of this tutorial of using an manual Newton solver but I am not able to understand the fundamentals of working of an iterative solver and its use in solving pde’s. I would highly appreciate resources to understand the use of this kind of solver. I tried understanding this tutorial on Custom Newton solvers by Dokken but it is complicated for me because I am unable to understand the mathematics of it.I would highly appreciate any inputs on the above.

4

I don’t have bandwidth to go through all this code. You would have to reduce the number of things happening in the code for me to have any chance of debugging it. For instance, please remove all post-processing that is not required to illustrate the behavior.

Please also remove all conditionals from the code (like the “plane_stress” option, only have one option).

5

Here is the modified code -

from dolfin import *
from mshr import *

domain= Rectangle(Point(0.,0.), Point(1., 1.))
mesh = generate_mesh(domain, 30)

# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
u, v = TrialFunction(W), TestFunction(W)
d = u.geometric_dimension()

# Introduce manually the material parameters
E = Constant(200e3) #E = 200 GPa 
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Constituive functions
def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

# Boundary conditions
top = CompiledSubDomain("near(x[1], 1.) && on_boundary")
bot = CompiledSubDomain("near(x[1], 0.) && on_boundary")

load = Expression("t", t = 0.0, degree=1)
bc_u = [DirichletBC(W, Constant((0.0,0.0)), bot), DirichletBC(W.sub(1), load, top)]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)

ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

# Variational form
unew, uold = Function(W), Function(W)

E_du = inner(sigma(u), epsilon(v))*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
solver_disp = LinearVariationalSolver(p_disp)

# Initialization of the iterative procedure and output requests
t = 0
u_r = 0.01
deltaT  = 0.005
iter = 0

s = sigma(unew) - (1./3)*tr(sigma(unew))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
von_Mises = project(von_Mises, WW)
vm = von_Mises.vector().max()

 # yield stress = 600 MPa
while vm < 600:
    t += deltaT    
    load.t=t*u_r
    iter += 1
    
    solver_disp.solve()
    err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
    print(err_u)
                    
    uold.assign(unew)

    s = sigma(unew) - (1./3)*tr(sigma(unew))*Identity(d)  # deviatoric stress
    von_Mises = sqrt(3./2*inner(s, s))
    von_Mises = project(von_Mises, WW)
    vm = von_Mises.vector().max()
    
    if vm > 600:
        break

I cut short it to the best of my knowledge.

6

Why would it decrease? You keep on increasing the load in a linear fashion. You are solving a linear problem, and you are getting a linear response.

7

Okay, got it. I think I’m mixing 2 things here, so in this tutorial the author is solving a nonlinear problem, that is why the error of l2 norm keeps on decreasing ? Pardon me if I am getting confused in the trivial things.

Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
    loading.t = t
    A, Res = assemble_system(a_Newton, res, bc)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    Du.interpolate(Constant((0, 0)))
    print("Increment:", str(i+1))
    niter = 0
    while nRes/nRes0 > tol and niter < Nitermax:
        solve(A, du.vector(), Res, "mumps")
        Du.assign(Du+du)
        deps = eps(Du)
        sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
        local_project(sig_, W, sig)
        local_project(n_elas_, W, n_elas)
        local_project(beta_, W0, beta)
        A, Res = assemble_system(a_Newton, res, bc)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1
    u.assign(u+Du)
    p.assign(p+local_project(dp_, W0))
    sig_old.assign(sig)

Can you please tell me what is going on in the above code and when it is useful to apply above algorithm ? I took this code from the same tutorial.
Thanks.

8

Because h e has set up a newton problem, as illustrated by:

a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)

which is solving the Newton problem.

9

What if I want to have a displacement-dependent load and I want to stop it when particular displacement is reached ? which is the case with the real material testing machine, we give displacement and note load corresponding it.

10

This post is now diverging from its original topic.

I would suggest you sit down and formulate a minimal example of what you’d like to do on a built in mesh (say a unit square or a cube), and ask a new question with a clear description of what you want to achieve

11

Ok, thanks for this information. I would probably go through this tutorial and your Custom Newton Solvers tutorial to understand how these solvers work. Thanks for the tutorial and your time.

12

Sure. Sorry but last question, I was thinking of starting CFD via FEniCS, are there more sources other than FEniCSx tutorials ? My end goal is replicating the results of Numerical Study of Diodicity Mechanism in Different Tesla-Type Microvalves.
Thanks, I will close this post now.

13

I’ve written up a bit on splitting schemes at: Splitting schemes — Oasisx
and you have divergence conforming discretizations at dolfinx/python/demo/demo_navier_stokes.py at cb806540f6a0e52430a7ee54f2399adbfcaa6ab1 · FEniCS/dolfinx · GitHub
and there are probably many other sources out there I’ve missed.

1 Like
14

Thanks. I will start from these. :+1: