Sacado automatic differentiation

[heikki]

Hi, I have try to solve a nonlinear problem which has a very unpleasant Jacobian matrix. So, I would like to use automatic differentiation. There is an example how to use Sacado, but I am not still absolutely sure how to use this. (step-33) Let's assume that I have a nonlinear Poisson problem, for example

div( grad phi) - exp(phi) phi = 0

with Dirichlet and Neumann boundary conditions. I try to solve this using Newton's method

F'(u_k,delta u) = - F(u_k)
u_k+1 = u_k+delta u

When I am assembling the system matrix and right hand side, I am doing the following:

std::vector<types::global_dof_index> dof_indices;
std::vector<Sacado::Fad::DFad<double> > independent_local_dof_values(dofs_per_cell);
std::vector<double> derivatives (dofs_per_cell);

for(; cell!=endc; ++cell){

         fe_values.reinit (cell);

         cell->get_dof_indices (dof_indices);

         for (unsigned int i=0; i<dofs_per_cell; ++i){

                independent_local_dof_values[i] = old_solution(dof_indices[i]);
         for (unsigned int i=0; i<dofs_per_cell; ++i){
                independent_local_dof_values[i].diff (i, dofs_per_cell);
      

         for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i){

                Sacado::Fad::DFad<double> F_i = 0;

                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point){

                           for (unsigned int d=0; d<dim; d++){

                                        F_i -= old_solution(dof_indices[i])*fe_v.shape_grad_component(i, point, component_i)[d] *
                                                  independent_local_dof_values[i]*fe_v.shape_value_component(i, point, component_i) -
                                                  std::exp(old_solution(dof_indices[i])*fe_v.shape_grad_component(i, point, component_i)[d]) *
                                                  old_solution(dof_indices[i])*fe_v.shape_value_component(i, point, component_i) *
                                                  independent_local_dof_values[i]*fe_v.shape_value_component(i, point, component_i) *
                                                  fe_v.JxW(point);

                            }
                 }       
                 for (unsigned int k=0; k<dofs_per_cell; ++k)

                            derivatives[k] = F_i.fastAccessDx(k);

                 system_matrix.add(dof_indices[i], dof_indices, derivatives);

                 right_hand_side(dof_indices[i]) -= F_i.val();

          }
}

I feel a little bit insecure how to do this, but this pretty much what I should do?
old_solution is a vector which contains dof's from the previous iteration.

-Heikki

 

[Wolfgang Bangerth]

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#aed056da8b247808a5ea995e2743d7083>;

> ++i){
> Sacado::Fad::DFad<double> F_i = 0;
> for (unsigned int point=0;
> point<fe_values.n_quadrature_points

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a8883dd58663474d7190e3a1cdb6070a0>;

> ++point){
> for (unsigned int d=0; d<dim; d++){
> F_i -=
> old_solution(dof_indices[i])*fe_v.shape_grad_component

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a65ad58f71e65b030018976eba21904b1>(i,

> point, component_i)[d] *
> independent_local_dof_values[i]*fe_v.shape_value_component

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a28c58c7ea9dca6da66e8ff528670232e>(i,

> point, component_i) -
>
> std::exp(old_solution(dof_indices[i])*fe_v.shape_grad_component

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a65ad58f71e65b030018976eba21904b1>(i,

> point, component_i)[d]) *
>
> old_solution(dof_indices[i])*fe_v.shape_value_component

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a28c58c7ea9dca6da66e8ff528670232e>(i,

> point, component_i) *
> independent_local_dof_values[i]*fe_v.shape_value_component

> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#a28c58c7ea9dca6da66e8ff528670232e>(i,
> point, component_i) *
> fe_v.JxW
> <http://www.dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#ab9832bac82d6a78d1592861faa59c665>(point);

> }
> }
> for (unsigned int k=0; k<dofs_per_cell; ++k)
> derivatives[k] = F_i.fastAccessDx(k);
> system_matrix.add(dof_indices[i], dof_indices,
> derivatives);
> right_hand_side(dof_indices[i]) -= F_i.val();
> }
> }
>
> I feel a little bit insecure how to do this, but this pretty much what I
> should do?

This looks basically correct to me, though that's of course not a
guarantee. Does it work?

Best
W.


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

[heikki]

No, it doesn't work. I have tried to modify step-15 such that its Jacobian can be solved also with Sacado. Unfortunately, something is wrong, because my right hand side is not matching with the original vector. (Then, of course, system matrix is also wrong too) I suspect that F_i is not correct.
Actually, independent variables are not used correctly in my previous post and I would guess that I am not using these correctly.

-Heikki

[Wolfgang Bangerth]

I don't think anyone on the mailing list has sufficient knowledge of
this issue any more. step-33 was written many years ago and the author
(to the best of my knowledge) is no longer on this mailing list.

I'm afraid you'll have to play around with simple cases yourself to
figure out where it is going wrong. I would try to use a simple
nonlinearity and compare what you're getting with what it should be.

[Praveen C]

Hi

I have played with sacado in the past and I changed step-15 to use sacado, see

https://code.google.com/p/cfdlab/source/browse/trunk/deal.ii/step-15-sacado/

It seems to be working but I have not thoroughly checked the results. The only changes are in assemble_system function.

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

[heikki]

Excellent! Thanks a lot! Your example was very useful Praveen..