[deal.II] FEValues pre-computed values

[Doug]

Hello again,

I am currently looking into the implementation FEValues. Here is my 
understanding

   1. FEValues takes some Mapping, FiniteElement, and Quadrature to 
   evaluate and store the values with the required UpdateFlags.
   2. FEValues::initialize through the FiniteElement::get_data will then 
   pre-compute all the required data that are possibly computable on the 
   reference cell.
   3. FEValues::reinit will then compute/fetch the data on a physical cell, 
   using the Jacobian information by calling the various fill_fe_values() 
   functions.

I see that most of the examples use implicit methods to solve the various 
problems. Therefore, the cost of initialize() and reinit() on each cell is 
quite low compared to solving the actual system at every step. However, for 
explicit methods such as RK, we wouldn't want to re-compute the Jacobian or 
reference gradients (\nabla_{physical} \phi = J^{-T} \nabla _{reference} 
\phi) at every time step. Therefore, the FEValues declaration would be 
outside the "assemble_system()" seen in most tutorials, and passed as an 
argument as "assemble_system(FEValues fe_v)". Unfortunately, I am getting 
lost in the code where things values are pre-computed versus fetching them.

My question is: What is being re-computed every time I call 
FEValues::reinit(cell)?

I am mainly worried about it re-computing Jacobians (its determinant and 
its inverse) and evaluating its product with the reference gradients 
(\nabla_{reference} phi) every time reinit(cell) is called. Otherwise, the 
solution will be to pre-compute and store those outside in some other 
container and avoid using FEValues.

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.

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.