Dolfin error message converted to PETSC form?

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

***

***

*** 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")
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(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'}})
(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

Do you still get the error when solving a much smaller problem? E.g: a 8^3 cube?

Hi Nate,

For 32^2 and 64^2 everything is going great.

How many degrees of freedom arise in the 128^3 problem (E.g. `V.dim()`), and do you think using a direct solver such as mumps is appropriate with the hardware you have available?

If not I’d recommend an iterative solver. It’s likely you’re running out of memory.

Hey Nate,

The DOF is 4293378. What would you suggest as a solver (I really have no clue which is the correct solver to call)?

Here is a piece of my code:

Jac = derivative(state, u, TrialFunction(V))
solve(state == 0, u, J=Jac,solver_parameters={‘newton_solver’: {‘linear_solver’: ‘mumps’}})
print(“saved state”)

Looks sesquilinear. Is this the discretisation of the complex valued modified helmholtz equation? Preconditioning may be challenging. My usual naive route is to try `gmres` with `hypre`.

Hi Nate

Yes it is! To do it do I how would I pass it correctly to

`solve(state == 0, u, J=Jac,solver_parameters={‘newton_solver’: {‘linear_solver’: ‘mumps’}})`?

I recommend taking control and setting the options directly rather than using the black-box method `solve()`. See here for example.

Thanks a lot Nate! Because my BC are in the weak form do I just ignore implementing/doing things with boundary conditions?

yes, or just an empty list.