Using Sacado with example 15

[Maxi Miller]

I would like to get a feeling in working with Sacado, and would like to try it on example 15. According to my understanding I now have to treat the boundary conditions on my own, and not by using the ConstraintMatrix. Is that still true, or was there an update since 2008?

Furthermore, now the function F() is depending on the gradient of u, and not as in example 33 depending on u, else it should look like exactly in example 33 (in both examples we have \nabla F(u)). Is that correct? Nevertheless I assume that I need the numerical flux function due to having to check the dangling nodes.

Are those assumptions correct? I would like to have the theory first before starting to modify/write programs...

Thanks!

[Wolfgang Bangerth]

> I would like to get a feeling in working with Sacado, and would like to try it
> on example 15. According to my understanding I now have to treat the boundary
> conditions on my own, and not by using the ConstraintMatrix. Is that still
> true, or was there an update since 2008?

I suspect that you are referencing a particular comment in step-33, but it's
been so long that I can't recall what that program does. Can you point out
what you have in mind?

In general, as long as you have Dirichlet boundary conditions, they can be
handled by ConstraintMatrix at the time when the local matrix is distributed
into the global matrix, and that happens after all of the Sacado magic happens
-- you get the local matrix out of Sacado and then copy it into the global matrix.

> Furthermore, now the function F() is depending on the gradient of u, and not
> as in example 33 depending on u, else it should look like exactly in example
> 33 (in both examples we have \nabla F(u)). Is that correct? Nevertheless I
> assume that I need the numerical flux function due to having to check the
> dangling nodes.

You can also handle hanging nodes as is always done, via ConstraintMatrix.
step-33 does it differently not because that is necessary but because the
author decided that he could.

I think that using Sacado (or in general, automatic differentiation) should be
relatively straightforward in step-15, in particular because the residual is
already computed there. I'd just give it a try and see what is necessary. You
may want to use the current developer version of deal.II -- there have been
numerous changes in recent weeks to make AD easier to use in the FEValues
class, in particular through functions such as
FEValues::get_function_{values,gradients}.

Using AD in step-15 is an interesting project. If you make it work, would you
mind sending the assembly function to the mailing list? I think it would make
for a great section in the "Possibilities for extensions" section of step-15!

Best
W.

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

[Praveen C]

I did this a while back, see here

https://bitbucket.org/cpraveen/deal_ii/src/865f1c963f1089a870054bebf7180c9712e93da8/step-15-sacado/?at=master

Best

praveen

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Maxi Miller]

Concerning the constraints: I was refering to this line:

In the first case, there is nothing we need to do: we are using a continuous finite element, and face terms do not appear in the bilinear form in this case. The second case usually does not lead to face terms either if we enforce hanging node constraints strongly (as in all previous tutorial programs so far whenever we used continuous finite elements – this enforcement is done by the ConstraintMatrix class together with DoFTools::make_hanging_node_constraints). In the current program, however, we opt to enforce continuity weakly at faces between cells of different refinement level, for two reasons: (i) because we can, and more importantly (ii) because we would have to thread the automatic differentiation we use to compute the elements of the Newton matrix from the residual through the operations of the ConstraintMatrix class. This would be possible, but is not trivial, and so we choose this alternative approach.

I will take a look at the code from Praveen, and try it, but it looks as if it needs some work (some includes are missing or shifted). But thank you very much!

[Wolfgang Bangerth]

On 08/24/2017 09:11 AM, Praveen C wrote:
> I did this a while back, see here
>
> https://bitbucket.org/cpraveen/deal_ii/src/865f1c963f1089a870054bebf7180c9712e93da8/step-15-sacado/?at=master

Ah, most excellent!

I suspect that this can be written even more concisely with the current
development version of deal.II. For example, I'm pretty sure you can write

std::vector<Tensor<1,dim,>Sacado::Fad::DFad<double>>
grad_u (n_q_points);

fe_values.get_function_gradients (solution, grad_u);

