Hello.
I’m running some code on a relatively fine mesh (approx 800000 cells) and I’m getting the following error:
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|| = -1.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 2
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 7e99df3564aea9dc05729ddb64498ae47b6bc15d
*** -------------------------------------------------------------------------
Now I understand from reading previous posts that this is related to the linear solvers and preconditioners I’m choosing for my problem. I have tried using different linear solvers / preconditioners and also changing the tolerances and nothing seems to work so I’d appreciate some help on this matter. I’ll attach below my code for reference:
import matplotlib.pyplot as plt
from fenics import *
import numpy as np
from tqdm import tqdm
import math
import ufl as uf
import csv
import os
from mpi4py import MPI
from lcf_funcs import *
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
handle = 'axi_lcf'
handle1 = 'axi_lcf'
if (rank == 0 and os.path.isfile(f'{handle1}/data.csv') == True):
os.remove(f'{handle1}/data.csv')
set_log_active(False)
parameters["std_out_all_processes"] = False;
parameters['ghost_mode'] = 'shared_facet'
parameters["mesh_partitioner"] = "ParMETIS"
file_phi = XDMFFile(f'{handle1}/phi.xdmf')
file_u = XDMFFile(f'{handle1}/u.xdmf')
file_p = XDMFFile(f'{handle1}/p.xdmf')
file_tau = XDMFFile(f'{handle1}/tau.xdmf')
file_phi.parameters['flush_output'] = True
file_u.parameters['flush_output'] = True
file_p.parameters['flush_output'] = True
file_tau.parameters['flush_output'] = True
rho_in = 1.204
eta_s_in = 0.00001825
eta_p_in = 0
rho_out = 1000.9
eta_s_out = 0.031
eta_p_out = 1.511
lamb_in = 0
lamb_out = 0.207
sig = 0.076
grav = 9.81
nx=100
ny=200
T=0.2
num_steps = 8000
qdiv=40
num_steps_rein = 3
rein_div=1
dt= T/num_steps
g = Constant((0, grav))
scaling = 0.1
toler = 1.0e-3
centrex=0
centrey=0.025
radius = 0.00288
x0 = 0
y0 = 0
x1 = 0.05
y1 = 0.1
mesh = RectangleMesh(Point(x0,y0),Point(x1,y1),nx,ny,'left')
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
for facet in facets(cell):
for vertex in vertices(facet):
if (0 <= vertex.point().array()[0] <= 0.02) and (0 <= vertex.point().array()[1] <= 0.075):
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
for facet in facets(cell):
for vertex in vertices(facet):
if (0 <= vertex.point().array()[0] <= 0.02) and (0 <= vertex.point().array()[1] <= 0.075):
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
for facet in facets(cell):
for vertex in vertices(facet):
if (0 <= vertex.point().array()[0] <= 0.02) and (0 <= vertex.point().array()[1] <= 0.075):
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
for facet in facets(cell):
for vertex in vertices(facet):
if (0 <= vertex.point().array()[0] <= 0.02) and (0 <= vertex.point().array()[1] <= 0.075):
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
h = 0.05*CellDiameter(mesh)
eps = Constant(2.5*mesh.hmin())
eps1 = Constant(2.5*mesh.hmin())
dtau = Constant(0.2*mesh.hmin())
dist = Expression('sqrt((pow((x[0]-A),2))+(pow((x[1]-B),2)))-r', degree=2, A=centrex, B=centrey, r=radius)
dist2 = Expression('(1/(1+exp((dist/eps))))',degree=2, eps=eps, dist=dist)
Qs = FunctionSpace(mesh, 'CG', 2)
phi = TrialFunction(Qs)
psi = TestFunction(Qs)
phi00 = Function(Qs)
phi0 = interpolate(dist2, Qs)
phi_rein = Function(Qs)
Vnormal = VectorFunctionSpace(mesh, 'CG', 1)
phigrad = Function(Vnormal)
vnorm = TestFunction(Vnormal)
Ps = FunctionSpace(mesh, 'CG', 1)
p = TrialFunction(Ps)
q = TestFunction(Ps)
p0 = Function(Ps)
p_ = Function(Ps)
Vs = VectorFunctionSpace(mesh, 'CG', 2)
u = TrialFunction(Vs)
v = TestFunction(Vs)
u0 = Function(Vs)
u_ = Function(Vs)
Ts = TensorFunctionSpace(mesh, 'CG', 2)
tau = TrialFunction(Ts)
S = TestFunction(Ts)
tau0 = Function(Ts)
Ttt = FunctionSpace(mesh, 'CG', 2)
tautt = TrialFunction(Ttt)
zetatt = TestFunction(Ttt)
tau0tt = Function(Ttt)
T_devssg = TensorFunctionSpace(mesh, 'DG', 2)
G = TrialFunction(T_devssg)
R = TestFunction(T_devssg)
G1 = Function(T_devssg)
T_tt_devssg = FunctionSpace(mesh, 'DG', 2)
G_tt = TrialFunction(T_tt_devssg)
R_tt = TestFunction(T_tt_devssg)
G1_tt = Function(T_tt_devssg)
walls = f'near(x[1], {y0}) || near(x[1], {y1})'
fswalls = f'near(x[0], {x0}) || near(x[0], {x1})'
bcu_noslip = DirichletBC(Vs, Constant((0, 0)), walls)
bcu_fslip = DirichletBC(Vs.sub(0), Constant(0), fswalls)
def flux(ui, ni, taui):
return (dot(ui,ni)*taui + abs(dot(ui,ni))*taui)/2
n=FacetNormal(mesh)
x = SpatialCoordinate(mesh)
r=x[0]
bcu = [bcu_noslip, bcu_fslip]
k = r
F_levelset = (phi/dt)*psi*k*dx - (phi0/dt)*psi*k*dx + inner(u0,grad(phi))*psi*k*dx \
+ ((phi - phi0)/dt + inner(u0, grad(phi))) \
* scaling*mesh.hmin()/uf.Max(2.0*sqrt(inner(u0, u0)),toler/mesh.hmin())*inner(u0, grad(psi))*k*dx
F_grad = inner((phigrad-ngamma(phi0)),vnorm)*k*dx
phi = Function(Qs)
"""new rein"""
F_rein = ((phi_rein - phi0)/dtau)*psi*k*dx \
- 4*phi_rein*(1.0 - phi_rein)*inner(grad(psi), phigrad)*k*dx \
+ eps1*inner(grad(phi_rein), grad(psi))*k*dx
supg = (h/magnitude(u0)+0.000001)
# F_devssg = inner(G - grad(u0),R)*r*dx
# F_tt_devssg = inner(G_tt - u0[0]/r,R_tt)*r*dx
if handle == 'axi_lcf':
tau_axi_pre = as_tensor([[tau0[0,0], tau0[0,1], 0],
[tau0[1,0], tau0[1,1] , 0],
[0, 0, tau0tt]])
F_stress = inner((1/dt)*(tau-tau0),S)*r*dx \
+ inner(reaction(tau_axi_pre, tau0, u0, lamb_out),S)*r*dx \
+ 0.5*inner(dot(u0,nabla_grad(tau)),S)*r*dx \
+ 0.5*inner(dot(u0,nabla_grad(tau0)),S)*r*dx \
+ 0.5*inner(dot(u0, nabla_grad(S)), supg*dot(u0,nabla_grad(tau)))*r*dx \
+ 0.5*inner(dot(u0, nabla_grad(S)), supg*dot(u0,nabla_grad(tau0)))*r*dx \
Ftt = (1/dt)*(tautt-tau0tt)*zetatt*r*dx \
+ reactt(tau_axi_pre, tau0, u0, lamb_out, r)*zetatt*r*dx \
+ 0.5*(u0[1]*tautt.dx(1) + u0[0]*tautt.dx(0))*zetatt*r*dx \
+ 0.5*(u0[1]*tau0tt.dx(1) + u0[0]*tau0tt.dx(0))*zetatt*r*dx \
+ 0.5*inner((u0[1]*zetatt.dx(1) + u0[0]*zetatt.dx(0)), supg*(u0[1]*tautt.dx(1) + u0[0]*tautt.dx(0)))*r*dx \
+ 0.5*inner((u0[1]*zetatt.dx(1) + u0[0]*zetatt.dx(0)), supg*(u0[1]*tau0tt.dx(1) + u0[0]*tau0tt.dx(0)))*r*dx \
elif handle == 'axi_norm':
ide = as_matrix([[1,0],[0,1]])
tau0.assign(project(ide,Ts))
gie = 0.1
F_stress = inner((1/dt)*(tau-tau0),S)*r*dx \
- inner(dot(tau0, trans(grad(u0))),S)*r*dx \
- inner(dot(grad(u0), tau0),S)*r*dx \
+ (1/lamb_out)*inner((exp(0.05*(tau0[0,0]+tau0[1,1]+tau0tt - 3)))*(tau0-Identity(2)),S)*r*dx \
+ 0.5*inner(dot(u0,nabla_grad(tau)), S)*r*dx \
+ 0.5*inner(dot(u0,nabla_grad(tau0)), S)*r*dx \
+ 0.5*inner(dot(u0, nabla_grad(S)), supg*dot(u0,nabla_grad(tau)))*r*dx \
+ 0.5*inner(dot(u0, nabla_grad(S)), supg*dot(u0,nabla_grad(tau0)))*r*dx \
ide_tt = 1
tau0tt.assign(project(ide_tt,Ttt))
Ftt = (1/dt)*(tautt-tau0tt)*zetatt*r*dx \
- 2*tau0tt*(u0[0]/r)*zetatt*r*dx \
+ (1/lamb_out)*inner((exp(0.05*(tau0[0,0]+tau0[1,1]+tau0tt - 3)))*(tau0tt-1),zetatt)*r*dx \
+ 0.5*(u0[1]*tautt.dx(1) + u0[0]*tautt.dx(0))*zetatt*r*dx \
+ 0.5*(u0[1]*tau0tt.dx(1) + u0[0]*tau0tt.dx(0))*zetatt*r*dx \
+ 0.5*inner((u0[1]*zetatt.dx(1) + u0[0]*zetatt.dx(0)), supg*(u0[1]*tautt.dx(1) + u0[0]*tautt.dx(0)))*r*dx \
+ 0.5*inner((u0[1]*zetatt.dx(1) + u0[0]*zetatt.dx(0)), supg*(u0[1]*tau0tt.dx(1) + u0[0]*tau0tt.dx(0)))*r*dx \
outn = as_vector([ngamma(phi0)[0],0,ngamma(phi0)[1]])
out = outer(outn, outn)
# G3D = as_tensor([[G1[0,0],0, G1[0,1]],
# [0,G1_tt,0],
# [G1[1,0],0,G1[1,1]]])
ns1 = (1/dt)*inner((rho(phi)*u - rho(phi)*u0), v)*r*dx \
+ inner(dot(rho(phi)*u0, nabla_grad(u)), v)*r*dx \
+ 2*inner(eta_s(phi)*epsilon_axi(u,x),epsilon_axi(v,x))*r*dx \
- p0*div_axi(v,x)*r*dx \
+ inner(rho(phi)*g,v)*r*dx \
+ sig*mgrad(phi0)*inner((Identity(3) - out), epsilon_axi(v,x))*r*dx \
# + 2*inner((1-eta_s(phi))*epsilon_axi(u0,x),epsilon_axi(v,x))*r*dx \
# - inner((1-eta_s(phi))*(G3D+trans(G3D)),epsilon_axi(v,x))*r*dx \
if handle == 'axi_lcf':
tau_axi = vmat3d(tau0)*as_matrix([[exp(e1(tau0)),0 ,0],
[0,exp(e2(tau0)) ,0],
[0, 0, exp(e3(tau_axi_pre,tau0))]])*vmat3d(tau0).T
ns1 += (eta_p(phi)/lamb_out)*inner(tau_axi - Identity(3), nga_flip(v, x))*r*dx \
elif handle == 'axi_norm':
tau_axi = as_tensor([[tau0[0,0], 0, tau0[0,1]],
[0, tau0tt, 0],
[tau0[1,0], 0, tau0[1,1]]])
ns1 += (eta_p(phi)/lamb_out)*inner(tau_axi - Identity(3), nabla_grad_axi(v, x))*r*dx \
ns2 = (1/rho(phi))*(dot(nabla_grad(p),nabla_grad(q)) - dot(nabla_grad(p0),nabla_grad(q)))*r*dx \
+ (1/dt)*div_axi(u_,x)*q*r*dx \
ns3 = inner(u,v)*r*dx - inner(u_,v)*r*dx \
+ (dt/rho(phi))*inner(nabla_grad(p_-p0),v)*r*dx \
bcs=[]
tol = 1.0e-4
tau =Function(Ts)
tautt =Function(Ttt)
t=0
q=0
# for k in range(int(8)):
# a = assemble(lhs(F_levelset))
# L = assemble(rhs(F_levelset))
# solve(a, phi.vector(), L, 'gmres', 'default')
# phi0.assign(phi)
for n in tqdm(range(num_steps)):
t+=dt
q+=1
a = assemble(lhs(F_levelset))
L = assemble(rhs(F_levelset))
solve(a, phi.vector(), L, 'gmres', 'default')
if n % rein_div == 0:
phi0.assign(phi)
solve(F_grad == 0, phigrad ,solver_parameters={"newton_solver": {"linear_solver": 'gmres', "preconditioner": 'default',\
"maximum_iterations": 200, "absolute_tolerance": 1e-25, "relative_tolerance": 1e-25}}, \
form_compiler_parameters={"optimize": True})
for n in range(num_steps_rein):
solve(F_rein == 0, phi_rein ,solver_parameters={"newton_solver": {"linear_solver": 'gmres', "preconditioner": 'default',\
"maximum_iterations": 20, "absolute_tolerance": 1e-8, "relative_tolerance": 1e-6}}, \
form_compiler_parameters={"optimize": True})
phi00.assign(phi0)
phi0.assign(phi_rein)
phi.assign(phi_rein)
else:
phi0.assign(phi)
A1=assemble(lhs(ns1))
M1=assemble(rhs(ns1))
[bc.apply(A1) for bc in bcu]
[bc.apply(M1) for bc in bcu]
solve(A1, u_.vector(), M1, 'gmres', 'default')
# u_, G_, G__tt = split(w_)
M2=assemble(rhs(ns2))
A2=assemble(lhs(ns2))
solve(A2, p_.vector(), M2, 'gmres', 'default')
A3=assemble(lhs(ns3))
M3=assemble(rhs(ns3))
[bc.apply(A3) for bc in bcu]
[bc.apply(M3) for bc in bcu]
solve(A3, u_.vector(), M3, 'cg', 'sor')
u0.assign(u_)
p0.assign(p_)
# A_devssg = assemble(lhs(F_devssg))
# M_devssg = assemble(rhs(F_devssg))
# solve(A_devssg, G1.vector(), M_devssg, 'gmres', 'default')
# A_tt_devssg = assemble(lhs(F_tt_devssg))
# M_tt_devssg = assemble(rhs(F_tt_devssg))
# solve(A_tt_devssg, G1_tt.vector(), M_tt_devssg, 'gmres', 'default')
A = assemble(lhs(F_stress))
M = assemble(rhs(F_stress))
solve(A, tau.vector(), M, 'gmres', 'default')
Att = assemble(lhs(Ftt))
Mtt = assemble(rhs(Ftt))
solve(Att, tautt.vector(), Mtt, 'gmres', 'default')
tau0.assign(tau)
tau0tt.assign(tautt)
if q % qdiv == 0:
file_phi.write(phi, t)
file_u.write(u0, t)
file_tau.write(tau0, t)
phi0.assign(phi)
area = assemble(phi0*r*dx)
area_x2 = 2*area
x_com = assemble(Expression("x[0]", degree = 1)*phi0*dx)/area
y_com = assemble(Expression("x[1]", degree = 1)*phi0*dx)/area
Pa = 2.0*sqrt(np.pi*area_x2)
Pb = 2*assemble(mgrad(phi0)*r*dx)
circ = Pa/Pb
u_rise = assemble(u0[0]*phi0*r*dx)/area
v_rise = assemble(u0[1]*phi0*r*dx)/area
timeseries = [t, area_x2, x_com, y_com, circ, u_rise, v_rise]
if rank == 0:
with open((f"{handle1}/data.csv"), 'a') as csvfile:
f = csv.writer(csvfile, delimiter='\t',lineterminator='\n',)
f.writerow(timeseries)
The line where the error is being caused is the following:
solve(F_grad == 0, phigrad ,solver_parameters={"newton_solver": {"linear_solver": 'gmres', "preconditioner": 'default',\
"maximum_iterations": 200, "absolute_tolerance": 1e-25, "relative_tolerance": 1e-25}}, \
form_compiler_parameters={"optimize": True})
and no time iterations are completed, the code fails instantly.
Thanks in advance. Will.