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.

1 Like

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

Hi! @xiaodong_zhang @dokken @dokken I have the same problem. Any solution yet?

Without a minimal example reproducing the error it is very hard to help you.

I’ve already suggested above that this issue has been fixed in the latest version of DOLFIN;

This is a sample of my code. I’m trying to see the evolution of a biological system consisting of a growth factor g and fibroblast. My original problem is way more complicated but this is a minimal example of what produces the error.

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

T = 20 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx = ny = 256
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)

# diffusivity contants
D_g = 0.0035 # m^2 per seconds
D_f = 0.0004 #m^2 per seconds

# Decay rates
lan_g = 0.2 
lan_f = 0.025 

# fraction of volume space ocuupied by each specie 
w_g = 1.0
w_f = 1.0

# growth factor production rate
p_g = 1.0


# rate of production of growth factor by 
pg_f = 0.2 # fibroblast

# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_g, v_f= TestFunctions(V)

# for the standard deviation of the gaussian initial condition
rho_1 = 0.1

# Define functions
u = Function(V)

# Initial condtion
u_0 = Expression(('0.0',
                  '0.4*exp(-1*(pow((x[0]), 2) + pow((x[1]), 2))/(2*rho_1*rho_1))'),
                  rho_1 = rho_1, degree=2)

u_n = interpolate(u_0, V)


# Split system functions to access components
g, f = split(u)
g_n, f_n = split(u_n)
# proliferation rate of fibroblast
pf_g = 5.0*g 

# define the Neumann boundary conditions for fibroblast and growth factors
Nbc_g = Constant(0) 
Nbc_f = Constant(0) 

# rate of the strength of adhesion junction
Sff_max = 0.5 

# starting time
t = 0

# extract the initial condtions from the vector space
g_0, f_0 = u_n.split()

plot(g_0, f_0)

# overcrowding term  
def positivity(func):
    cond=conditional(lt(func,0.0),0.0,func)
    return cond

# volume fraction occupied
rho_u = w_g*g + w_f*f 

# The LHS of the trapezoidal rule
a_b = (g - g_n)*v_g*dx + (f - f_n)*v_f*dx 

# growth factor production term
W_fm = pg_f*f


# Defining the adhesion function of fibroblast-fibroblast adhesion junction 'Sff'
EXP_g = exp(1 - 1/(2*g - g**2))
Sff = Sff_max*EXP_g

# Defining the weak and including the taxis term
P_1 = ((Sff_max*(2-2*g)*EXP_g)/((2*g - g**2)**2)) # derivative of Sff wrt. g

grad_P1 = Sff_max*((((2*((2-2*g)**2)*EXP_g)/((1-(1-g)**2)**3)) -((2*EXP_g)/((1-(1-g)**2)**2)) + ((((2 - 2*g)**2)*EXP_g)/((1-(1-g)**2)**4)))*grad(g))

grad_Sff = P_1*grad(g) # Derivative of Sff

# defining fibroblast-fibroblast adhesion term 
Af = 0.5*((Sff*positivity(1-rho_u) - w_f*Sff*f)*grad(f) + (f*positivity(1-rho_u)*P_1 - w_g*Sff*f)*grad(g))

a_b  = (g - g_n)*v_g*dx + (f - f_n)*v_f*dx

# defining the functions
func_g =   dot(D_g*grad(g), grad(v_g))*dx + lan_g*g*v_g*dx - p_g*W_fm*v_g*dx- D_g*Nbc_g*v_g*ds

func_f = dot(D_f*grad(f), grad(v_f))*dx + lan_f*f*v_f*dx - D_f*Nbc_f*v_f*ds - pf_g*f*positivity(1-rho_u)*v_f*dx - dot(grad(f), Af)*v_f*dx 

# Adding the functions together
f_u = func_g + func_f 

# Euler 
F = a_b + dt*f_u

du = TrialFunction(V)
J = derivative(F, u, du)
Problem = NonlinearVariationalProblem(F, u, J = J)
Solver  = NonlinearVariationalSolver(Problem)
Solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'
Solver.parameters['newton_solver']['linear_solver'] = 'mumps'
Solver.parameters['newton_solver']['relative_tolerance'] = 1.e-20
Solver.parameters['newton_solver']['absolute_tolerance'] = 1.e-11
Solver.parameters['newton_solver']['maximum_iterations'] = 1000
vtkfile_g = File('Output/growth_factor.pvd')
vtkfile_f = File('Output/fibroblast.pvd')
vtkfile_m = File('Output/macrophage.pvd')
vtkfile_e = File('Output/ECM.pvd')
# time step

for n in range(num_steps):
	
	# Update current time
	t += dt
	
	# Compute solution
	Solver.solve()
	
	# Extract the solutions from u
	gn_sol, fn_sol = u.split()
    
	vtkfile_g << (gn_sol,t)
	vtkfile_f << (fn_sol,t)
	
	# Update previous solution
	u_n.assign(u)

I think the issue is from the exponential function used to define Sff but I don’t know how to go about it.

@dokken I just posted my sample code above.

This is very far from a minimal example. I also do not get the error message that is in this issue, but:

Solving nonlinear variational problem.
Traceback (most recent call last):
  File "mwe_old.py", line 139, in <module>
    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.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

I would strongly suggest that you try to reduce the complexity of your problem, and have a look at Default absolute tolerance and relative tolerance - #4 by nate