Dolfin error message converted to PETSC form?


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)

    left = Left()
    top = Top()
    right = Right()
    bottom = Bottom()
    topz = zTop()
    obstacle = Obstacle()
    obstacletwo = Obstacletwo()
    obstaclefour= ObstacleFour()
    circle = Circle()


    p1 = Point(-1.0,-1.0,-1.0)
    p2 = Point(1.0,1.0,1.0)
    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())

    obstacle.mark(domains, 1)
    obstacletwo.mark(domains, 3)

    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)

    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_r,u_i) = split(u)
    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")
    all_states = split(u)
    u_r = all_states[0]
    u_i = all_states[1]
    lolz = u_r*dx
    multiplers = Function(V)
    tests = TestFunction(V)
    v_r = tests[0]
    v_i =tests[1]
    print("solving adjoint")
    (lambda_1,lambda_2) = split(multiplers)
    weight_con = Constant(weights[theta_counter])
    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)
    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)

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)


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:

state=inner(grad(u_r),grad(v_r))dx + inner(grad(u_i),grad(v_i))dx +kou_iv_rds -kou_rv_ids -v_r*(kosquared*(1+qw)u_r +kosquaredqwcos_f)dx -v_i(kosquared(1+qw)u_i +kosquaredqw*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”)

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.

List of KSP types
List of PC types