Poisson equation with pure Neumann BC and NonLinear solver, PETSc error code is: 63 (Argument out of range)

General
1

Hi everyone,

I am trying to solve the Poisson problem with pure Neumann boundary conditions in FEniCS 2019.1.

I looked at the example here:

https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/neumann-poisson/demo_neumann-poisson.py.html

My code is mostly copied from there, but I made a few changes:

from dolfin import *

# Create mesh
#mesh = UnitSquareMesh.create(64, 64, CellType.Type.quadrilateral)
mesh = UnitSquareMesh(64,64)

# Build function space with Lagrange multiplier
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)

# Define variational problem
#s1 = Function(W)
#u, c = split(s1)
(u, c) = Function(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
u_mean = Constant(0.32)
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + u_mean*d*dx
F = a - L

# Compute solution
w = Function(W)
solve(F == 0, w)#, solver_parameters={"newton_solver":{"relative_tolerance":1e-6}})
(u, c) = w.split()

# Save solution in VTK format
file = File("neumann_poisson_mod_NL.pvd")
u.rename("f","f")
file << u

Unfortunately the program fails during the first newton iteration with the following error:

Newton iteration 0: r (abs) = 3.212e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "/home/pbi/CASES/tutorials/onlyNeumannMod_NL/onlyNeumannMod_NL.py", line 26, in <module>
    solve(F == 0, w)#, solver_parameters={"newton_solver":{"relative_tolerance":1e-6}})
    ^^^^^^^^^^^^^^^^
  File "/home/pbi/miniconda3/envs/t02/lib/python3.11/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/pbi/miniconda3/envs/t02/lib/python3.11/site-packages/dolfin/fem/solving.py", line 266, 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 /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1668433332839/work/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  3ea2183fbfe7277de9f16cbe1a9ffaab133ba1fa
*** -------------------------------------------------------------------------

I also looked this “PETSc error code is: 63 (Argument out of range).” up in this forum and tried a few solutions (like using the split(…) function), but nothing helped so far. Does someone know a solution to this problem?

2

This should be

w = Function(W)
(u, c) = split(w)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
u_mean = Constant(0.32)
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + u_mean*d*dx
F = a - L

# Compute solution
solve(F == 0, w)
1 Like
3

Hi, I am trying to calculate the three-bit Poisson problem using the UnitHyperCube(divisions) function with the following code:

from __future__ import print_function
from fenics import *
import numpy as np

def solver(f, u_D, divisions, degree=1):
    # Create mesh and define function space
    # mesh = UnitCubeMesh(Nx, Ny, Nz)
    mesh = UnitHyperCube(divisions)
    V = FunctionSpace(mesh, 'P', degree)

    # Define boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

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

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

    return u

def UnitHyperCube(divisions):
    mesh_classes = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
    d = len(divisions)
    mesh = mesh_classes[d - 1](*divisions)
    # mesh = UnitSquareMesh(divisions[0], divisions[1])
    return mesh

def run_solver():
    "Run solver to compute and post-process solution"

    # Set up problem parameters and call solver
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1] + 2*x[2]*x[2]', degree=2)
    f = Constant(-6.0)
    u = solver(f, u_D, '3', 1)

    # Save solution to file in VTK format
    vtkfile = File('poisson_solver222/solution.pvd')
    vtkfile << u

if __name__ == '__main__':
    run_solver()

But when I run it, I get an error as below:

TypeError                                 Traceback (most recent call last)
Cell In[2], line 49
     46     vtkfile << u
     48 if __name__ == '__main__':
---> 49     run_solver()

Cell In[2], line 42, in run_solver()
     40 u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1] + 2*x[2]*x[2]', degree=2)
     41 f = Constant(-6.0)
---> 42 u = solver(f, u_D, '3', 1)
     44 # Save solution to file in VTK format
     45 vtkfile = File('poisson_solver222/solution.pvd')

Cell In[2], line 8, in solver(f, u_D, divisions, degree)
      5 def solver(f, u_D, divisions, degree=1):
      6     # Create mesh and define function space
      7     # mesh = UnitCubeMesh(Nx, Ny, Nz)
----> 8     mesh = UnitHyperCube(divisions)
      9     V = FunctionSpace(mesh, 'P', degree)
     11     # Define boundary condition

Cell In[2], line 32, in UnitHyperCube(divisions)
     30 mesh_classes = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
     31 d = len(divisions)
---> 32 mesh = mesh_classes[d - 1](*divisions)
     33 # mesh = UnitSquareMesh(divisions[0], divisions[1])
     34 return mesh

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.generation.UnitIntervalMesh(arg0: int)
    2. dolfin.cpp.generation.UnitIntervalMesh(arg0: MPICommWrapper, arg1: int)

Invoked with: '3'

May I ask if the type of divisions is string? And how can I solve the above error?

4

It should be an integer, ie 3 should work. Please note that your post is off-topic to the original post.

5

I’m so sorry. But when I change ‘3’ to 3, it shows the following error：

TypeError                                 Traceback (most recent call last)
Cell In[3], line 49
     46     vtkfile << u
     48 if __name__ == '__main__':
---> 49     run_solver()

Cell In[3], line 42, in run_solver()
     40 u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1] + 2*x[2]*x[2]', degree=2)
     41 f = Constant(-6.0)
---> 42 u = solver(f, u_D, 3, 1)
     44 # Save solution to file in VTK format
     45 vtkfile = File('poisson_solver222/solution.pvd')

Cell In[3], line 8, in solver(f, u_D, divisions, degree)
      5 def solver(f, u_D, divisions, degree=1):
      6     # Create mesh and define function space
      7     # mesh = UnitCubeMesh(Nx, Ny, Nz)
----> 8     mesh = UnitHyperCube(divisions)
      9     V = FunctionSpace(mesh, 'P', degree)
     11     # Define boundary condition

Cell In[3], line 31, in UnitHyperCube(divisions)
     29 def UnitHyperCube(divisions):
     30     mesh_classes = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
---> 31     d = len(divisions)
     32     mesh = mesh_classes[d - 1](*divisions)
     33     # mesh = UnitSquareMesh(divisions[0], divisions[1])

TypeError: object of type 'int' has no len()
6

Wait, you wrote the unit hypercube function.

It should be quite clear from your own code that divisions should be a list of integers with number of elements in each direction.