Hi,

I am using a nonlinear solver NonlinearVariationalSolver(problem) and honestly I don’t understand why it works for a certain scheme. Could someone point out maybe what I am misunderstanding?

Let’s say both u, v and p, q belong to a degree 1 Lagrangian space, so their basis are [\phi_1 = 1, \phi_2 = x, \phi_3 = y],

u = \alpha_1 + \alpha_2 x + \alpha_3 y, and p = \alpha_4 + \alpha_5 x + \alpha_6 y,

v = \alpha'_1 + \alpha'_2 x + \alpha'_3 y, and q = \alpha'_4 + \alpha'_5 x + \alpha'_6 y

The scheme (which contains other terms not including v) depends on v exclusively through

p * v.dx(1) * dx,

i.e. we do not have v in the scheme, only its derivative wrt one variable. Then two rows of the Jacobian J = derivative(F, …) will be fully zeros:

F_1 = p * \alpha’_1 * \phi_1.dx(1) = 0, and

F_2 = p * \alpha’_2 * \phi_2.dx(1) = 0,

as both \phi_1 = 1 and \phi_2 = x disappear as we differentiate them wrt y. So any partial derivatives of these zeros will be zero.

Thus the Jacobian according to this has at least 2 fully zero lines / columns. How is it possible that the Newton solver still works?

I am using http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html,

line (28) to understand how the Jacobian is created.

Thanks a lot for any hints!