Dear Richard, You are absolutely right in saying that if you have the Stokes equation your problem should be linear, even with the PSPG stabilization term. Laplacian of u is only going to be zero for Q1-Q1 elements if your mesh is perfectly aligned with the axis of your domain. However, this is what you mentioned so I do not think your error stems from there. If it I were you, I would do as Wolfgang suggested and try to look at what happens when you use a direct solver to solve the system (say UMFPACK). Which procedure are you using the define the "pressure constant" if you domain is enclosed?
Best Bruno On Wednesday, 19 February 2020 03:11:20 UTC-5, Richard Schussnig wrote: > > Hi Bruno, > Thank you very much for taking the time to read my post & giving me some > ideas! > I am already familiar with lethe, thats where i compared my implementation > of the residual to, but for Stokes, > we actually do not introduce a nonlinearity, since we do not have the > (nonlinear) convective term as in Navier Stokes, so it should actuallly be > fine, correct? > As far as I know, the Laplacian(u) in the residual does not lead to > nonlinearities, since the residual in PSPG for Stokes(!) is only r := > -nu*laplace(u) + grad(p) - f > and we simply add on element level + int( r * grad(q) )dx > Nonetheless, I solve it using Newton OR Picard+Aitken, giving me identical > results. > > That it takes longer to reach the asymptotic convergence I can confirm > (using my matlab code), > but using deal.ii and an iterative solver it just does not reach the point > (I tested up to 3 mio cells in 2d on just a rectangular grid with > poisuille flow ... didnt help) > thanks also for that hint! > > Kind regards, > Richard > > Am Mittwoch, 19. Februar 2020 02:49:26 UTC+1 schrieb Bruno Blais: >> >> Dear Richard, >> We implement a quasi identical scheme in lethe (PSPG for continuity and >> SUPG for momentum). You can find the code on the github, it might help you : >> >> https://github.com/lethe-cfd/lethe/blob/master/source/solvers/gls_navier_stokes.cc >> >> I can confirm that you should be getting order p+1 for velocity. To be >> honest I have never really taken the time to look at pressure convergence, >> but I recall that it can take very fine mesh to reach the asymptotic >> convergence. >> >> Please keep in mind that the introduction of the PSPG term introduces a >> non-linearity in the solver except in the case where you rmesh is made of >> perfect rectangles that are aligned with the axis of the system. >> Are you solving it using Newton's method? >> >> Best >> Bruno >> >> >> On Friday, 14 February 2020 07:34:09 UTC-5, Richard Schussnig wrote: >>> >>> Hi everyone! >>> I am currently implementing some stationary Stokes solvers based on >>> step-55. >>> Therein, Taylor-Hood elements are being used. One can check the optimal >>> order of >>> convergence easily comparing to the Kovasznay or Poisuille flow >>> solution, the first one being already implemented in this step. >>> I went ahead and added PSPG and Bochev-Dohrmann stabilizations to check >>> the errors using those. >>> >>> For the latter, i get the expected suboptimal (!) convergence rates >>> (vel:2 and p~1.5-2). >>> Using PSPG, one adds >>> tau * residual momentum dot gradient(q) >>> (q being the pressure test function) >>> per element like: >>> >>> local_matrix (i,j) += fe_values.JxW(q) * tau_sups * phi_grads_p [j] * >>> phi_grads_p [i]; >>> >>> if (ADD_THE_LAPLACIAN) >>> { >>> local_matrix (i,j) += fe_values.JxW(q) * tau * ( >>> - viscosity * phi_laplacians_v [j] >>> ) * phi_grads_p [i]; >>> } >>> >>> Using the laplacian of the velocity u defined as below: >>> >>> if (get_laplace_residual) >>> { >>> phi_hessians_v[k] = fe_hessians[velocities].hessian(k,q); >>> for (int d = 0; d < dim; ++d) >>> { phi_laplacians_v[k][d] = trace(phi_hessians_v[k][d]); >>> } >>> >>> // Check laplacian, both give zero for rectangular grid and identical >>> values for distorted grids. >>> unsigned int comp_k = fe.system_to_component_index(k).first; >>> >>> Tensor<1,dim> vector_laplace_v; >>> vector_laplace_v = 0; >>> if (comp_k<dim) >>> { >>> Tensor<2,dim> nonzero_hessian_k = >>> fe_hessians.shape_hessian_component(k,q,comp_k); >>> vector_laplace_v [comp_k] += trace(nonzero_hessian_k); >>> } >>> Tensor<1,dim> diff; >>> diff = phi_laplacians_v[k] - vector_laplace_v; >>> >>> if (diff.norm() > 1e-6 || vector_laplace_v.norm() > 1e-16) >>> std::cout << "|divgrad_v| = " << std::setw(15) << >>> vector_laplace_v.norm() << " , diff = " << diff.norm() << std::endl; >>> >>> } >>> >>> Which also includes a second option to compute the wanted laplacian. >>> >>> Using a perfectly rectangular grid, the laplacian is actually zero for >>> all elements, thus everything reduces to just adding a pressure stiffness >>> in the pressure-pressure block. >>> Nonetheless, i observe convergence rates of vel:2 and p:1.5-1.8 which is >>> suboptimal and not(!) expected. >>> I implemented the same methods in a matlab inhouse code & observed >>> perfect convergence rates of 2&2 for v&p using a direct solver. >>> >>> The stabilization parameter is in both cases (matlab or deal.ii) chosen >>> as >>> tau = 0.1*h^2/viscosity >>> with >>> double h = std::pow(cell->measure(), 1.0/ static_cast<double>(dim)) >>> or h = cell->diameter() giving similar results. >>> >>> Could the observed behaviour be caused by applying an iterative solver*? >>> (using ReductionControl and a reduction factor of 1e-13 which is really low) >>> * Block-triangular Schur-complement-based approach, similar as in >>> step-55 suggested in the further possibilities for extension, basically >>> using S~-(1/viscosity)*Mass_pressure giving 20-40 iterations up to 3mio >>> dofs tested. >>> >>> I have been checking the code for a week now, and i really doubt, that >>> the integration of just tau*grad(p)*grad(p) is wrong (again, the laplace >>> drops out for a rectangular grid). >>> Just to check, i used a manufactured solution from Poisuille flow and >>> added the laplacian from PSPG to the rhs, giving me the exact (linear) >>> solution in the pressure and quadratic convergence in the velocity, thus I >>> assume, that the grad(p)grad(q) term added is correctly implemeted. (exact >>> solution meaning, that i get an L2 error of err<1e-9 for any number of >>> elements used) >>> >>> Anyone has ever experienced this or has anyone some tipps for further >>> debugging? >>> Thanks for taking the time to read the post, any help or comment would >>> be greatly appreciated! >>> >>> Greetings from Austria, >>> Richard >>> >> -- 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 dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7a7a6380-37cd-4d72-b9db-81ca67663f0f%40googlegroups.com.
