Bug in newtonsolver for mixed-dimensional problem

[Ingeborg Gjerde]

Hi,

The following code gives different results depending on if i run it using ceciledc/fenics_mixed_dimensional:21-06-22 or ceciledc/fenics_mixed_dimensional:16-01-23. The updated version does not converge.

If I use the older version the Newton solver converges within 6 iterations at each mesh refinement and the errors are decrease as the mesh is refined.

If I use the newer version the Newton solver does not converge.

from dolfin import *

errors =[]
for N in [2, 4, 8, 16]:

    # Create mesh
    mesh2 = UnitSquareMesh(N,N)

    # Create analytic solution
    x, y = SpatialCoordinate(mesh2)
    u2a_2 = sin(pi*x)*sin(pi*y)

    f = -div((1+u2a_2**2)*grad(u2a_2))

    # Define boundaries
    marker = MeshFunction("size_t", mesh2, mesh2.topology().dim() - 1, 0)

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

    b_left = BoundaryLeft()
    b_left.mark(marker, 1)
    mesh1 = MeshView.create(marker, 1)

    x, y = SpatialCoordinate(mesh1)
    u2a_1 = sin(pi*x)*sin(pi*y)


    dx2 = Measure("dx", domain=mesh2)
    dx1 = Measure("dx", domain=mesh1)


    # Define function space
    V2 = FunctionSpace(mesh2, 'CG', 1)
    V1 = FunctionSpace(mesh1, 'CG', 1)

    W = MixedFunctionSpace(V2, V1)


    # Define trial and test functions
    sol = Function(W)
    u2, u1 = sol.split()

    v2, v1 = TestFunctions(W)


    # Weak form
    a = inner((1+u2**2)*grad(u2) , grad(v2)) * dx2 
    a += inner(u1, v2) * dx1
    a += inner(u2, v1) * dx1 


    L = + inner(v2, project(f, V2)) * dx2
    L += inner(v1, project(u2a_1, V1)) * dx1

    R = a-L

    solve(R==0, sol)    
    v2_sol, v1_sol = sol.split()

    u2a_2 = project(u2a_2, V2)
    errors.append(errornorm(u2a_2, v2_sol))

print(errors)

This issue was first flagged in https://fenicsproject.discourse.group/t/mixed-dimensional-branch-problem-with-time-integration/10428/12. I changed the test case to correspond to the nonlinear Poisson demo with the boundary conditions imposed using a lower-dimensional Lagrange multiplier.

[Ingeborg Gjerde]

To expand a bit, using the newer image causes the Newton solver to bounce around

Solving mixed nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.470e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 9.055e-01 (tol = 1.000e-10) r (rel) = 2.026e-01 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 3.753e-01 (tol = 1.000e-10) r (rel) = 8.396e-02 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 3.284e-01 (tol = 1.000e-10) r (rel) = 7.347e-02 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 4.897e-01 (tol = 1.000e-10) r (rel) = 1.096e-01 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 1.179e+00 (tol = 1.000e-10) r (rel) = 2.638e-01 (tol = 1.000e-09)
  Newton iteration 6: r (abs) = 2.702e+00 (tol = 1.000e-10) r (rel) = 6.045e-01 (tol = 1.000e-09)

My best guess is that this is due to a faulty function evaluation on the lower-dimensional mesh, I reported this here.