petsc4py.PETSc.Error: error code 55

Errors
1

Hello,
I have a problem similar to the one in Error setting PETSc options repeatedly. I want to solve a non-linear problem with an iterative scheme, where the system matrix (not the right-hand side as in the similar problem) changes in each iteration. After a couple of iterations, I obtain the error

  problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 594, in __init__
    opts[k] = v
  File "PETSc/Options.pyx", line 23, in petsc4py.PETSc.Options.__setitem__
  File "PETSc/Options.pyx", line 91, in petsc4py.PETSc.Options.setValue
petsc4py.PETSc.Error: error code 55

I think I should not initialize a new LinearSolver in each iteration, but how do I update the bilinear form a otherwise? Here is some code that leads to the error:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 16, 16, cell_type=dolfinx.mesh.CellType.triangle)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
u_h = dolfinx.fem.Function(V)
p = dolfinx.fem.Constant(mesh,1.5)
f = dolfinx.fem.Constant(mesh,1.)
    
def u_exact(x):
    return x[0]*x[1]
u_h.interpolate(u_exact)

def KacanovStep(u_h):
    from ufl import conditional          
    #from petsc4py import PETSc  
    def ufl_max(a, b):
        return conditional(a > b, a, b)   
    def ufl_min(a, b):
        return conditional(a < b, a, b)
    V = u_h.function_space
    mesh = V.mesh
    eps_min_sq = dolfinx.fem.Constant(mesh,0.01)
    eps_max_sq = dolfinx.fem.Constant(mesh,100.)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    # Define bilinear form and right-hand side
    a = (ufl_min(eps_max_sq,ufl_max(eps_min_sq,ufl.dot(ufl.grad(u_h),ufl.grad(u_h))))**(p/2-1)*ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx     
    L = f * v * ufl.dx
    # Apply homogeneous Dirichlet boundary conditions
    def boundary(x):
        return np.full((1, x.shape[1]), True)
    boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1,boundary)
    bc = dolfinx.fem.dirichletbc(u_h, dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets))        
    problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    u_h = problem.solve()
    return u_h
    
for curIter in range(0,1000):
    u_h = KacanovStep(u_h) 
    print('nrIter =', curIter,'J(u_h) =',dolfinx.fem.assemble_scalar(dolfinx.fem.form((1/p*ufl.inner(ufl.grad(u_h),ufl.grad(u_h))**(p/2))*ufl.dx)))
2

You can have a look at: Error setting PETSc options repeatedly - #2 by nate
as I currently do not have time to parse your code in detail. Note that you can create a separate function for the previous solution, and assign data to it after each iteration (similar to any time stepping demo), to avoid recompilation. See for instance the dolfinx tutorial.

3

Thanks for the response. I adjusted the Code, now it works:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
from ufl import conditional   

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 16, 16, cell_type=dolfinx.mesh.CellType.triangle)
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
u_h = dolfinx.fem.Function(V)
p = dolfinx.fem.Constant(mesh,1.5)
f = dolfinx.fem.Constant(mesh,1.)
    
u_h.interpolate(lambda x: x[0]*x[1])

# KacanovStep:       
def ufl_max(a, b):
    return conditional(a > b, a, b)   
def ufl_min(a, b):
    return conditional(a < b, a, b)
V = u_h.function_space
mesh = V.mesh
eps_min_sq = dolfinx.fem.Constant(mesh,0.01)
eps_max_sq = dolfinx.fem.Constant(mesh,100.)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define bilinear form and right-hand side
a = (ufl_min(eps_max_sq,ufl_max(eps_min_sq,ufl.dot(ufl.grad(u_h),ufl.grad(u_h))))**(p/2-1)*ufl.dot(ufl.grad(u), ufl.grad(v))) * ufl.dx     
L = f * v * ufl.dx
# Apply homogeneous Dirichlet boundary conditions
def boundary(x):
    return np.full((1, x.shape[1]), True)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1,boundary)
bc = dolfinx.fem.dirichletbc(u_h, dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets))        
KacanovStep = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    
for curIter in range(0,1000):
    u_new = KacanovStep.solve()
    u_h.vector.array = u_new.vector.array
    print('nrIter =', curIter,'J(u_h) =',dolfinx.fem.assemble_scalar(dolfinx.fem.form((1/p*ufl.inner(ufl.grad(u_h),ufl.grad(u_h))**(p/2))*ufl.dx)))

Best regards,
Johannes