Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

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

Best
  W.

-- 

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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

[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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

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

Best
  W.

-- 

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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

[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


On Wednesday, February 6, 2019 at 9:22:17 AM UTC-5, Wolfgang Bangerth wrote:
>
>
> 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/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

[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: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Automatic Differentiation of the Finite Element Jacobian

[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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] 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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.