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.

Re: [deal.II] Mistake in Step 16?

[Wolfgang Bangerth]

On 2/5/19 3:29 AM, ky...@math.uh.edu wrote:
> 
> Hi All,
> 
> I was working through example step 16 and I think I found a mistake.
> 
> In|
> LaplaceProblem::setup_system()
> |
> 
> |We have the following|
> 
> ||
> |
> |
> std::set         dirichlet_boundary_ids;
> Functions::ZeroFunction 
> 
>   
>         homogeneous_dirichlet_bc;
> consttypenameFunctionMap::type 
> dirichlet_boundary_functions
> ={{types::boundary_id 
> (0),_dirichlet_bc
>  
> }};
> VectorTools::interpolate_boundary_values 
> (static_cast
>  
> &>(dof_handler),
>                                             dirichlet_boundary_functions,
>                                             constraints);
> // ...
> mg_constrained_dofs.clear();
> mg_constrained_dofs.initialize(dof_handler);
> mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,dirichlet_boundary_ids);
> 
> 
> |
> The problem I see is that dirichlet_boundary_ids is never initialized with 
> anything, and can be checked that before this last line it has size 0.|


Yes, you are correct that this is a bug. I suppose you are working with the 
9.0 release? This has been fixed here already:
   https://github.com/dealii/dealii/pull/6636

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.

[deal.II] Mistake in Step 16?

[kylew]

Hi All,

I was working through example step 16 and I think I found a mistake.

In
LaplaceProblem::setup_system()

We have the following

std::set  dirichlet_boundary_ids;
Functions::ZeroFunction 

 
 homogeneous_dirichlet_bc;
const typename FunctionMap::type 

 
dirichlet_boundary_functions
= { { types::boundary_id 

(0), _dirichlet_bc } };
VectorTools::interpolate_boundary_values 

 
(static_cast 
&>(
dof_handler),
   dirichlet_boundary_functions,
   constraints);
// ...
mg_constrained_dofs.clear();
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, 
dirichlet_boundary_ids);


The problem I see is that dirichlet_boundary_ids is never initialized with 
anything, and can be checked that before this last line it has size 0.
If this was meant to be empty, then the tutorial is misleading when it says 

"The multigrid constraints have to be initialized. They need to know about 
> the boundary values as well, so we pass the dirichlet_boundary here as 
> well."


I believe it should have something like

std::set dirichlet_boundary_ids;
dirichlet_boundary_ids.insert(types::boundary_id(0));
//...


Is this correct, or is it supposed to be empty and simply need a 
placeholder?

-Kyle Williams

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