[dune-pdelab] Element-Jacobians by Automatic Differentiation

[Christian Waluga]

Hi,

On 08/10/2012 02:39 PM, Peter Bastian wrote:
> this is really cool. I thought about this from the beginning. We have a
> Master student (computer science) here who is interested to do this as 
> a master thesis (I think he wants to start in a few months). 
> So there is definitively interest in having this in PDELab (!) and
> we should see that the student takes off from where you are...
> BTW, which automatic differentiation software do you work with?
I didn't want to use code transformation, because PDELab is heavily
templated and I'm not really sure if today's differentiation software
can handle that. I think if somebody wants to generate code for
Jacobians, then it should be done e.g. in the Fenics-FFC and interfaced
to PDELab. If I remember correctly, Steffen has experimented in this
direction.

Currently I use a simple class which implements the AD by operator
overloading, i.e., I have a type that stores function values (e.g. v)
and their derivatives (vector dv). Everytime some operation is applied
to the variable, the derivatives are also updated, e.g. operator *=
would do the following (pseudocode):

AutoDiff& operator *= (AutoDiff const& rhs)
{
   for(int i = 0; i < rhs.dv.size(); ++i)
     dv[i] = dv[i] * rhs.v + v * rhs.dv[i];
   v *= rhs.v;
}

I could have also used a library for this but I figured out that this
would be overkill for demonstration purposes. However, I had to make
some sacrifices concerning the performance of the code in order to
incorporate this into PDELab. Also the implementation of the AD was done
in a quick and dirty fashion, so either this has to be improved or one
should consider interfacing a stable and optimized library for this
(e.g. from the list here:
http://en.wikipedia.org/wiki/Automatic_differentiation, but I have no
experience with any of those).

The part where the AD comes into contact with PDELab then is very
similar to the Numerical-Jacobian backend. Except now, one has to hand
over a coefficient vector of AutoDiffs to the alpha-Routines with
appropriate 'seeds' for the derivatives. The residual vector also
contains elements of type AutoDiff and while computing the alpha-terms,
the code automatically also computes the Jacobian.

Meanwhile I have added a download of the code on my homepage (section
"Software" at bottom)
https://www-m2.ma.tum.de/bin/view/Allgemeines/ChristianWalugaEN

The relevant code is in the directory dune-ad/dune/pdelab/localoperator.
There you can find the implementation of the Jacobian-Backend, the
AD-Datatype and the (slightly) modified Poisson-Localoperator.

Best, Christian