[firedrake] difference between matrix-action and form-action

[Lawrence Mitchell]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 27/04/15 13:14, Eike Mueller wrote:
> Hi,
> 
> I think I just saw a rounding error, and the results are the same.
> I hadn’t realised that the numbers I compare are very large, so
> that an absolute difference of 1.E-5 is rounding error. I will
> simply normalise the numbers by the maximum before doing the
> comparison.

Aha, I think this will depend on what your boundary condition /value/
is.  Here's a bit more explanation:

Consider splitting the function space V = V_0 \osum V_g

where V_0 are the basis functions that vanish on the boundary nodes
and V_g are the basis functions that vanish on the interior.

So when I assemble my operator it splits into:

[ A_00, A_0g;
  A_g0, A_gg ]

Now, A_g0 is zero (due to choice of basis functions) and A_gg is the
identity, giving:


[ A_00, A_0g;
  0   , I    ]

When we assemble the operator for solving a system of equations, we
/also/ throw away A_0g giving:

[ A_00, 0;
  0   , I ]

And forward substitute A_0g onto the right hand side, since it acts on
the boundary condition nodes, for which we know the correct solution
value.

However, consider some (known) Function u, we wish to compute the
action of the operator on u.  If we have the assembled operator (with
bcs) then we end up with:


Au = A_00 u_0 + I u_g

where u_g is the value on the bc nodes, and u_0 on the complement.

So to have Au satisfy the boundary conditions we can /either/ apply
the bcs before or after, since the application of the bc operator is
the identity.

If you're computing the action:

assemble(action(a, u), bcs=bcs)

Then something different happens:

Au|bcs = (A_00 u_0 + A_0g u_0 + A_gg u_g)|bcs

So we first assemble the action of the operator /without/ bcs and then
apply the bc values.  Now the answer will be correct on the bc nodes
(it will contain the bc values) but will contain an additional A_0g
u_0 term on some of the non-bc nodes (that overlap with the bc basis
functions).

To get the correct answer, we need to kill the A_0g u_0 term, to do
this, we can apply homogeneous bcs to u, compute the action, and then
apply the correct bc values to the result.  Like so:

hbc = homogenize(bc)

hbc.apply(u)
Au = assemble(action(a, u), bcs=bc)

Cheers,

Lawrence
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1

iQEcBAEBAgAGBQJVPjEHAAoJECOc1kQ8PEYvDEUIAL6+30Xer/JQsPJpHf+UjCQS
ZbaIdfoRc63Gd7WyQr3MVqqgQpOqyNABEWDXunS7vWsLIV1l+5lNHtkhG0Fe58NO
L0skC+tAetOiIRswev1FZgrhvIBcLcvYePqbhjr7bCGGZGEYfGajErwv/87iydx9
fQL1TO27bAAH3CRaQHJ1mjLOYsfw1PeEpmh+lFQtKpCDKyuvNeQfGtV49kGYl869
2ERZJngUBaVlBIqLaO67mS094Un+6XHAeTcD3Ok5D20Q8pKd4vGxUT7wuXRSP2ji
ZhDZY40xuAXvpXQBQbGxoOy59lQzxMF0GSLpCxrtB4L+0C97SMf5N8lnPsOlY3o=
=Yxo8
-----END PGP SIGNATURE-----