PETSc Krylov Solver RuntimeError with fine mesh

Hi everyone,

I currently have a fine cylindrical mesh that I built in Gmsh of the form:

R = 1.5;
lc  = .15;
Point(1) = {0,0,0,lc};
Point(2) = {R,0,0,lc};
Point(3) = {0,R,0,lc};
Point(4) = {-R,0,0,lc};
Point(5) = {0,-R,0,lc};

Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};

Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};

Extrude {0,0,-8} {
Surface{6};
}

Field[1] = Distance;
Field[1].NodesList = {1};

Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc / 30;
Field[2].LcMax = lc;
Field[2].DistMin = 0.01;
Field[2].DistMax = 4;

Field[3] = Min;
Field[3].FieldsList = {2};
Background Field = 3;

Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;

I have then imported this into my code via

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2) 

But when I run it, I get this errror:

fenics@aef3e8f96120:~/shared$ python3 target3Dmsh.py 
Time 0.0
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 6.335e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)
Traceback (most recent call last):
  File "target3Dmsh.py", line 86, in <module>
    solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
  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 solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

I think that this is an issue with the solver because I am using a very fine mesh, but I am unsure on how I can fix this. My full code is shown:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import as_tensor
from ufl import Index
import math
import mshr

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(20)

# SOLVE THE HEAT TRANSFER PROBLEM FOR A 3D TARGET

#----------------
t_start = 0.0
t_end = 20
nstep = 40
dtime = (t_end - t_start)/ nstep  #time step
#-----------------
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172         #conductivity [W/mm K] 
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.01925         #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag
#alpha = 4.3943e+4 #absorption coefficient (mm^-1)
alpha = 1
#---------------------------------------


pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.25  # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 500 #peak power of laser in W
w = 0.5     #mm
R = 1.5    #mm
area = pi*R*R  #area of target (mm^2)
depth = 8.0 #thickness of target (mm)
volume = area * depth #volume of target
Iden = Pmax / (pi*w*w*volume*rho)    #Intensity per unit mass
absorb_depth = 1/alpha

Laser = Expression('2*A*alpha*Iden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2))',degree=3, A=0.25, alpha=alpha, Iden=Iden,w= 0.56, w2= absorb_depth) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium


# Create mesh
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)    
    
Space = FunctionSpace((mesh), 'P', 1)      #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets, metadata = {'quadrature_degree': 2})  #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells, metadata = {'quadrature_degree': 2})   #volume element

initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space)         # Temperature
V = TestFunction(Space)     # Test Function used for FEA
dT = TrialFunction(Space)   # Temperature Derivative
q0 = Function(VectorSpace)  # heat flux at previous time step
i = Index()
G = as_tensor(T.dx(i), (i))  #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step 

q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da   #final form to solve
Gain = derivative(F, T, dT)    # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system 

file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
	print( "Time", t)
	solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
	file_T << (T,t)
	q_tmp = project(q, VectorSpace)
	q0.assign(q_tmp)
	T0.assign(T)      #change so that T0 is equal to T for the next time step

This means that the residual is zero at the first iteration, which usually indicates that the initial guess solves your problem. This was also reported for another case at the old Q&A forum.

Hmm thats odd, I doubt that the first guess taken would solve this system as this system should rapidly evolve with time (using a coarse mesh it typically took around 4 iterations to reach a solution per time step). Therefore I suspect that perhaps my implementation of the XDMF file into the code has lead to this error?

There is nothing wrong with the input. As far as I can tell, the issue seems to stem from the usage of mumps as the linear solver. If you use “umfpack” or “lu” it seems to be working. Could be memory related as you are running this on a large problem.
See for instance: Dolfin error message converted to PETSC form?

Ok I see, I will also try reducing the degrees of freedom of my mesh in order to make the problem less computationally intensive. I also have tried using ‘lu’ and ‘umfpack’ as the linear solver but I get no iteration and an error message stating ‘out of memory’:

fenics@a6ee3f384b9a:~/shared$ python3 target3Dmsh.py
Time 0.0
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 6.335e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-03)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "target3Dmsh.py", line 86, in <module>
    solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "umfpack", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
  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 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Have you tried to run your program in parallel using mpi? If you are running this on your own laptop, it is limited how big system you will be able to solve.

No, I haven’t tried running it in parallel, mainly because I haven’t ever run a simulation in parallel before. Any recommendations on where to start?

Most FEniCS codes can be executed in parallel with mpi, I.e. mpirun -n 4 python3 file.py to run on four processors.

When I try to run this, I get

Traceback (most recent call last):
  File "target3Dmsh.py", line 70, in <module>
    T0 = interpolate(initial_T, Space)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/interpolation.py", line 71, in interpolate
    Pv.interpolate(v._cpp_object)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 365, in interpolate
    self._cpp_object.interpolate(u)
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 interpolate function into function space.
*** Reason:  Wrong size of vector.
*** Where:   This error was encountered inside FunctionSpace.cpp.
*** Process: 3
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Try using quay.io/fenicsproject/dev:latest as I believe this issue is fixed in the development release of FEniCS.

Additionally, if you have compiled PETSc and SLEPc with MUMPS you can choose that as your linear solver for a original fine mesh, it should handle it pretty well.

I have tried using MUMPS as the linear solver, but this has lead to the errors above.

Hello, is your problem fixed? I meet the same problem and stuck here