PETSc Krylov Solver RuntimeError with fine mesh

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?

1 Like

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

Hi,dokken
I’m using “Incompressible Navier-Stokes equations for flow around a cylinder” from the fenics tutorial and modifying the program to generate mesh with about 35,000 vertices. I used mpirun for parallel computing, but encountered the following error

*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 3

When using a splitting scheme, you need to satisfy the CFL condition. This is dependent on the mesh size. See for instance:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html

the dolfinx is the newest fenics project,but the version of fenics is v2019.1.0 on my machine.
How should this be resolved in this version?

You would need to bring your CFL number down.

For instance using an implicit scheme helps.
Also reducing the time step size helps. See Courant–Friedrichs–Lewy condition - Wikipedia