Threadsafe traceless forward mode with 2nd order derivatives

[cgraeser]

We have bindings for AdolC in dune. Unfortunately the lack of thread-safety of the tape-based mode prohibits its use in thread-parallel assembly. On the other hand traceless forward mode with adtl::adouble has several drawbacks that cannot easily be fixed by patching it. As a remedy I implemented my own forward mode class. If there is interest, I'd happily propose it for inclusion in AdolC.

Here's a rough sketch on the approach:

• The core is a class template ADValue<T, maxOrder, dim> for storing a AD-aware value. T is the raw type (fixed to double in adtl::adouble), maxOrder is the order of computed derivatives, dim is the domain dimension.
• Currently maxOrder<=2 is supported (adtl::adouble supports a single derivative only), but the interface is allows for later extension.
• Importantly, the dimension can be a compile time constant, such that all computations happen on the stack.
• By using the special value dynamicDim as template parameter (i.e. ADValue<T, maxOrder, dynamicDim>), one can switch to a dynamic mode as used in adtl::adouble. As far the latter this supports simple allocation (via std::vector) and boost::pool to manage the storage.
• The used approach provides some thread safety guarantees: With compile time size all computations happen on the stack. With dynamic size an can isolate evaluation contexts. Strictly speaking this does not imply that the ADValue is threadsafe, but it allows to do evaluations (with local variables) in concurrent threads without race conditions due to global variables.

In adtl::adouble is roughly equivalent to ADValue<double,1,dynamicDim>.

In my tests (with complicated nested function but small dimension), using a compile time dimension is about twice as fast as adtl::adouble with boost::pool. In another test I used this to compute the 1st and 2nd order derivatives for a Newton step of a quasilinear PDE (minimal surface). Here the performance is exactly the same as if the derivatives are handwritten.

Notice that the above obviously cannot be achieved by patching adtl::adouble. Hence adding this would mean to add a new templated class.

[cgraeser]

Unfortunately I cannot attach a source file. So I'll link to the current version:
https://gitlab.dune-project.org/fufem/dune-fufem/-/blob/feature/forward-ad/dune/fufem/functions/advalue.hh?ref_type=heads.
Notice that this is almost independent of the rest of Dune. The remaining uses of Dune::range() can easily be replaced.

[cgraeser]

I just noticed that nesting ADValue works. So you could even use ADValue<ADValue<T, 1, dim>, 1, dim> to compute second order derivatives or, with one level more, third order ones.

Since this is just a lucky coincidence, it there may be cases where this fails. I will also be less efficient, since one computes several derivatives multiple times. But it's probably worth to be investigated.

[awalther1]

Hi Carsten,

that is indeed very interesting. Thanks a lot for the opportunity to discuss the integration into ADOL-C! There is actually a variant of the traceless mode that allows the calculation of higher order derivatives. So maybe these two approaches could be merged to some extend. Should we coordinate via email how to proceed?

Best

Andrea

[cgraeser]

I also noticed adtl_hov.h today, but did not understand how to use it. Since my approach is templated, it can probably not be integrated into one of the other solutions. Dropping the template would also destroy the performance benefit. But I understand that it should be avoided to have functionally similar solutions in parallel.

Should we coordinate via email how to proceed?
That's probably a good idea. I'll contact you.

[cgraeser]

See #64 for a Dune-free implementation.