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


Best,
Jean-Paul


On Tue, 9 Jul 2024 at 05:17, Wolfgang Bangerth <[email protected]> 
wrote:

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

-- 
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 [email protected].
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/fcf9ad8a-952d-4272-af18-e26715ba07d6n%40googlegroups.com.

Reply via email to