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/278e2f76-155a-4de2-b1bf-970585c822a0%40googlegroups.com.
