Dear Felix, I don’t claim to have read through your comprehensive documentation of your issue in any detail, but one observation of your code is that you have a couple of suspect terms in your local matrix assembly. In particular, the terms of this form > 1/2/phiscale*(u_j_y-v_j_x) do not compute to what you think they do. The compiler will read this from left to right, and the first operation it sees is an integer divided by an integer. Integer division in C++ will return the value of the quotient rounded down (towards zero, if positive). There’s a more detailed explanation here: https://stackoverflow.com/a/3602857 <https://stackoverflow.com/a/3602857> So these terms always compute zero! Here’s a short example if you need further proof (looks at the debug output in the central column): https://godbolt.org/z/XJZ59h <https://godbolt.org/z/XJZ59h>
So what you need to do to fix this is to use operator precedence to make sure that you’re always performing division with floating point numbers, or (IMO the preferable solution) explicitly specify the factor 1/2 as a floating point value, either by writing it as 1.0/2.0 or just 0.5. I’ve got no idea if this is the only source of error in your code, but its certainly not helping you :-) Good luck in sorting this out! Best, Jean-Paul > On 25 Jul 2019, at 16:10, Félix Bunel <[email protected]> wrote: > > One final remark for today, if I remove the two asymmetrical terms in my > equation so it becomes : > > <Auto Generated Inline Image 1.png> > > Then all three codes gives back the same solution... > > So I thought it would come from my assembly of this part of the system but I > have checked and checked again... It should not come from my weak formulation > since the fenics and dealii have the same one. > Something to do with the matrix no longer being symmetrical ? But I used > UMFPack as a solver so it should not be an issue... > > Here is my assembly process for the stokes problem in case you see something > wrong > void Elfe::stokes_assemble_system () > { //Fonction qui assemble le problème de Stokes > > std::cout << BLUE << "Assemblage Stokes " << RESET > << std::endl; > > //On réinitialise la matrice et le rhs > stokes_system_matrix = 0; > stokes_system_rhs = 0; > > //On charge de quoi extraire les valeurs > QGauss<2> quadrature_formula(4); > const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; > const unsigned int n_q_points = quadrature_formula.size(); > > std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell); > > FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell); > Vector<double> local_rhs (dofs_per_cell); > > //Stokes values > FEValues<2> stokes_fe_values(stokes_fe, quadrature_formula, > update_values | update_gradients | > update_JxW_values); > > double u_i, v_i; > double u_i_x, u_i_y, v_i_x, v_i_y; > double u_j, v_j; > double u_j_x, u_j_y, v_j_x, v_j_y; > double P_i, P_j; > > const FEValuesExtractors::Vector Velocities(0); > const FEValuesExtractors::Scalar Pressure(2); > > //Phi values > FEValues<2> phi_fe_values(phi_fe, quadrature_formula, > update_values | update_gradients | > update_JxW_values); > > double phi, phi_x, phi_y; > double DX; > > double philin = 2*pi*n; > double phiscale = 2*pi*s; > > std::vector<Tensor<1, 2> > phi_solution_gradients(n_q_points); > std::vector<double> phi_solution_values(n_q_points); > > //Il nous faut 2 cell iterator, un pour phi, un pour stokes > typename DoFHandler<2>::active_cell_iterator > stokes_cell = stokes_dof_handler.begin_active(), > stokes_endc = stokes_dof_handler.end(); > > typename DoFHandler<2>::active_cell_iterator > phi_cell = phi_dof_handler.begin_active(); > > for (; stokes_cell!=stokes_endc; ++stokes_cell, ++phi_cell) > { //Pour chaque cell > stokes_fe_values.reinit(stokes_cell); > phi_fe_values.reinit (phi_cell); > > local_matrix = 0; > local_rhs = 0; > > //On récupère les valeurs pour les angles phis > phi_fe_values.get_function_values(phi_solution, > phi_solution_values); > phi_fe_values.get_function_gradients(phi_solution, > phi_solution_gradients); > > for (unsigned int q=0; q<n_q_points; ++q) > { //Pour chaque point de quadrature > DX = stokes_fe_values.JxW(q); > > phi = phi_solution_values[q]; > phi_x = phi_solution_gradients[q][0]; > phi_y = phi_solution_gradients[q][1]; > > for (unsigned int i=0; i<dofs_per_cell; ++i) > { > u_i = stokes_fe_values[Velocities].value (i, q)[0]; > v_i = stokes_fe_values[Velocities].value (i, q)[1]; > u_i_x = stokes_fe_values[Velocities].gradient(i, q)[0][0]; > u_i_y = stokes_fe_values[Velocities].gradient(i, q)[0][1]; > v_i_x = stokes_fe_values[Velocities].gradient(i, q)[1][0]; > v_i_y = stokes_fe_values[Velocities].gradient(i, q)[1][1]; > P_i = stokes_fe_values[Pressure] .value (i, q); > > local_rhs(i) += ( > > u_i*( > 8*philin*phiscale*X > *(cos(2*phiscale*phi)*phi_x + > sin(2*phiscale*phi)*phi_y) > > + 2*phiscale*phiscale*phi_x*(-4*philin/phiscale) > ) > > + v_i*( > 8*philin*phiscale*X > *(sin(2*phiscale*phi)*phi_x - > cos(2*phiscale*phi)*phi_y) > > + 2*phiscale*phiscale*phi_y*(-4*philin/phiscale) > ) > > )*DX; > > for (unsigned int j=0; j<dofs_per_cell; ++j) > { > u_j = stokes_fe_values[Velocities].value (j, q)[0]; > v_j = stokes_fe_values[Velocities].value (j, q)[1]; > u_j_x = stokes_fe_values[Velocities].gradient(j, q)[0][0]; > u_j_y = stokes_fe_values[Velocities].gradient(j, q)[0][1]; > v_j_x = stokes_fe_values[Velocities].gradient(j, q)[1][0]; > v_j_y = stokes_fe_values[Velocities].gradient(j, q)[1][1]; > P_j = stokes_fe_values[Pressure] .value (j, q); > > local_matrix(i, j) += ( > //Classical Stokes terms > - u_i_x*u_j_x - u_i_y*u_j_y - v_i_x*v_j_x - > v_i_y*v_j_y > + (u_i_x+v_i_y)*P_j > + P_i*(u_j_x+v_j_y) > > //Leslie Asymetrical terms > + phiscale/a > *(v_i_x-u_i_y) > *(u_j*phi_x + v_j*phi_y + > 1/2/phiscale*(u_j_y-v_j_x)) > > //Terme provenant du laplacien phi > -2*phiscale*phiscale/a > *(u_i*phi_x + v_i*phi_y) > *(u_j*phi_x + v_j*phi_y + > 1/2/phiscale*(u_j_y-v_j_x)) > > )*DX; > } > } > } > > stokes_cell->get_dof_indices (local_dof_indices); > > //On applique les boundary conditions > stokes_constraints.distribute_local_to_global(local_matrix, > local_rhs, > local_dof_indices, > stokes_system_matrix, > stokes_system_rhs); > } > > } > > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/a984ed95-fa81-4758-bc46-1357ec5ce2b5%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/a984ed95-fa81-4758-bc46-1357ec5ce2b5%40googlegroups.com?utm_medium=email&utm_source=footer>. > <Auto Generated Inline Image 1.png> -- 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/03D394E6-AF2F-43AA-B2F5-C73E61D05193%40gmail.com.
