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/fad25991-1271-cc8a-7752-554ed5aea3d8%40colostate.edu.
