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/2ded867b-a2e7-4041-b90c-2dd774f1f5dd%40googlegroups.com.