PETSc function error

Errors
1

Hello! me and my group mates are trying to run a nonlinear structural problem to solve for the variable ho0 and are encountering this particular PETSc function 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 /build/dolfin-DneRIZ/dolfin-2019.2.0~git20201207.b495043/dolfin/la/PETScMatrix.cpp.
*** Process: 0

This is our minimal runing code(Example-ho0-possion ):

Example-ho0-possion:

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


tol=1E-12

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)



Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
plt.show()



class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1)


top = Top()

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

ds = ds(subdomain_data=boundarie) 

ho0=TrialFunction(V)
s=TestFunction(V)
t = TestFunction(V)

ho0= Function(V)

F0 = dot((-1*(ho0**2)/0.001),s)* ds(1) + dot((3-u),t)*ds(1)

ho0J=derivative(F0,ho0)



Ho0problem = NonlinearVariationalProblem(F0, ho0, None , ho0J)

Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0prm = Ho0solver.parameters

Ho0solver.solve()

ho0plot=plot(ho0, title='show ho0')
plt.colorbar(ho0plot)
plt.show()

print ('rho0 calculated')

Any help is much appreciated!

2

I believe the error is related to you having no integral over the domain in your definition of F0, i.e. you only integrate over the boundary. This was also discussed here.

Changing

F0 = dot((-1*(ho0**2)/0.001),s)* ds(1) + dot((3-u),t)*ds(1)

to

F0 = dot((-1*(ho0**2)/0.001),s)* ds(1) + dot((3-u),t)*ds(1) + Constant(0.0) * ho0 * s * dx

removes the error. However, the Newton solver doesn’t converge for the system you have specified in your code.

1 Like
4

Yes, thank you very much, with your tips we solved the problem of error63, but we did face the problem of Newton solver not converging.
We are currently taking these measures：

  1. Scaling up the dimension of the function space to 2.
  2. Increasing the tolerance value
  3. and Trying out SNES solver

But the result is not ideal.

Do you have any other suggestions？
Thank you so much!

6

Great!

Another generic thing you could try is the mesh resolution. However, it could also be that your problem isn’t well posed? Therefore, it might be worth looking into the actual equation you’re trying to solve. Unfortunately, I don’t have too much experience with this to be able to help you prove if it is or not. I’ve also not personally come across any equations which consist solely of surface integral terms, perhaps it could be related to that as well.

7

For problems on the surface of your mesh one should use MeshView developed by @cdaversin. There are plenty of posts regarding this subject on the forum

8

Hi,
Thank you for your help,
I’m in the same group as FFEIEI and we have been trying to implement the MeshView method suggested but we keep running into two problems depending on how we modify the code.

Here is a code for our current code

from dolfin import *
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt


tol=1E-12

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 2) # upped dimension

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
#plt.show()

###########################################################
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
    for k in facets(mesh):
        marker[k] = 0.5 - DOLFIN_EPS < k.midpoint().x() < 0.5 + DOLFIN_EPS

## Create submesh ##
submesh = MeshView.create(marker, 1)

Q = FunctionSpace(mesh, "Lagrange", 1)
LM = FunctionSpace(submesh, "Lagrange", 1)
W = MixedFunctionSpace(Q,LM)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1)

top = Top()

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

dBound = ds(subdomain_data=boundarie) 
dSurf = Measure("ds", domain=W.sub_space(1).mesh())

ho0,NotUsed0=TrialFunctions(W)
(s,s_lm)=TestFunctions(W)

#F0 = dot((-1*(ho0**2)/0.001),s)* dBound(1) + dot((3-u),s_lm)*dBound(1)
F0 = dot((-1*(ho0**2)/0.001),s)* dSurf + dot((3-u),s_lm)*dSurf

ho0J=derivative(F0,ho0)
print("ho0J = ", ho0J)
  
ho0= Function(W)

Ho0problem = NonlinearVariationalProblem(F0, ho0, None , ho0J)

Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0solver.solve()

ho0plot=plot(ho0, title='show ho0')
plt.colorbar(ho0plot)
plt.show()

print ('ho0 calculated')

The problems we keep running into are either

Computing derivative of form w.r.t. 'v_1^0'. Supply Function as a Coefficient
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/formmanipulations.py", line 85, in derivative
    raise RuntimeError("Computing derivative of form w.r.t. '{}'. Supply Function as a Coefficient".format(u))
  File "/home/ucarsten/Desktop/FEniCS/Example-ho0-poisson_mod.py", line 87, in <module>
    ho0J=derivative(F0,ho0)

Which seems to be a problem with the derivative code probably because we now got the two FunctionSpaces and the two TrialFunctions.

When we use the other F0 function

F0 = dot((-1*(ho0**2)/0.001),s)* dBound(1) + dot((3-u),s_lm)*dBound(1)
We get

'<' not supported between instances of 'Mesh' and 'Mesh'
  File "/home/ucarsten/Desktop/FEniCS/Example-ho0-poisson_mod.py", line 84, in <module>
    F0 = dot((-1*(ho0**2)/0.001),s)* dBound(1) + dot((3-u),s_lm)*dBound(1)

Sadly we have not been able to resolve either problem yet.
Again thank you for your help any further suggestions would be appreciated.