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
El jueves, 7 de octubre de 2021 a la(s) 23:04:15 UTC+2, Wolfgang Bangerth
escribió:
> 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.
>
> 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/b526b942-7ef7-42ca-b6d5-d8ae25849d94n%40googlegroups.com.
