-----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-----