Current distribution problem

Hello everyone! I want to simulate the current distribution. I use the continuity equation and the Neumann boundary conditions.

The current is set at the boundary as a function of potential. However, the task diverges as soon as the total current is not zero. What condition needs to be added to fix the sum of the currents?

(abs(int1d(Th, 1)(ISteel) + int1d(Th, 2)(IMg)) < 1e-1)

func ISteel = 14225 * (-phi)^6 + 93891 * (-phi)^5 + 249128 * (-phi)^4 + 341543 * (-phi)^3 + 255613 * (-phi)^2 + 99090 * (-phi) + 15552;
func IMg = -1764.2 * (-phi)^4 - 10108 * (-phi)^3 - 20773 * (-phi)^2 - 17653 * (-phi) - 4910.6; 

problem Current(phi, v) =  
  int2d(Th)(Sigma * grad(phi)' * grad(v)) 
  - int1d(Th, 1)(ISteel * v)
  - int1d(Th, 2)(IMg * v);

Hello,
About your nonlinear equation, it is well-posed if f (or both f on the two pieces of boundary) is increasing in the range of \phi that you expect. Then it can be solved by the Newton method. If f is not everywhere increasing, it does not mean that the problem is ill-posed, but it could be difficult to find a stable scheme to solve it. Eventually you could add to f a stabilizing term c*(\phi^{k+1}-\phi^k) where k is the iteration index, and c is a sufficiently large constant (so that c+f'\geq 0), or in the Newton method replace f'(\phi^k) by \max(0,f'(\phi^k)).
Anyway it seems that if f is ISteel or IMg, the sign is wrong in your int1d