Automatic Differentiation together with boundary_worker (step-72 + step-74)

[Nik Markus Leuenberger]

Dear all,

I'm opening a new thread because the question is no longer whether or how to use deal.ii for DG/finite volume implementations, but what I hope is a more specific question regarding a combination of step-74 and step-72.

I want to use automatic differentiation to obtain the Jacobian for a (non)-linear residual of a Poisson equation. (Non) is in brackets because I would like to solve the linear case first in a single Newton-Raphson step.

The question described in more detail in the attached .pdf is how to use automatic differentiation for the boundary_worker function of step-74.

Thank you very much in advance for any help.

Best regards,
Nik

[Wolfgang Bangerth]

On 7/8/24 08:51, Nik Markus Leuenberger wrote:
>
> I'm opening a new thread because the question is no longer whether or how to
> use deal.ii for DG/finite volume implementations, but what I hope is a more
> specific question regarding a combination of step-74 and step-72.
>
> I want to use automatic differentiation to obtain the Jacobian for a
> (non)-linear residual of a Poisson equation. (Non) is in brackets because I
> would like to solve the linear case first in a single Newton-Raphson step.
>
> The question described in more detail in the attached .pdf is how to use
> automatic differentiation for the boundary_worker function of step-74.

Nik:
I'm not sure anyone can help with this problem short of debugging it in detail
-- which I think you'll have to do yourself. But it may be useful to see how
others have done something similar, for example in this recent code gallery
program:

https://dealii.org/developer/doxygen/deal.II/code_gallery_nonlinear_heat_transfer_with_AD_NOX.html

Best
W.

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

[Jean-Paul Pelteret]

Hi Nik,

> My question is: for boundary_worker. What am I missing that the calculated residual contributions
don’t show up in the global RHS of the system?

I think that I might have an idea as to what's going on here. I feel that the problem likely lies with what these lines

Vector<double> &cell_rhs = copy_data.cell_rhs;
...
ad_helper.compute_residual(cell_rhs);

do in unison, versus the

copy_data.cell_rhs(i) +=

that you have in the hand written assembly code. Normally you'd be accumulating into the cell RHS/matrix (as you've demonstrated), but the compute_residual/linearization functions don't do this. Instead, they overwrite the local matrix/vector (although its a little deep into the AD framework, you can see this here https://github.com/dealii/dealii/blob/master/source/differentiation/ad/ad_drivers.cc#L1836-L1854). This is incompatible with the meshloop() framework as it reuses the same local matrices and vectors on a cell and its faces before scattering their contributions to the global linear system.

What you could do as a workaround is to simply create an intermediate local matrix and RHS vector in each of the workers, have the ADHelper write into those, and then accumulate the result into the data structures in copy_data.

I hope that my thoughts were on the right track, and that this helps. If it does, then this would motivate the addition of a set of functions along the lines of

ad_helper.accumulate_residual(copy_data.cell_rhs, -1.0 /* factor */);
ad_helper.accumulate_linearization(copy_data.cell_matrix);

that would work better with mesh_loop(). We'd be happy to accept any patches that help to realize this functionality!

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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/15834e9c-a00c-4f1a-a30e-602c2aaca296%40colostate.edu.

[Nik Markus Leuenberger]

Dear Jean-Paul,

On Wednesday, July 10, 2024 at 10:42:57 PM UTC+2 Jean-Paul Pelteret wrote:

Hi Nik,

> My question is: for boundary_worker. What am I missing that the calculated residual contributions
don’t show up in the global RHS of the system?

I think that I might have an idea as to what's going on here. I feel that the problem likely lies with what these linesVector<double> &cell_rhs = copy_data.cell_rhs;
... ad_helper.compute_residual(cell_rhs);

do in unison, versus the

copy_data.cell_rhs(i) +=

that you have in the hand written assembly code. Normally you'd be accumulating into the cell RHS/matrix (as you've demonstrated), but the compute_residual/linearization functions don't do this. Instead, they overwrite the local matrix/vector (although its a little deep into the AD framework, you can see this here https://github.com/dealii/dealii/blob/master/source/differentiation/ad/ad_drivers.cc#L1836-L1854). This is incompatible with the meshloop() framework as it reuses the same local matrices and vectors on a cell and its faces before scattering their contributions to the global linear system.

This is exactly what was happening!

What you could do as a workaround is to simply create an intermediate local matrix and RHS vector in each of the workers, have the ADHelper write into those, and then accumulate the result into the data structures in copy_data.

Yes, as a workaround, I did create an intermediate local matrix called local_cell_matrix and RHS vector called local_rhs that then does something like this. This seems to fix the problem. 

Thank you very much.

// compute residual and linearization on the local matrix and vector

            ad_helper.register_residual_vector(residual_ad);

            ad_helper.compute_residual(local_rhs);

            local_rhs *= -1.0;

            ad_helper.compute_linearization(local_cell_matrix);

            // add contributions from local matrix and vector to copy object

            for (int i = 0; i < fe_iv.dof_indices().size(); i++)

            {

                for (int j = 0; j < fe_iv.dof_indices().size(); j++)

                {

                    copy_data_face.cell_matrix(i, j) += local_cell_matrix(i, j);

                }

                copy_data_face.cell_rhs(i) += local_rhs(i);

            }

I hope that my thoughts were on the right track, and that this helps. If it does, then this would motivate the addition of a set of functions along the lines of ad_helper.accumulate_residual(copy_data.cell_rhs, -1.0 /* factor */); ad_helper.accumulate_linearization(copy_data.cell_matrix);

that would work better with mesh_loop(). We'd be happy to accept any patches that help to realize this functionality!

Thank you also for pointing me to the location of the source code where the = instead of += is located.

Your two suggested functions would indeed make the ad_helper work much more conveniently with the mesh_loop framework (especially the face_worker and boundary_worker).

I will try to see how one could implement these functions and in case I'm successful, I will follow the process described here (https://github.com/dealii/dealii/wiki/Contributing) to contribute the patch.

Thank you very much once again and best regards,

Nik