[firedrake] Pass approximate Jacobian to derivative() call

[Lawrence Mitchell]

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