Why is my program not converging for this system of coupled PDEs?

variational formulation
1

Dear all,

I am trying to solve the following coupled system of PDEs using the mixed-dimensional branch:

- \Delta_\Gamma u = v|_\Gamma - 2u + g on \Gamma = \partial \Omega,

-\nabla v \cdot n = v|_\Gamma - u on \Gamma,

- \Delta v = 0 in \Omega = \{x \in \mathbb{R}^2 | \lVert x \rVert < 1\}

With g = 17x_1x_2, this system has the classical unique solution (u, v) with u(x) = 3x_1x_2 and v(x) = x_1x_2.

The weak form is given by

\int_\Gamma \nabla_\Gamma u_h \nabla_\Gamma \varphi_{u_h} d\sigma = \int_\Gamma (v|_\Gamma -u)\varphi_{u_h} + (-u + g) \varphi_{u_h} d\sigma,

\int_\Omega \nabla v_h \nabla_\Gamma \varphi_{v_h} dx = - \int_\Gamma (v|_\Gamma -u)\varphi_{v_h} d\sigma.

I was able to solve both of the decoupled systems. When i try to solve the coupled system, however, the solver does not converge, even with the exact solution as initial guess.

Any help is greatly appreciated!

Minimum working example:

from dolfin import *

# Exact solutions
exact_v = Expression("x[0]*x[1]", degree=2)
exact_u = Expression("3*x[0]*x[1]", degree=2)


class Bound(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


def coupled_bulk_surface():

    mesh = Mesh("unitcircle.xml")
    sub_domains = MeshFunction("size_t", mesh, "unitcircle_physical_region.xml")
    boundary_parts = MeshFunction("size_t", mesh, "unitcircle_facet_region.xml")
        
    marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
    marker.set_all(0)
    bound = Bound()
    bound.mark(marker, 1)

    # Create submesh 
    submesh = MeshView.create(marker, 1)

    # Setup formulation 
    V = FunctionSpace(mesh, "CG", 1)
    V_b = FunctionSpace(submesh, "CG", 1)
    W = MixedFunctionSpace(V, V_b)
    
    # Define expressions used in variational forms
    g = Expression('17*x[0]*x[1]', degree=2)    
    z = Function(W)
    v = z.sub(0)
    u = z.sub(1)
    (p, k) = TestFunctions(W)    
            
    # Define measures
    dV = Measure("dx", domain=W.sub_space(0).mesh())
    dL = Measure("dx", domain=W.sub_space(1).mesh())
        
    # Variational formulation
    F = dot(grad(u) , grad(k)) * dL + dot(grad(v), grad(p)) * dV - dot((v - u), k) * dL - dot((-u + g), k) * dL + dot((v - u), p) * dL

    # Set exact solution as inital guess
    v0 = interpolate(exact_v, V)
    assign(z.sub(0), v0)
    u0 = interpolate(exact_u, V_b)
    assign(z.sub(1), u0)
   

    # solve(F == 0, z) 
    # not working
    # Error:   Unable to successfully call PETSc function 'KSPSolve'.
    # Reason:  PETSc error code is: 56 (No support for this operation for this object type).

    solve(F == 0, z, solver_parameters={"nonlinear_solver": "newton", "newton_solver": {"linear_solver": "gmres"}})
    
    approx_v = z.sub(0)
    approx_u = z.sub(1)
        
    return (approx_v, approx_u)


# Main function

# Compute solution    
(approx_v, approx_u) = coupled_bulk_surface()

Output:

python3 mwe.py
No Jacobian form specified for nonlinear mixed variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
[problem] create list of residual forms OK
[problem] create list of jacobian forms OK, J_list size =  4
[MixedNonlinearVariationalProblem]::check_forms - To be implemented
Solving mixed nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.179e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "mwe.py", line 69, in <module>
    (approx_v, approx_u) = coupled_bulk_surface()
  File "mwe.py", line 58, in coupled_bulk_surface
    solve(F == 0, z, solver_parameters={"nonlinear_solver": "newton", "newton_solver": {"linear_solver": "gmres"}})
  File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 298, in _solve_varproblem
    solver.solve()
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 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 3.496623e-06).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  d7350469ca8a7a2d3ba3b5efd3cb34a5bea91853
*** -------------------------------------------------------------------------