Automatic Differentiation of the Finite Element Jacobian

[Doug]

Hello,

We had chosen to use deal.ii since most of it looked templated in terms of accepting the float type. However, I noticed that the FEValuesBase returns a 'double', instead of a generic 'value_type' for the Jacobian term.

Our goal would be to differentiate through the entire residual (and possibly the time step) with respect to the solution variables and mesh points. Therefore, the general element Jacobian (and its determinant) would need to be differentiated. I see that the Jacobian entries already have differentiation functions, however, that would entail manually assembling the derivative of the residual with respect to the mesh point. We would like to instead have a function that assembles the residual right-hand side and pass the AD types to that high-level function to obtain dRdX.

Looking at the current implementation in fe_values, mapping_q_generic, and local integrators, I see that it looks relatively easy to template the value type since it is already partially templated.

Is there a reason why those functions were not templated? Is there anything we should be careful about before we start templating those functions?

Doug

[Jean-Paul Pelteret]

Dear Doug,

Looking at the current implementation in fe_values, mapping_q_generic, and local integrators, I see that it looks relatively easy to template the value type since it is already partially templated.

Is there a reason why those functions were not templates?

No, its likely that a more generic use-case was just never considered — this was previously found to be the case from the “pre-AD supported” implementation of FEValues and FEValuesViews classes. We’d be happy to accept any patches that adds AD support to the classes that you find to be lacking it. 

Is there anything we should be careful about before we start templating those functions?

Yes, off the top of my head I can think of at least one interesting issue that you may face. There are occasional checks / shortcuts that are used to speed up computations (using floating-point types), but would invalidate AD-based evaluation. See one such example here. You’d need to check carefully for similar cases and make sure that the entire computational graph is traversed in order for the result of differentiation to be correct.

Best,

Jean-Paul

[Wolfgang Bangerth]

Doug,

in addition to what J-P already wrote:

> We had chosen to use deal.ii since most of it looked templated in terms of
> accepting the float type. However, I noticed that the FEValuesBase returns a
> 'double', instead of a generic 'value_type' for the Jacobian term.

That depends on the context.

Shape functions are indeed defined as phi(x) returning simply a double. That's
because of an underlying assumption that x is a location parameterized as
doubles. So functions such as shape_value(), shape_grad(), etc really just
return doubles. (This includes the Jacobian of the shape functions, i.e.,
their gradients.) In other words, it's currently not possible to do
autodifferentiation of the shape functions with regard to the position x.

But functions such as get_function_gradient() which compute quantities such as

u_h(x) = \sum_j U_j phi(x_j)

are templatized on the data type of the vector U of unknowns. So they should
allow computing Jacobians of u_h with regard to the coefficients U_j. This is
what one needs to do to compute residual vectors from energy functionals, and
Hessian matrices from residual vectors.

The latter functionality is implemented and routinely used. In other words...

> Our goal would be to differentiate through the entire residual (and possibly

> the time step) with respect to the solution variables and *mesh points*.

...the former of this works, whereas the latter probably does not.


As J-P already mentioned, the reason for all of this is not a conscious design
decision, but that nobody who needed this in the past has implemented it. We'd
be quite happy to accept patches in this direction. It is often useful if you
posted an example of what you want to do and discussed implementation
strategies before starting the work -- this makes merging this functionality
later on easier. I see that you've read through a lot of code already, though,
so you'll have a good sense of what it'll take.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Doug]

Thank you both for the amazingly quick reply.

We will very likely implement this work at some point in our development in a few months. We have quite a bit of other basics to implement first since we have just started reading through the documentation. We plan on using the deal.ii library as the workhorse behind an hp-adaptive, parallel flow solver used for aerodynamic shape optimization. The derivatives with respect to the mesh points will be needed for optimization purposes.

