### PETSc error: Argument out of range (error code 63)

Hi everyone,

I'm facing an error code I don't really understand. I don't seem to be the first person to encounter this, but I couldn't apply any of the solutions I found to my code. This is it:

```
from __future__ import print_function
from fenics import *
from dolfin import *
from mshr import *
from math import sin, cos, pi
import matplotlib.pyplot as plt
import numpy as np
# some constants
mu_ri = 1.0e6
mu_rc = 1
n = 50
# coil currents
I_left = 10
I_right = 10
# following current densities
J_left = I_left*n/(109.4e-3*27.88e-3)
J_right = I_right*n/(109.4e-3*27.88e-3)
mesh = Mesh("Feldgenerator.xml")
subdomains = MeshFunction("size_t", mesh, "Feldgenerator_physical_region.xml")
V = FunctionSpace(mesh, "CG", 1)
bc = DirichletBC(V, Constant(0.0), 'on_boundary')
dx = Measure("dx")(subdomain_data=subdomains)
# define relative permeability for subdomains
class Permeability(Expression):
def __init__(self, subdomains, **kwargs):
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 0 or 18 or 19:
values[0] = 1.0 # 0 background, 18/19 isolations
elif self.subdomains[cell.index] == 1:
values[0] = mu_ri # iron yoke
else:
values[0] = mu_rc
# permeability
mu = 4*pi*1e-7*Permeability(subdomains, degree=1)
# define variational problem for x, y, z
Ax = TrialFunction(V)
Ay = TrialFunction(V)
Az = TrialFunction(V)
Ax = project(1.0, V) # non-zero initial guess
Ay = project(1.0, V)
Az = project(1.0, V)
vx = TestFunction(V)
vy = TestFunction(V)
vz = TestFunction(V)
# left hand sides of variational problem
ax = dot(nabla_grad(Ax), nabla_grad(vx))*dx
ay = dot(nabla_grad(Ay), nabla_grad(vy))*dx
az = dot(nabla_grad(Ax), nabla_grad(vx))*dx
# right hand sides of variational problem
# left coil with subdomains
Lxl1 = mu*J_left*vx*dx(2) # top
Lxl2 = 1/sqrt(2)*mu*J_left*vx*dx(3) # top left
Lxl3 = -1/sqrt(2)*mu*J_left*vx*dx(5) # bottom left
Lxl4 = -mu*J_left*vx*dx(6) # bottom
Lxl5 = -1/sqrt(2)*mu*J_left*vx*dx(7) # bottom right
Lxl6 = 1/sqrt(2)*mu*J_left*vx*dx(9) # top right
Lzl1 = -1/sqrt(2)*mu*J_left*vz*dx(3) # top left
Lzl2 = -mu*J_left*vz*dx(4) # left
Lzl3 = -1/sqrt(2)*mu*J_left*vz*dx(5) # bottom left
Lzl4 = 1/sqrt(2)*mu*J_left*vz*dx(7) # bottom right
Lzl5 = mu*J_left*vz*dx(8) # right
Lzl6 = 1/sqrt(2)*mu*J_left*vz*dx(9) # top right
Lxl = Lxl1 + Lxl2 + Lxl3 + Lxl4 + Lxl5 + Lxl6
Lzl = Lzl1 + Lzl2 + Lzl3 + Lzl4 + Lzl5 + Lzl6
#right coil with subdomains
Lxr1 = -mu*J_right*vx*dx(10) # top
Lxr2 = -1/sqrt(2)*mu*J_right*vx*dx(11) # top left
Lxr3 = 1/sqrt(2)*mu*J_right*vx*dx(13) # bottom left
Lxr4 = mu*J_right*vx*dx(14) # bottom
Lxr5 = 1/sqrt(2)*mu*J_right*vx*dx(15) # bottom right
Lxr6 = -1/sqrt(2)*mu*J_right*vx*dx(17) # top right
Lzr1 = -1/sqrt(2)*mu*J_right*vz*dx(11) # top left
Lzr2 = -mu*J_right*vz*dx(12) # left
Lzr3 = -1/sqrt(2)*mu*J_right*vz*dx(13) # bottom left
Lzr4 = 1/sqrt(2)*mu*J_right*vz*dx(15) # bottom right
Lzr5 = mu*J_right*vz*dx(16) # right
Lzr6 = 1/sqrt(2)*mu*J_right*vz*dx(17) # top right
Lxr = Lxr1 + Lxr2 + Lxr3 + Lxr4 + Lxr5 + Lxr6
Lzr = Lzr1 + Lzr2 + Lzr3 + Lzr4 + Lzr5 + Lzr6
Lx = Lxl + Lxr
Ly = 0
Lz = Lzl + Lzr
Fx = ax - Lx
Fy = ay - Ly
Fz = az - Lz
# solve variational problem
solve(Fx == 0, Ax, bc)
solve(Fy == 0, Ay, bc)
solve(Fz == 0, Az, bc)
# compute magnetic field (B = curl(A))
W = VectorFunctionSpace(mesh, "CG", 1)
Bx = Az.dx(1) - Ay.dx(2)
By = Ax.dx(2) - Az.dx(0)
Bz = Ay.dx(0) - Ax.dx(1)
B = project(as_vector((Bx, By, Bz)), W)
# plot & show
plot(B, mode = "glyphs")
plt.show()
```

It would have been prettier with vectors, but I'm still a complete novice to FEniCS I had to keep the syntax simple. Anyway, when I try to run it, I get the following message:

```
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.725e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.893e-13 (tol = 1.000e-10) r (rel) = 1.677e-16 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.725e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.893e-13 (tol = 1.000e-10) r (rel) = 1.677e-16 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.698e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File "/home/rebecca/Documents/Bachelorarbeit/3D_model.py", line 117, in <module>
solve(Fz == 0, Az, bc)
File "/home/rebecca/anaconda3/lib/python2.7/site-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/home/rebecca/anaconda3/lib/python2.7/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 /feedstock_root/build_artefacts/fenics_1519993789708/work/dolfin-2017.2.0.post0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: 774aa9b05f5a21fcf3d1bd632e722933a05fdb45
*** -------------------------------------------------------------------------
[Finished in 11.198s]
```

Looks like it's calculating the x and y components just fine, but has problems with z? I'd be happy for any idea how to fix this. Thanks in advance.

### 1 Answer

Ax = TrialFunction(V)

Ay = TrialFunction(V)

Az = TrialFunction(V)

Ax = project(1.0, V) # non-zero initial guess

Ay = project(1.0, V)

Az = project(1.0, V)

such that your system does not have any bilinear form anymore. But the error is very confusing and not very helpful.

Best, Emek

Abali, B. E., & Reich, F. A. (2017). Thermodynamically consistent derivation and computation of electro-thermo-mechanical systems for solid bodies. Computer Methods in Applied Mechanics and Engineering, 319, 567-595.

Have a look at my website http://bilenemek.abali.org

and the code in FEniCS is here: http://www.lkm.tu-berlin.de/ComputationalReality/

If you have specific questions, do not hesitate to ask here or via e-mail directly.

All the best, Emek

First, I do not think so that solving each component in an order that you have chosen is equivalent to solve them together. Consider B=curl A and you want to solve div B = 0 (one of the Maxwell equations). The components of A are clearly coupled, you cannot solve each component separately.

Second, you are mixing solving method of linear weak form with the solving method of nonlinear weak form. That is why I have put a link for you to see FEniCS code in Python (as you do) for digging in. You can also have a look at the official tutorial under https://fenicsproject.org/tutorial/ and see a linear pde's solution method and then a nonlinear pde's solution method.

Best, Emek

Hi Emek,

I solved this issue and tried to correct the things you mentioned. Hopefully, my current problem is a bit more specific, so perhaps you'd like to have a look at it: https://www.allanswered.com/post/kgbwm/vector-potential-and-magnetic-field-calculated-as-zero/

Would be great if you can help.