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


109
views
0
4 weeks ago by

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.

Community: FEniCS Project

1 Answer


0
4 weeks ago by
Emek  
You overwrite your trial functions:

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
Thanks for your answer. I assumed that the problem is somewhere there. Unfortunately, when I delete the last three lines (so my initial guess is zero), I only get zero iterations and a zero solution. How do I correctly define those functions in the given context so that doesn't happen?
written 4 weeks ago by Rebecca Pflanze  
I cannot find out from the code if your problem is stated correctly or not. If you want to solve electromagnetism, why do you try to solve every component independently? You can solve a vector for the magnetic potential and a scalar for the electric potential out of Maxwell equations. You may find this paper interesting:

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
written 4 weeks ago by Emek  
I'm totally new to FEniCS and unfamiliar with the vector syntax. That's why I'm trying to get the code running like this before starting the research all over again, and I would really like to keep that for know. Mathematically, it shouldn't make a difference. Can you tell me how to define the functions Ax, At, Az properly for this problem, so that the solver won't start with an initial guess of zero? Any idea you might have would be a great help to me.
written 4 weeks ago by Rebecca Pflanze  
Hi, I'm trying to help, but the formulation is wrong in many ways.

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
written 4 weeks ago by Emek  
Ok, I'll have a look at those things. Thanks a lot for your time!
written 4 weeks ago by Rebecca Pflanze  

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.

written 4 weeks ago by Rebecca Pflanze  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »