Error in Newton solver

Linear Algebra
1

I have written the following code to calculate Magnetic potential in 3D.The following is my code.

from fenics import *
from mshr import *
from math import *
from dolfin import *      

# Create Geometry

Rair = 5
cyltot = Cylinder(Point(0,0,10),Point(0,0,0),Rair,Rair)
mesh_tot = generate_mesh(cyltot,64)
vtkfile = File('/home/fenics/shared/Cylinder/Mesh_tot.pvd')
vtkfile << mesh_tot
V = FunctionSpace(mesh_tot,'P',1)`


# Subdomains

def boundary(x,on_boundary):
    return on_boundary

import numpy as np
tol = 1E-14

class Copper(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) <= 1 + tol

    
class Air(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) >= 1 - tol


materials = MeshFunction('size_t',mesh_tot,mesh_tot.topology().dim())
Subdom_cu = Copper()
Subdom_air = Air()
materials.set_all(0)
Subdom_cu.mark(materials, 1)

file = File("/home/fenics/shared/Cylinder/Materials.pvd")
file << materials

# Defining the parameters
RCu = 1/1.256637E-6
RAir = 1/1.256629E-6
JCu = 6/(np.pi)
JAir = 0

class Rel(UserExpression):
    def __init__(self,materials,**kwargs):
        super().__init__(**kwargs)
        self._materials = materials
    def eval_cell(self, values, x, cell):
        if self._materials[cell.index] == 0:
            values[0] = RAir
        else:
            values[0] = RCu

   
Reluc = Rel(materials,degree=0)

class Jd(UserExpression):
    def __init__(self,materials,**kwargs):
        super().__init__(**kwargs)
        self.materials = materials
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = JAir
        else:
            values[0] = JCu


Jz = Jd(materials,degree=0)

# Boundary Conditions
A_bound = Constant(0.0)
bc = DirichletBC(V,A_bound,boundary)

# Defining the Function
Ax = Function(V)
v = TestFunction(V)
F = inner(grad(Ax*Reluc),grad(v))*dx 
Ax = Function(V)
solve(F==0,Ax,bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-6}})

Ay = Function(V)
F = inner(grad(Ay*Reluc),grad(v))*dx
Ay = Function(V)
solve(F==0,Ay,bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-6}})

Az = Function(V)
F = inner(grad(Az*Reluc),grad(v))*dx - Jz*v*dx
Az = Function(V)
solve(F==0,Az,bc,solver_parameters={"newton_solver":{"relative_tolerance":1e-6}})

# Plotting
Vp = VectorFunctionSpace(mesh_tot,'P',1)
A_tot = Function(V)
A_tot = Expression(('Ax','Ay','Az'),Ax=Ax,Ay=Ay,Az=Az)
vtkfile = File('Vectorpotential.pvd')
vtkfile<<A_to

There are three Poisson equations here that needs to be solved.The first two were solved
,but the third Poisson equation shows a RuntimeError with regards to PETsc,which I am not able to understand.

fenics@d2cb9f28de4f:~$ python3 ~/shared/Cylinder/cylinder.py
Generating mesh with CGAL 3D mesh generator
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
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.
Calling FFC just-in-time (JIT) compiler, this may take some time.
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) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-06)
Newton solver finished in 0 iterations and 0 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) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-06)
Newton solver finished in 0 iterations and 0 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) = 1.060e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Traceback (most recent call last):
File “/home/fenics/shared/Cylinder/cylinder.py”, line 92, in
solve(F==0,Az,bc,solver_parameters={“newton_solver”:{“relative_tolerance”:1e-6}})
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-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 /tmp/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

fenics@d2cb9f28de4f:~$

Can somebody please help me in understanding the solution?

2

Please format your code using ``` such that it runnable by others without having to do manual indentation and formatting

3

I am sorry for the inconvenience.I have made the edit as required.

4

There are several things in your code you need to have a look at.

as far as I can tell, the solution to this equation is 0. The same holds for the second equation.

Here is a working version of your code using linear solvers (and not user expressions, as I do not like them).

from fenics import *
from mshr import *
from math import *
# Create Geometry`
Rair = 5
cyltot = Cylinder(Point(0,0,10),Point(0,0,0),Rair,Rair)
mesh_tot = generate_mesh(cyltot,64)
vtkfile = File('/home/fenics/shared/Cylinder/Mesh_tot.pvd')
vtkfile << mesh_tot
V = FunctionSpace(mesh_tot,'P',1)


# Subdomains

import numpy as np
tol = 1E-14

class Copper(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) <= 1 + tol


class Air(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) >= 1 - tol


materials = MeshFunction('size_t',mesh_tot,mesh_tot.topology().dim())
Subdom_cu = Copper()
Subdom_air = Air()
materials.set_all(0)
Subdom_cu.mark(materials, 1)

# Defining the parameters
RCu = 1/1.256637E-6
RAir = 1/1.256629E-6
JCu = 6/(np.pi)
JAir = 0

Q = FunctionSpace(mesh_tot, "DG", 0)
Jz = Function(Q)
x = Q.tabulate_dof_coordinates()
for i in range(x.shape[0]):
    if materials.array()[i] == 0:
        Jz.vector().vec().setValueLocal(i, JAir)
    else:

        Jz.vector().vec().setValueLocal(i, JCu)

Reluc = Function(Q)
for i in range(x.shape[0]):
    if materials.array()[i] == 0:
        Reluc.vector().vec().setValueLocal(i, RAir)
    else:
        Reluc.vector().vec().setValueLocal(i, RCu)

# Boundary Conditions
A_bound = Constant(0.0)
bc = DirichletBC(V,A_bound,"on_boundary")

# Defining the Function
Ax = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(Ax*Reluc),grad(v))*dx
l = inner(Constant(0), v)*dx
Ax = Function(V)
solve(a==l,Ax,bc)
Ay = TrialFunction(V)
a2 = inner(grad(Ay*Reluc),grad(v))*dx
Ay = Function(V)
solve(a2==l,Ay)
Az = TrialFunction(V)
a3 = inner(grad(Az*Reluc),grad(v))*dx - Jz*v*dx
l3 = rhs(a3)
Az = Function(V)
solve(lhs(a3)==l3,Az,bc)
5

Thanks! The linear solver worked.