Howdy,
Just curious - I’m talking to the PETSC folks, and shared this error message
*** -------------------------------------------------------------------------
*** 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: f3de96341c54e434c48115de5a3975b482aee933
*** -------------------------------
Here is an example code that causes the isssue:
from dolfin import *
import math
import collections
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -1.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -1.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class zBottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -1.0)
class zTop(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (-(5.0/8.0), (5.0/8.0))) and between(x[0], (-(5.0/8.0), (5.0/8.0))) and between(x[2], (-(5.0/8.0), (5.0/8.0))))
class Obstacletwo(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (-.5, .5)) and between(x[1], (.7, 1.0)) and between(x[2], (.7, 1.0)))
def direction_three(psi,theta):
return (math.sin(psi)*math.cos(theta),math.sin(psi)*math.sin(theta),math.cos(psi))
class ObstacleFour(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (.7, 1.0)) and between(x[1], (.7, 1.0)) and between(x[1], (.7, 1.0)))
class Circle(SubDomain):
def inside(self, x, on_boundary):
return (x[0]-.85)**2 + (x[1]-.85)**2 + (x[2]-.85)**2 < pow(.1,2)
set_log_level(LogLevel.ERROR)
#set_log_active(False)
left = Left()
top = Top()
right = Right()
bottom = Bottom()
topz = zTop()
bottomz=zBottom()
obstacle = Obstacle()
obstacletwo = Obstacletwo()
obstaclefour= ObstacleFour()
circle = Circle()
weights=[1.0]
p1 = Point(-1.0,-1.0,-1.0)
p2 = Point(1.0,1.0,1.0)
size=128
mesh = BoxMesh(p1,p2,size,size,size)
print("building function space")
element = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V = FunctionSpace(mesh, element)
objfunc=FunctionSpace(mesh, "CG", 1)
dg = FunctionSpace(mesh,'DG',0)
domains = MeshFunction("size_t", mesh,mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
obstacletwo.mark(domains, 3)
obstaclefour.mark(domains,4)
circle.mark(domains,5)
domain_number=5
dx = Measure("dx")(subdomain_data=domains)
kosquared= Constant(36*pi*pi)
ko= Constant(6*pi)
q= Constant(.75)
psi = math.pi/4.0
theta = math.pi/4.0
direction = direction_three(psi,theta)
w=Constant(1.0)
k0=6.0*math.pi
exp = str(k0)+"*(" + str(direction[0]) + "*x[0]+" + str(direction[1]) + "*x[1]+" + str(direction[2]) + "*x[2])"
exp = "6*pi*(" + str(direction[0]) + "*x[0] + " + str(direction[1]) + "*x[1])"
uo_r = Expression("36*pi*pi*cos(" +str(exp) + ")",degree=2)
uo_i = Expression("36*pi*pi*sin(" +str(exp) +")",degree=2)
cos_f = Expression("cos("+str(exp) + ")",degree=2)
sin_f = Expression("sin("+str(exp) + ")",degree=2)
u=Function(V)
(u_r,u_i) = split(u)
print(type(u))
print((type(u_r),type(u_i)))
tests = TestFunction(V)
v_r = tests[0]
v_i =tests[1]
print("solving state")
state=inner(grad(u_r),grad(v_r))*dx + inner(grad(u_i),grad(v_i))*dx +ko*u_i*v_r*ds -ko*u_r*v_i*ds -v_r*(kosquared*(1+q*w)*u_r +kosquared*q*w*cos_f)*dx -v_i*(kosquared*(1+q*w)*u_i +kosquared*q*w*sin_f)*dx
Jac = derivative(state, u, TrialFunction(V))
solve(state == 0, u, J=Jac,solver_parameters={'newton_solver': {'linear_solver': 'mumps'}})
print("saved state")
#u_list.append((u.sub(0,True),u.sub(1,True)))
all_states = split(u)
u_r = all_states[0]
u_i = all_states[1]
#print((type(u_r),type(u_i)))
#print(collections.Counter(u.sub(0,True).vector()[:]))
#print(collections.Counter(u.sub(1,True).vector()[:]))
lolz = u_r*dx
print(assemble(lolz))
theta_counter=0
multiplers = Function(V)
tests = TestFunction(V)
v_r = tests[0]
v_i =tests[1]
print("solving adjoint")
#print(weights[theta_counter])
#print(domain_number)
(lambda_1,lambda_2) = split(multiplers)
weight_con = Constant(weights[theta_counter])
#print(weight_con)
adjoint_eq = weight_con*(u_r - cos_f)*v_r*dx(domain_number) + weight_con*(u_i - sin_f)*v_i*dx(domain_number) +inner(grad(v_r),grad(lambda_1))*dx +inner(grad(v_i),grad(lambda_2))*dx - lambda_1*kosquared*(1.0+q*w)*v_r*dx - lambda_2*kosquared*(1+q*w)*v_i*dx - ko*v_r*lambda_2*ds + ko*v_i*lambda_1*ds
t1 = weight_con*(u_r - cos_f)*dx(domain_number) + weight_con*(u_i - sin_f)*dx(domain_number)
#print(assemble(t1))
Jac = derivative(adjoint_eq, multiplers, TrialFunction(V))
solve(adjoint_eq == 0, multiplers, J=Jac,solver_parameters={'newton_solver': {'linear_solver': 'mumps'}})
print("solved adjoint")
(lambda_1,lambda_2) = split(multiplers)
#print(collections.Counter(multiplers.sub(0,True).vector()[:]))
#print(collections.Counter(multiplers.sub(1,True).vector()[:]))
They told me that this is dolfin generated and not the “true” petsc error message - I was trying to solve a 3d pde with 128^3 mesh - could this possibly be a running out of memory issue (this is what they are thinking)
Thanks