Loss of symmetry by changing parameter

[Sylvain Mathonnière]

Dear all,

I am experiencing a symmetry issue in my calculation with certain set of parameters.
In the 3 figures below, I used different sets of parameters (same equation) and plotted the result once the simulation reaches steady state.
The two first figures look fairly symmetric (as it should be) but the 3rd one is obviously not symmetric.

Upon closer inspection, one can see that the bottom left corner looks a bit "shaky" on the first 2 figures.

I tried to get to the bottom of why this happens numerically ?
I found that at some point during my assembly of the RHS, there is a difference in what I obtained even though it should be perfectly symmetric. The solver takes also significantly longer to solve it.

I represented the RHS on a coarse grid to highlight the difference.

By digging even more I realised that the discrepancy steams from the boundary_worker (similar to the one of step-12).
On the fig. 4, the bottom right cell has its form function different than 0 for dof =0 (angle dof) (i.e fe_face.shape_value(dof = 0; point_gauss) =/= 0 )

whereas on the fig. 5, the bottom left cell has fe_face.shape_value(dof = 0; point_gauss) = 0 for all gauss point.

For the record, I am solving several time the radiative transport equation I  using the discontinuous galerkin method using several directions.

So the fig.4 is the RHS for one direction beta=(0.14; 0.14) and fig 5. is the RHS for the other direction beta=(0.14;-0.14). The BC

My question is the following:

Is this a normal behaviour that I am experiencing and the lack of symmetry in fig. 3 is linked with the fact that my mesh is not fine enough ?

or could it be the result of something wrong in my code (I attached a snipset of it below) ?

I am not experienced enough to know the difference. I know the code works well because I compared it to scientific reference + using the method of manufactured solution, so if there is a mistake, it has to be "small".

Best regards,

Sylvain Mathonnière

[Sylvain Mathonnière]

Always some errors are left after I wrote the message... errare humanum est.

I inversed the figure in the text, I meant :

"On the fig. 5, the bottom right cell has its form function different than 0 for dof =0 (angle dof) (i.e fe_face.shape_value(dof = 0; point_gauss) =/= 0 )

whereas on the fig. 4, the bottom left cell has fe_face.shape_value(dof = 0; point_gauss) = 0 for all gauss point.

Also when I wrote "The BC", I wanted to say that the boundary conditions are the same in both cases and are actually implemented directly in the assembly (like in step-12).

[Wolfgang Bangerth]

> Upon closer inspection, one can see that the bottom left corner looks a bit
> "shaky" on the first 2 figures.

This is an artifact of plotting. What you are trying to visualize is a
bilinear function on the bottom left cell, but Paraview/Visit actually show is
a subdivision of the square cell into two triangles, and then using a linear
interpolation on each triangle. The way the square is subdivided into
triangles uses the bottom right -> top left diagonal. As a consequence, the
bottom left and the bottom right cell are not shown the same way, even their
values are perfectly symmetric.

> I tried to get to the bottom of why this happens numerically ?
> I found that at some point during my assembly of the RHS, there is a
> difference in what I obtained even though it should be perfectly symmetric.

Are you observing this visually, or are you looking at actual numbers?

> By digging even more I realised that the discrepancy steams from the
> boundary_worker (similar to the one of step-12).
> On the fig. 4, the bottom right cell has its form function different than 0

> for dof =0 (angle dof) (i.e *fe_face.shape_value(dof = 0; point_gauss) =/= 0* )
> whereas on the fig. 5, the bottom left cell has *fe_face.shape_value(dof = 0;
> point_gauss) *= 0 for all gauss point.

It is difficult to imagine that something this basic doesn't work given how
widely these functions are used. But since you think you know what is
happening, can you strip down your program to as-simple-as-possible a program
that illustrates the issue?

> So the fig.4 is the RHS for one direction beta=(0.14; 0.14) and fig 5. is the
> RHS for the other direction beta=(0.14;-0.14). The BC

This too does not actually look wrong to me, if you consider how Paraview
plots the rhs vector. You can only trust the values in the vertices of that
bottom row of cells, and the two pictures look perfectly symmetric to me.

Best
W.

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

[Sylvain Mathonnière]

Dear Wolfgang Bangerth,

I am looking at actual numbers because I figured out I could not trust completely paraview images.

I managed to make a minimun working example. It is still 400 lines of code though.
Almost most of it is very standard deal.II "routine" that I borrowed from the tutorials.

I do the following

-generate the mesh,
-setup the system
-initialize my variable "old_T" using VectorTools::project. It looks like this :

Then for each direction I do :
 - assemble my RHS and I export it immediately.


Since my problem occurs at the RHS definition, I did not assemble the matrix and solve to keep it as simple as possible.

Since I have two directions, the program output 2 vtk files that I can plot using paraview and I get this :

I also recorded the values of the corner dof in the code :

