Simulation works in 1D but not in 2D

I think this may be a similar issue to my answer here. Try setting the absolute_tolerance for the Newton solver to 0:

    solve(F==0, u1, bc, solver_parameters={"newton_solver":{"absolute_tolerance":0.0}})