Hi Felix, > Okay. When I read you e-mail I sincerely hoped it wasn't the problem. But it > was. Now the results match.
That’s great that this was an easy fix! I’m glad that you can now move along with your problem :-) > I feel so so stupid. Such a rookie mistake. I already made this mistake 3 > years ago and forgot about it. Well now, I bet each one of us has who has been in the C++ game for a while carries around a list of “obvious” mistakes that they first look for when evaluating their own or someone else’s code. Don’t be too upset about this. I’d be embarrassed to admit how many times I’ve been tripped up by (what turn out to be) trivial errors. The important thing is to learn from it and try your hardest not to do it again! > Thank you very much for taking the time to look through this part of the > code. I own you big time. > > Thank you also Daniel for helping me. I will still try a manufactured > solution to learn about this process. You’re very welcome. I only had a couple of minutes to help, so its entirely coincidental that my contribution happened to be the one that assisted you the most. I pass the owed favour over to Daniel, who spent a lot more brain power thinking about your problem than I did. And I think that there’s certainly something to be learned about implementing a manufactured solution, so go for it! > Thanks also to everyone that contributed to dealii. Its always a delight for us to hear that the users of the library are happy with it. Thank you for your nice compliment! We hope that in the future you’ll also be on the contributors list, with either a little patch or hopefully something more! > Its the perfect library for people that are not as stupid as me. Seriously, don’t be too hard on yourself. We’ve all made (and will continue to make) these sorts of mistakes from time to time :-) Best, Jean-Paul > On Thursday, July 25, 2019 at 7:28:13 PM UTC+2, Jean-Paul Pelteret wrote: > 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 <bunel...@ <>gmail.com >> <http://gmail.com/>> 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 dea...@ <>googlegroups.com <http://googlegroups.com/>. >> 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/ > <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/f887ca46-262f-4805-8b94-ef68341cb455%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/f887ca46-262f-4805-8b94-ef68341cb455%40googlegroups.com?utm_medium=email&utm_source=footer>. -- 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/88C52C56-D027-40E1-8B7A-0DA37F38911A%40gmail.com.