The points 0; 21; 42 and 63 correspond to the corner dof of my mesh as can be seen here :



One can see that the corner 21 and 63 are indeed symetric when I change direction but it is not the case for corner 0 and 42. Since old_T is symetric it should be the same.

The main issue in the code is the assembly, more precisely on line 305 : copy_data.cell_rhs_RTE(i) +=  pow(T_old_data[point],4) * 1e-8 * fe_face.shape_value(i, point)*JxW[point];
if I remove the pow(T_old_data[point],4), then it gets perfectly symmetric but otherwise no. I cannot get my mind around the reason why.

Also, later when I solve the system, the solver takes way longer for the 2nd direction whereas it should be perfectly symetric.
Any help would be greatly appreciated.

Best,

Sylvain Mathonnière

[Wolfgang Bangerth]

On 10/7/21 8:36 AM, Sylvain Mathonnière wrote:
>
> The main issue in the code is the assembly, more precisely on line 305 :
> copy_data.cell_rhs_RTE(i) +=  pow(T_old_data[point],4) * 1e-8 *
> fe_face.shape_value(i, point)*JxW[point];
> if I remove the pow(T_old_data[point],4), then it gets perfectly
> symmetric but otherwise no. I cannot get my mind around the reason why.

Let's do some more remote debugging then :-) Since you've already
identified that the problem appears to be the old values T_old_data,
output those as well. You should get ... actually, I found your bug: You
are looping over the quadrature points of the FEFaceValues object here:

const auto &q_points = fe_face.get_quadrature_points();
...
for (unsigned int point = 0; point < q_points.size(); ++point)
{
double beta_dot_n = direction_temp * normals[point];

if (beta_dot_n < 0)
{
for (unsigned int i = 0; i < n_facet_dofs; ++i)
{
// problem is here

copy_data.cell_rhs_RTE(i) += pow(T_old_data[point],4) * 1e-8
* fe_face.shape_value(i, point)*JxW[point];

}

So one would expect that T_old_data[point] is also indexed by the
quadrature points on the face. But you initialize

fe_values_HE.reinit(cell_2);
fe_values_HE.get_function_values(old_T,T_old_data);

that is, T_old_data is indexed by the quadrature points on the *cell*,
not the face. This can't be right.

Separately, you use the same T_old_data object for all face_workers. But
these may be called in parallel, and you will get race conditions. You
need to make sure you have a separate T_old_data object on every thread,
for example by making the variable local to the face_worker lambda function.

[Sylvain Mathonnière]

Oh thank you very much ! That was indeed the problem, I changed the code by calling FEValues using the face quadrature formula and used it locally to avoid those race conditions. It looks like this :

    // new cell iterator for other dof_handler

    typename DoFHandler<dim>::active_cell_iterator cell_2 (&cell->get_triangulation(), cell->level(), cell->index(), &dof_handler_HE);

    // Initialise FEValues with FACE quadrature formula and reinitialise with current cell

    FEValues<dim> fe_values_HE(mappingQ, fe_HE, q_points, update_values);

    fe_values_HE.reinit(cell_2);

   

    // get data at quadrature point ON THE FACE

    std::vector<double> T_old_data(q_points.size());

    fe_values_HE.get_function_values(old_T,T_old_data);

You managed to solve 2 of the problem I was having. I believe the code works as intended now.

Those kind of errors are quite tricky to catch. It did not show much in the results I was getting and I am pretty sure the solution I had (except the set parameter 3 which was obviously wrong) would have converged by decreasing the size of the mesh.

Moreover the method of manufactured solution could not catch this error at all since it was in the RHS.

How can one suspect there is such problem with the code in those cases ? Just with experience I guess and doing multiple simulations ?

Regardless, thank you very much again ! I hope one day I will be able to help back this community as much as it helped me so far.

Best,

Sylvain Mathonnière

[Wolfgang Bangerth]

On 10/8/21 5:13 AM, Sylvain Mathonnière wrote:
> How can one suspect there is such problem with the code in those cases ? Just
> with experience I guess and doing multiple simulations ?

Making every mistake myself at least once in the last 20 years has certainly
helped me have an eye for things :-)

[Sylvain Mathonnière]

For the record, the lines of code in my last post were wrong. It should be for the boundary worker :

    // new cell iterator for other dof_handler

    typename DoFHandler<dim>::active_cell_iterator cell_2 (&cell->get_triangulation(), cell->level(), cell->index(), &dof_handler_HE);

    // Initialise FEValues with FACE quadrature formula and reinitialise with current cell

    FEFaceValues<dim> fe_values_HE_bound(mappingQ, fe_HE, quadrature_face_formula, update_values);

    fe_values_HE_bound.reinit(cell_2, face_no);

    

    std::vector<double> T_old_data_bound(q_points.size());

    fe_values_HE_bound.get_function_values(old_T,T_old_data_bound);

Best,

Sylvain M.