I see that in a previous thread (https://groups.google.com/forum/#!searchin/dealii/tpetra%7Csort:date/dealii/I8aLPu5anrM/gBzaaqivDAAJ), you plan on switching over to Tpetra, which will also be useful for us since we might instantiate those vectors with AD types. Do you have an idea of the time horizon on which this will be developped?

Doug

[Wolfgang Bangerth]

On 2/8/19 1:55 PM, Doug wrote:
>
> We will very likely implement this work at some point in our development in a
> few months. We have quite a bit of other basics to implement first since we
> have just started reading through the documentation. We plan on using the
> deal.ii library as the workhorse behind an hp-adaptive, parallel flow solver
> used for aerodynamic shape optimization. The derivatives with respect to the
> mesh points will be needed for optimization purposes.

Ah, you want to do shape optimization. I understand now why you need these
derivatives.

> I see that in a previous thread
> (https://groups.google.com/forum/#!searchin/dealii/tpetra%7Csort:date/dealii/I8aLPu5anrM/gBzaaqivDAAJ),
> you plan on switching over to Tpetra, which will also be useful for us since
> we might instantiate those vectors with AD types. Do you have an idea of the
> time horizon on which this will be developped?

There is a patch here:
https://github.com/dealii/dealii/pull/7180
I suspect that it will be merged in the next few weeks or months, but the
question is what one can do with these vectors. I'll let Daniel Arndt answer
that question, but suspect that additional work is necessary to use them with
linear solvers, multiply with matrices, etc.

You should already be able to instantiate the other deal.II vectors with AD
types, though!

[Jean-Paul Pelteret]

The derivatives with respect to the 

mesh points will be needed for optimization purposes.

I’m not familiar with the algorithms used in shape optimisation, but one *possible* work-around if you need derivatives wrt. the mesh points on the finite element cell level would be to do the following:

- Add an extra component to an FESystem, specifically a FE_Q(1)^dim. You would interpret as the current position of the mesh vertices because the support points for this FE coincide with the grid vertices.

- You can initialise this field using VectorTools::interpolate. Specifically you would simply populate the relevant entries of the result returned by this Function with the coordinates given by the input Point<dim>.

- Now that you have these “extra DoFs”, you can now compute the sensitivities on any cell with respect to them.

This, of course, would not really help if you compute some objective function based on the entire solution, and you in fact need your entire program to de differentiable.

Best,

Jean-Paul

[Daniel Arndt]

> I see that in a previous thread
> (https://groups.google.com/forum/#!searchin/dealii/tpetra%7Csort:date/dealii/I8aLPu5anrM/gBzaaqivDAAJ),
> you plan on switching over to Tpetra, which will also be useful for us since
> we might instantiate those vectors with AD types. Do you have an idea of the
> time horizon on which this will be developped?

There is a patch here:
https://github.com/dealii/dealii/pull/7180
I suspect that it will be merged in the next few weeks or months, but the
question is what one can do with these vectors. I'll let Daniel Arndt answer
that question, but suspect that additional work is necessary to use them with
linear solvers, multiply with matrices, etc.

The pull request mentioned above just introduces a wrapper class for Tpetra vectors

and is a first step towards Tpetra support. This vector class can be used with TrilinosWrappers::SparseMatrix if

the Scalar template parameter is double. The next step is to implement a Tpetra::CrsMatrix wrapper

that is templated with respect to the scalar type as well. Using these wrappers with deal.II`s iterative solvers should then work out-of-the-box.

For preconditioners, we likely want to use wrappers for Trilinos classes as well, though. This will likely require some more work.

Of course, every contribution in that direction is welcome!

Best,

Daniel

[Wolfgang Bangerth]

Doug,
let me add another thought I had: While the vertex locations of the mesh are
hard-coded as doubles (and you won't be able to change that), there is a class
MappingEulerian that allows you to provide a displacement field. This field is
parameterized by a vector of some sort, and it is possible that you can make
it work if that field is a vector of some AD type for which you can thread the
remainder of the computations through AD evaluations.

I don't know whether that's going to work out any easier. but I did want to
point it out as an option worth reading through.