Pass approximate Jacobian to derivative() call
Dear all, is it possible to just pass the approximate Jacobian to the derivative call? Since it's so complicated in my case to get the exact Jacobian I would like to do sth. like du = [(F(u + h) - F(u - h))/(2*h)] Thank you! 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 On 11/11/15 07:33, Buesing, Henrik wrote:
Dear all,
is it possible to just pass the approximate Jacobian to the derivative call? Since it’s so complicated in my case to get the exact Jacobian I would like to do sth. like
du = [(F(u + h) – F(u - h))/(2*h)]
So you can't do this by passing something to derivative (which instead computes the jacobian for you given a residual). However, you can ask PETSc to compute a finite difference approximation to your Jacobian using colouring, it, however, uses one-sided, rather than centred differences: J = (F(u + h) - F(u))/h if you pass: solver_parameters={'snes_fd_color': True} Then PETSc will compute the Jacobian using finite differencing with colouring. This requires an additional Ncolours residual computations per Newton iteration, where Ncolours depends on the sparsity structure of the matrix: for 3D it will typically be around 20. You can instead of course build an approximate, symbolic, Jacobian "by hand" and provide that to the solver: you'll then end up doing quasi-newton, rather than newton. Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQw/4AAoJECOc1kQ8PEYvF1AH/3XxCtckj5TIP5BPdhnX+IFz G9fXIniTnTa3Jy6jR1gDLhQf/1YTRCt9j9tfTAvCLBQ9yPpLRd6FNpTNHyOsx3dN GX0DWmYzmswDZUe+Awz1CSaxUuu4vn3uPpCcDL6WBS+XcKSdVgfSQ3kI+BRXUYtC VNH0E0842dK0NdBXNFKR+c88QOp2MkdldCtRqgzS9LfAdYUYWYpB7xp+MSURo5QB f1ewh0dk9umagssaIEDLXygh6a+VOG8MuaPXjvRzHK/hRAWsHQppj4f5G/VpPkrv 0sLQeptGQVP1Ob8fWMTIPD754OcGdIpv0bguJr8gECGnR8kL2IsUBzcEUCexngk= =izop -----END PGP SIGNATURE-----
is it possible to just pass the approximate Jacobian to the derivative call? Since it's so complicated in my case to get the exact Jacobian I would like to do sth. like
du = [(F(u + h) - F(u - h))/(2*h)]
So you can't do this by passing something to derivative (which instead computes the jacobian for you given a residual).
[Buesing, Henrik] Yes, true. I mixed this up. I want to do solve(F == 0, u, J=approximate symbolic Jacobian)
You can instead of course build an approximate, symbolic, Jacobian "by hand" and provide that to the solver: you'll then end up doing quasi- newton, rather than newton.
[Buesing, Henrik] That's what I would like to do. In my multi-physics case (see attachement) the Jacobian would have a 2x2 block structure. I'm unsure how to provide this... I tried to write sth. down for the blocks (see test.py), but this is obviously wrong. Thank you for your help! Henrik
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 11/11/15 12:04, Buesing, Henrik wrote:
is it possible to just pass the approximate Jacobian to the derivative
call? Since it's so complicated in my case to get the exact Jacobian I would like to do sth. like
du = [(F(u + h) - F(u - h))/(2*h)]
So you can't do this by passing something to derivative (which instead computes the jacobian for you given a residual).
[Buesing, Henrik] Yes, true. I mixed this up. I want to do
solve(F == 0, u, J=approximate symbolic Jacobian)
You can instead of course build an approximate, symbolic, Jacobian "by hand" and provide that to the solver: you'll then end up doing quasi- newton, rather than newton.
[Buesing, Henrik] That's what I would like to do. In my multi-physics case (see attachement) the Jacobian would have a 2x2 block structure.
I'm unsure how to provide this... I tried to write sth. down for the blocks (see test.py), but this is obviously wrong.
It looks like you getting mixed up between symbolic and actual computation. You look like you want to do the finite differencing sort of by hand? For example here: delta = Constant(1e-8) tmp = pw eps = delta + delta*tmp eps = assemble(eps) pw = tmp + eps val1 = F pw = tmp - eps val2 = F pw = tmp tmp00 = Constant(0.5)*(val1-val2)/eps it looks like you want val1 to be F evaluated at tmp + eps and val2 to be F evaluated a tmp - eps, but in fact, this is just symbolic assignment so both val1 and val2 are just different names for F. To write your approximate symbolic jacobian, you need to write J in terms of functions of the coefficients and test and trial functions (just like you did F). Then it will have the appropriate block structure, I think. I don't quite understand what all the bits are doing, but it looks like, for example the dF/dpw block is just zero (since F doesn't depend on pw?) Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQzU2AAoJECOc1kQ8PEYvpFIIAL6fgCYZnDr6/utEN5ANbfbA jXLPikM9FvStXwryCQOBrdt7iXuOkrz/ACTqjKvGJtFpqfwdJR/MuVLSfB8+rofr Ic4tOQV+wXkWQQa7Z7kRdmH0hhFRia8zFKqXvBoGkOPQNjKfXMaHfMDsgPhu7ehV xsjZt9zsrBzuZQxuY9fOSoABvgrkVk3xSyLaBML+f0rYJTJZ4BLjUbcXuIl9fJEX 940PUx58ZO3mYqkLCwzQQt+vqchptuKZgaaf/Bg1kZVK+ImlfgeVZrXidlTCH8VA bvwLSl387se2z7lQQOeMGH/yW+NeM48biz+W/EPskkcjLE4QTtbu7JBqN7HxcVw= =PWld -----END PGP SIGNATURE-----
For example here:
delta = Constant(1e-8) tmp = pw eps = delta + delta*tmp eps = assemble(eps) pw = tmp + eps val1 = F pw = tmp - eps val2 = F pw = tmp tmp00 = Constant(0.5)*(val1-val2)/eps
it looks like you want val1 to be F evaluated at tmp + eps and val2 to be F evaluated a tmp - eps, but in fact, this is just symbolic assignment so both val1 and val2 are just different names for F.
To write your approximate symbolic jacobian, you need to write J in terms of functions of the coefficients and test and trial functions (just like you did F). Then it will have the appropriate block structure, I think.
[Buesing, Henrik] Ok. Let's assume u=(pw,enth) are my unknowns. F(u)=0 and G(u)=0 are my two pdes formulated as weak forms. Now my Jacobian would have the form J=[dF(u)/dpw, dF(u)/denth dG(u)/dpw, dG(u)/denth] I guess, I will manage to write down the individual blocks in functional form. But how do I write down the Jacobian in terms of these blocks? Or maybe I define H = F + G. Then J = [dH(u)/dpw, dH(u)/d_enth] and I get the appropriate blocks automatically. Or maybe even J = [dH(u)/du], but here I do not know anymore how to write this down in functional form... Is there a multi-physics example where the user supplies a Jacobian to solve()? I think if I knew how this Jacobian should look like in my 2x2 block case, I would manage with the rest. Thank you!
I don't quite understand what all the bits are doing, but it looks like, for example the dF/dpw block is just zero (since F doesn't depend on pw?)
[Buesing, Henrik] This is just a truncated example. In the end there will be non-zeros in every block.
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 11/11/15 13:30, Buesing, Henrik wrote:
For example here:
delta = Constant(1e-8) tmp = pw eps = delta + delta*tmp eps = assemble(eps) pw = tmp + eps val1 = F pw = tmp - eps val2 = F pw = tmp tmp00 = Constant(0.5)*(val1-val2)/eps
it looks like you want val1 to be F evaluated at tmp + eps and val2 to be F evaluated a tmp - eps, but in fact, this is just symbolic assignment so both val1 and val2 are just different names for F.
To write your approximate symbolic jacobian, you need to write J in terms of functions of the coefficients and test and trial functions (just like you did F). Then it will have the appropriate block structure, I think.
[Buesing, Henrik] Ok. Let's assume u=(pw,enth) are my unknowns. F(u)=0 and G(u)=0 are my two pdes formulated as weak forms. Now my Jacobian would have the form
J=[dF(u)/dpw, dF(u)/denth
dG(u)/dpw, dG(u)/denth]
I guess, I will manage to write down the individual blocks in functional form. But how do I write down the Jacobian in terms of these blocks?
Or maybe I define H = F + G. Then J = [dH(u)/dpw, dH(u)/d_enth] and I get the appropriate blocks automatically.
Or maybe even J = [dH(u)/du], but here I do not know anymore how to write this down in functional form...
Is there a multi-physics example where the user supplies a Jacobian to solve()? I think if I knew how this Jacobian should look like in my 2x2 block case, I would manage with the rest.
firedrake/tests/regression/test_vector_laplace_on_quadrilaterals.py does something similar. It solves with one specified Jacobian, and uses a different one as the preconditioner. The idea is that the blocks are "implicit" in the fact that your test and trial functions are from a mixed space. For example: W = DG1*DG2 u, p = TrialFunctions(W) v, q = TestFunctions(W) J = u*v*dx + p*q*dx Builds a jacobian that is, block-wise: J = [DG1-mass-matrix, 0; 0, DG2-mass-matrix] So you just need to write down the symbolic form in terms of things that have come from the mixed space. Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQ0TFAAoJECOc1kQ8PEYvumMH/Ag8LSVLA0CrhYKk5RDjlvGt 9rCx2BpRbFnWWdoj17bFKNpVgTHN6EnTJ05rnwzVo5WI1ox4y4xiX/YA61djeiq4 WafttNE61gspaWenxN4LiYAn8E3bPtNNIJwkwDYr7eiae/xGPxKNLMGM0Ij8QqXZ GoTA7kDwx5RIwrARKskBa42c5HiAZaROAFVWwigUKsXfl0zfkTN6IvTfktaQDRKH Pr72knYSUTgdH/hcx5AVm9Z1GUVeSyzlKDel0U1UaLsNhUU2XHf1eci90mp0xnAi R1bA8VfqvL6GeUjSUchAqIdlaLK5x0aIgLFea0tXUtLubjq/D244yCNT0uvs5go= =gBWk -----END PGP SIGNATURE-----
Is there a multi-physics example where the user supplies a Jacobian to solve()? I think if I knew how this Jacobian should look like in my 2x2 block case, I would manage with the rest.
firedrake/tests/regression/test_vector_laplace_on_quadrilaterals.py does something similar. It solves with one specified Jacobian, and uses a different one as the preconditioner.
The idea is that the blocks are "implicit" in the fact that your test and trial functions are from a mixed space.
For example:
W = DG1*DG2
u, p = TrialFunctions(W) v, q = TestFunctions(W)
J = u*v*dx + p*q*dx
Builds a jacobian that is, block-wise:
J = [DG1-mass-matrix, 0; 0, DG2-mass-matrix]
So you just need to write down the symbolic form in terms of things that have come from the mixed space.
[Buesing, Henrik] So in my nonlinear case I would use u=Function(W) pw, enth=split(u) and write down my J in terms of pw,enth. Nevertheless I will have four contributions (for the four blocks) in my Jacobian, which I just add together as in your example. Every contribution has a (pw + eps) and (pw - eps), or with enth. Can I write down my form as a function dependent on the primary variables. Then I do not need to define all these different +/- eps... Henrik
participants (2)
-
Buesing, Henrik
-
Lawrence Mitchell