Modifying diagonal entries of Jacobian (preconditioning)
Dear Firedrakers, if I wanted to modify the diagonal entries of my Jacobian in a Newton step during a time step... How would I do this? Do I have to define my own preconditioner here? Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de ------------------------------------------------------
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Dear Henrik, On 09/10/15 14:16, Buesing, Henrik wrote:
Dear Firedrakers,
if I wanted to modify the diagonal entries of my Jacobian in a Newton step during a time step... How would I do this? Do I have to define my own preconditioner here?
Do you actually want to modify the Jacobian, or do you only want to precondition the Newton step with an operator that is different from the true Jacobian? In the former case, you'll be doing some kind of quasi-Newton. Recall the Newton update is: du = -J^{-1} F(u) (1) where J is the linearisation of F about u (the Jacobian). If you directly modify J, then you're solving for a different update: du = -J'^{-1} F(u) (2) Which may or may not converge to the solution you want, depending on the choice of J'. If you want to do this, then when you build your solver you can explicitly provide a bilinear form for the Jacobian (rather than asking firedrake to derive it) by saying: solve(F == 0, u, J=bilinear_form, ...) Alternately, you may be wanting to solve (1), i.e. the problem with the exact Jacobian, but you have a choice of how to approximate J^{-1}. The system: J du = -F(u) is solved for du using an iterative method, by default, the preconditioner for this is constructed using J, effectively solving: P(J)^{-1} J du = -P(J)^{-1} F(u) However, you need not use J to construct the preconditioner, you can instead provide an approximate Jacobian that captures the structure of your problem but may be "easier" to invert, call it Jp. Then you end up solving: P(Jp)^{-1} J du = -P(Jp)^{-1} F(u) To do this, you can provide a bilinear form for Firedrake to construct Jp like so: solve(F == 0, u, Jp=approximate_jacobian, ...) Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWJMWhAAoJECOc1kQ8PEYvMGoH/iSoNwXTlgZL1BmPj5jr9QYB BWG8JsUn+v0yO/CHHx7MY1HRFELjtzUCyfmLHzOQcJYb7pYJIiE9UtKtMKYktR1g lo9O/GKuw8oIEED9usx+jFtw5tZvV9KB45JyXtOHrPY9jroRubXC6o5ZTNtlCRyG G85/zKFmnRZB+sF7PFkFUjxewr95WpQ+E/O4uvyuII4FfUYMwVTmxNc7mg6iXJnd 0LKwf7etueG1VnCkTnNA8vPjPzhRbltpXFoZJy3tLzhRRkjFjdKBhDGUz7ZjWXrj 35ozjxfC1jSG8VgH8IkaN0BDy4NhCw4o85eZLMCsAQwXSWbJIcTRVVDSLYKBikw= =w1/i -----END PGP SIGNATURE-----
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell Gesendet: 19 October 2015 12:28 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Modifying diagonal entries of Jacobian (preconditioning)
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Dear Henrik,
On 09/10/15 14:16, Buesing, Henrik wrote:
Dear Firedrakers,
if I wanted to modify the diagonal entries of my Jacobian in a Newton step during a time step... How would I do this? Do I have to define my own preconditioner here?
Do you actually want to modify the Jacobian, or do you only want to precondition the Newton step with an operator that is different from the true Jacobian?
[Buesing, Henrik] Well, in certain problems my formulation has zero entries in the diagonal of the Jacobian. What I did in my Fortran Code was just adding a small value to these diagonal entries, and then solve for a different update as you say. This solved problems with LU decomposition. I understand that I can just provide a different J for the solve: NonlinearVariationalProblem(F, u, J= bilinear_form, ...) But how would I get the modified J? I can calculate J and assemble it, then change the values in the assembled matrix J = derivative(F,u) M = assemble(J).M for i in range(0,2) : for j in range(0,2) : print "Block: ", i,j print M[i,j].values # change M[i,j].values[l,k] If I then just do solver.solve() in every timestep ... is the modified J used? Or do I have to call NonlinearVariationalProblem in every timestep? At the moment I only change constants t and dt via assign. But I also want to be able to change coefficients, e.g. changing density with pressure/temperature... And sometimes during a phase change, even a coefficient turns into a primary unknown... The former I can realize with Constant() and assign, I guess. But for the latter, do I have to reassemble the form? Thank you! Henrik
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWJMWhAAoJECOc1kQ8PEYvMGoH/iSoNwXTlgZL1BmPj5jr9QYB BWG8JsUn+v0yO/CHHx7MY1HRFELjtzUCyfmLHzOQcJYb7pYJIiE9UtKtMKYktR1g lo9O/GKuw8oIEED9usx+jFtw5tZvV9KB45JyXtOHrPY9jroRubXC6o5ZTNtlCRyG G85/zKFmnRZB+sF7PFkFUjxewr95WpQ+E/O4uvyuII4FfUYMwVTmxNc7mg6iXJnd 0LKwf7etueG1VnCkTnNA8vPjPzhRbltpXFoZJy3tLzhRRkjFjdKBhDGUz7ZjWXrj 35ozjxfC1jSG8VgH8IkaN0BDy4NhCw4o85eZLMCsAQwXSWbJIcTRVVDSLYKBikw= =w1/i -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 27/10/15 08:51, Buesing, Henrik wrote:
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell Gesendet: 19 October 2015 12:28 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Modifying diagonal entries of Jacobian (preconditioning)
Dear Henrik,
On 09/10/15 14:16, Buesing, Henrik wrote:
Dear Firedrakers,
if I wanted to modify the diagonal entries of my Jacobian in a Newton step during a time step... How would I do this? Do I have to define my own preconditioner here?
Do you actually want to modify the Jacobian, or do you only want to precondition the Newton step with an operator that is different from the true Jacobian?
[Buesing, Henrik]
Well, in certain problems my formulation has zero entries in the diagonal of the Jacobian. What I did in my Fortran Code was just adding a small value to these diagonal entries, and then solve for a different update as you say. This solved problems with LU decomposition.
If you're using a factorisation-based preconditioner (e.g. LU, Cholesky, Incomplete LU (or Incomplete Cholesky)), then I think you can ask PETSc to shift the zero diagonals. Run with the solver paramete r: 'pc_factor_shift_amount': epsilon To shift the matrix diagonal by epsilon.
I understand that I can just provide a different J for the solve:
NonlinearVariationalProblem(F, u, J= bilinear_form, ...)
But how would I get the modified J? I can calculate J and assemble it, then change the values in the assembled matrix
Under the hood what happens is that the PETSc SNES (nonlinear solver object) gets a callback to firedrake to assemble the Jacobian at each Newton step. This callback is written to just assemble the form and not do anything else. However, in theory you can override it to do something different. In practice, this isn't part of the exposed API, so it's a little tricky. If you really need to do this, we can point in the right direction.
J = derivative(F,u) M = assemble(J).M for i in range(0,2) : for j in range(0,2) : print "Block: ", i,j print M[i,j].values # change M[i,j].values[l,k]
FWIW, this won't change the values at all (because the values array is a dense (copied) view of the matrix).
If I then just do solver.solve() in every timestep ... is the modified J used?
Or do I have to call NonlinearVariationalProblem in every timestep? At the moment I only change constants t and dt via assign.
But I also want to be able to change coefficients, e.g. changing density with pressure/temperature... And sometimes during a phase change, even a coefficient turns into a primary unknown...
So we don't have good facilities for modifying the system in between Newton steps. If the only thing that changes in the form are the values of the coefficients, you can just write things in the natural way. However, you can't take a newton step, then (say) update an out-of-band equation of state and take another newton step easily. If you just need to update the coefficients in between solves, the that's easy.
The former I can realize with Constant() and assign, I guess. But for the latter, do I have to reassemble the form?
You need to do somewhat more than that, since you're solving for a new unknown, you'll need to use a different form entirely, no? Or do I misunderstand? Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWL0mVAAoJECOc1kQ8PEYvMEMIAOt1HRpkc+1HmUGhw+xpnVIC BT0B7tgrSRGj3i8/hR1sPKPMJwNFWiiNFITK2KSe2N6uqrvWmH17cOKvru0FAs8L OfWAdOFkeWHCl2eBuuVMyEnCy9LOO5SLwttTl6jePLFM+ROgh+/4VOkPhPtJ6Bjg 3JwOl15ACIWQSK2zxbjJvyZIYSVjAyKFpMPJmLEXkKdDkkCeouIREnfD/tWMrynz hQabk0RHJ6hg3JNPz3vq/aC7Gv7P0KDkVwSTtWWsuYYjfGkkR7WQ3wxGfFtG0C0M OVHE12HP0S1ibPS709nCMdgeWSqjHFoVSKgHCOBeUfG9JyFc6V1zZ75UJzwu3S4= =2xqt -----END PGP SIGNATURE-----
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell Gesendet: 27 October 2015 10:53 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Modifying diagonal entries of Jacobian (preconditioning)
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
On 27/10/15 08:51, Buesing, Henrik wrote:
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell Gesendet: 19 October 2015 12:28 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Modifying diagonal entries of Jacobian (preconditioning)
Dear Henrik,
On 09/10/15 14:16, Buesing, Henrik wrote:
Dear Firedrakers,
if I wanted to modify the diagonal entries of my Jacobian in a Newton step during a time step... How would I do this? Do I have to define my own preconditioner here?
Do you actually want to modify the Jacobian, or do you only want to precondition the Newton step with an operator that is different from the true Jacobian?
[Buesing, Henrik]
Well, in certain problems my formulation has zero entries in the diagonal of the Jacobian. What I did in my Fortran Code was just adding a small value to these diagonal entries, and then solve for a different update as you say. This solved problems with LU decomposition.
If you're using a factorisation-based preconditioner (e.g. LU, Cholesky, Incomplete LU (or Incomplete Cholesky)), then I think you can ask PETSc to shift the zero diagonals. Run with the solver paramete r:
'pc_factor_shift_amount': epsilon
To shift the matrix diagonal by epsilon. [Buesing, Henrik] Yes, I didn't know this option. That's exactly what I need. I will use that one.
But I also want to be able to change coefficients, e.g. changing density with pressure/temperature... And sometimes during a phase change, even a coefficient turns into a primary unknown...
So we don't have good facilities for modifying the system in between Newton steps. If the only thing that changes in the form are the values of the coefficients, you can just write things in the natural way. However, you can't take a newton step, then (say) update an out-of-band equation of state and take another newton step easily. If you just need to update the coefficients in between solves, the that's easy.
The former I can realize with Constant() and assign, I guess. But for the latter, do I have to reassemble the form?
You need to do somewhat more than that, since you're solving for a new unknown, you'll need to use a different form entirely, no? Or do I misunderstand?
[Buesing, Henrik] I think, I need a new form yes. The unknowns stay the same. I'm solving for pressure and enthalpy, whether in the two-phase or single-phase region. But in the single-phase region one of the coefficients turns into an unknown.
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWL0mVAAoJECOc1kQ8PEYvMEMIAOt1HRpkc+1HmUGhw+xpnVIC BT0B7tgrSRGj3i8/hR1sPKPMJwNFWiiNFITK2KSe2N6uqrvWmH17cOKvru0FAs8L OfWAdOFkeWHCl2eBuuVMyEnCy9LOO5SLwttTl6jePLFM+ROgh+/4VOkPhPtJ6Bjg 3JwOl15ACIWQSK2zxbjJvyZIYSVjAyKFpMPJmLEXkKdDkkCeouIREnfD/tWMrynz hQabk0RHJ6hg3JNPz3vq/aC7Gv7P0KDkVwSTtWWsuYYjfGkkR7WQ3wxGfFtG0C0M OVHE12HP0S1ibPS709nCMdgeWSqjHFoVSKgHCOBeUfG9JyFc6V1zZ75UJzwu3S4= =2xqt -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
So we don't have good facilities for modifying the system in between Newton steps. If the only thing that changes in the form are the values of the coefficients, you can just write things in the natural
way.
However, you can't take a newton step, then (say) update an out-of-band equation of state and take another newton step easily. If you just need to update the coefficients in between solves, the that's easy.
The former I can realize with Constant() and assign, I guess. But for the latter, do I have to reassemble the form?
You need to do somewhat more than that, since you're solving for a new unknown, you'll need to use a different form entirely, no? Or do I misunderstand?
[Buesing, Henrik] I think, I need a new form yes. The unknowns stay the same. I'm solving for pressure and enthalpy, whether in the two-phase or single-phase region. But in the single-phase region one of the coefficients turns into an unknown.
[Buesing, Henrik] So, can I have a coefficient c(p,h), which is a nonlinear function of my primary variables and then for some h<=hmin or h>=hmax it turns into c(p,h)=h, where h is an unknown? Henrik
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWL0mVAAoJECOc1kQ8PEYvMEMIAOt1HRpkc+1HmUGhw+xpnVIC BT0B7tgrSRGj3i8/hR1sPKPMJwNFWiiNFITK2KSe2N6uqrvWmH17cOKvru0FAs8L OfWAdOFkeWHCl2eBuuVMyEnCy9LOO5SLwttTl6jePLFM+ROgh+/4VOkPhPtJ6Bjg 3JwOl15ACIWQSK2zxbjJvyZIYSVjAyKFpMPJmLEXkKdDkkCeouIREnfD/tWMrynz hQabk0RHJ6hg3JNPz3vq/aC7Gv7P0KDkVwSTtWWsuYYjfGkkR7WQ3wxGfFtG0C0M OVHE12HP0S1ibPS709nCMdgeWSqjHFoVSKgHCOBeUfG9JyFc6V1zZ75UJzwu3S4= =2xqt -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Buesing, Henrik
-
Lawrence Mitchell