or something similar in place of what you do in lines 199, 207-217. Do
you want to see if something like this works now? I think it would be
excellent to have an example like this in step-15.

@jppelteret -- as the AD master, do you have any other suggestion?

[Wolfgang Bangerth]

On 08/24/2017 09:22 AM, Maxi Miller wrote:
>
> In the first case, there is nothing we need to do: we are using a
> continuous finite element, and face terms do not appear in the
> bilinear form in this case. The second case usually does not lead to
> face terms either if we enforce hanging node constraints strongly
> (as in all previous tutorial programs so far whenever we used
> continuous finite elements – this enforcement is done by the
> ConstraintMatrix

> <https://www.dealii.org/8.4.0/doxygen/deal.II/classConstraintMatrix.html>
> class together with DoFTools::make_hanging_node_constraints
> <https://www.dealii.org/8.4.0/doxygen/deal.II/group__constraints.html#ga3eaa31a679484e80c193e74e8a967dc8>).

> In the current program, however, we opt to enforce continuity weakly
> at faces between cells of different refinement level, for two
> reasons: (i) because we can, and more importantly (ii) because we
> would have to thread the automatic differentiation we use to compute
> the elements of the Newton matrix from the residual through the
> operations of the ConstraintMatrix

> <https://www.dealii.org/8.4.0/doxygen/deal.II/classConstraintMatrix.html>

> class. This would be possible, but is not trivial, and so we choose
> this alternative approach.

Hm. I'm not sure I completely understand the issue any more. The second
part may or may not be correct. I think what it refers to is that not
all variables on a cell are actually degrees of freedom (e.g., if a
variable is defined at a hanging node) but may be tied to other
variables. Taking the derivative with regard to this variable may
therefore not immediately be what one had in mind. It's a different
question whether the result is actually wrong, or whether everything
works out when one distributes the local matrix into the global one with
ConstraintMatrix::distribute_local_to_global(). I think it's entirely
possible that the result turns out to be correct after all.

One would have to think about this some more, I'm afraid.

[Jean-Paul Pelteret]

I suspect that this can be written even more concisely with the current
development version of deal.II. For example, I'm pretty sure you can write

   std::vector<Tensor<1,dim,>Sacado::Fad::DFad<double>>
      grad_u (n_q_points);

   fe_values.get_function_gradients (solution, grad_u);

or something similar in place of what you do in lines 199, 207-217.

Yes, we now have the functions get_function_*_from_local_dof_values() in the FEValuesView classes, which can definitely help in this scenario. One would use them in the following way:

std::vector<double> local_dof_values(cell->get_fe().dofs_per_cell);
cell->get_dof_values(solution, local_dof_values.begin(), local_dof_values.end());
std::vector<Sacado::Fad::DFad<double>> local_dof_values_AD(cell->get_fe().dofs_per_cell);

// ... Setup local_dof_values_AD using local_dof_values, as is done in step-33 ...

std::vector<Tensor<1,dim,Sacado::Fad::DFad<double>>> qp_gradients_AD (n_q_points_cell);
fe_values[scalar_extractor].get_function_gradients_from_local_dof_values(local_dof_values_AD, qp_gradients_AD);
// ... use qp_gradients_AD ...

Regards,

Jean-Paul

[Maxi Miller]

How would I use the function get_function_gradients_from_local_dof_values with a vector_extractor? When trying to compile it with the vector shown above (including changing the first value) it fails with a rather long list of template errors.

Thanks!

[Daniel Arndt]

Maxi,

How would I use the function get_function_gradients_from_local_dof_values with a vector_extractor? When trying to compile it with the vector shown above (including changing the first value) it fails with a rather long list of template errors.

Thanks!

what are the template types you are using? How do the error messages look like?

You can find examples in which this function is used in

- tests/fe/fe_values_view_get_function_from_local_dof_values_01.cc [1]

- tests/fe/fe_values_view_get_function_from_local_dof_values_02.cc [2]

 

Best,

Daniel

[1] https://github.com/dealii/dealii/blob/master/tests/fe/fe_values_view_get_function_from_local_dof_values_01.cc

[2] https://github.com/dealii/dealii/blob/master/tests/fe/fe_values_view_get_function_from_local_dof_values_02.cc

[Maxi Miller]

I added the errors I get while compiling in the log file. My code for producing it is:

std::vector<Table<3, Sacado::Fad::DFad<double>>> grad_N_AD(n_q_points); std::vector<Sacado::Fad::DFad<double>> local_dof_values_AD(cell->get_fe().dofs_per_cell);
 

 for (unsigned int i=0; i<dofs_per_cell; ++i)
 {
 local_dof_values_AD[i] = present_solution(local_dof_indices[i]);
 local_dof_values_AD[i].diff(i, dofs_per_cell);
 }
 fe_values[surface_TE].get_function_gradients_from_local_dof_values(local_dof_values_AD, grad_N_AD);

Is that approach correct so far?

[Wolfgang Bangerth]

On 10/26/2017 02:18 AM, 'Maxi Miller' via deal.II User Group wrote:
> I added the errors I get while compiling in the log file. My code for
> producing it is:
>
> |

> std::vector<Table<3,Sacado::Fad::DFad<double>>>grad_N_AD(n_q_points);std::vector<Sacado::Fad::DFad<double>>local_dof_values_AD(cell->get_fe().dofs_per_cell);
>
>
> for(unsignedinti=0;i<dofs_per_cell;++i)

> {
> local_dof_values_AD[i]=present_solution(local_dof_indices[i]);
> local_dof_values_AD[i].diff(i,dofs_per_cell);
> }
> fe_values[surface_TE].get_function_gradients_from_local_dof_values(local_dof_values_AD,grad_N_AD);
> |
>
>
> Is that approach correct so far?

This looks reasonable. The code does not compile because of this error:

/opt/dealII/include/deal.II/base/template_constraints.h:388:40: error: no
match for 'operator*' (operand types are 'Sacado::Fad::DFad<double>' and
'dealii::Tensor<1, 2, double>')
typedef decltype(std::declval<T>() * std::declval<U>()) type;
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~

This should, in principle, be covered by the following operator in tensor.h:

template <int rank, int dim,
typename Number,
typename OtherNumber>
inline
Tensor<rank,dim,typename ProductType<typename EnableIfScalar<Number>::type,
OtherNumber>::type>
operator * (const Number &factor,
const Tensor<rank,dim,OtherNumber> &t)
{
// simply forward to the operator above
return t * factor;
}


with Number=DFad<double> and OtherNumber=double. You need to investigate
whether either the ProductType or EnableIfScalar is not defined for this pair
of operands. If you can figure this out, you should be in business.

[Wolfgang Bangerth]

On 10/26/2017 07:27 AM, Wolfgang Bangerth wrote:
>
> This looks reasonable. The code does not compile because of this error:
>
> /opt/dealII/include/deal.II/base/template_constraints.h:388:40: error:
> no match for 'operator*' (operand types are 'Sacado::Fad::DFad<double>'
> and 'dealii::Tensor<1, 2, double>')
>      typedef decltype(std::declval<T>() * std::declval<U>()) type;
>                       ~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~

By the way, did you try to include
#include <deal.II/differentiation/ad/sacado_product_types.h>
at the top of your code?

[Maxi Miller]

No, but including it removes the compilation errors, thanks! 

How do I know for the future which includes I should use?

[Wolfgang Bangerth]

On 10/26/2017 09:19 AM, 'Maxi Miller' via deal.II User Group wrote:
> No, but including it removes the compilation errors, thanks!
> How do I know for the future which includes I should use?

This one was difficult, but in general when you get a compiler error for
a function that the compiler can't find, you have to look up where that
function is declared. The manual pages online can help with that as they
list the file in question for each function.