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!