Dear Richard,
Thanks for your message! It is very interesting. Your results made me doubt
our own results, so I re-ran Error convergence analysis on a manufactured
solution ( an infinitely continuous one) where the domain is bounded by
purely Dirichlet boundary condition.
I did the simulations with two approaches:
Q1-Q1 using the GLS (SUPG +PSPG) stabilized solver of Lethe using a
monolithic iterative solver and Q2-Q1 using the Grad-Div solver with a
Schur completement solution strategy (similar to step-57 but using
Trilinos).
*Here are the results I obtain:*
*Q1-Q1 - Velocity error and convergence*
cells error
64 1.3282e-01 -
256 3.4363e-02 1.95
1024 8.7362e-03 1.98
4096 2.1969e-03 1.99
16384 5.5029e-04 2.00
65536 1.3767e-04 2.00
262144 3.4426e-05 2.00
1048576 8.6075e-06 2.00
*Conclusion : Order = p+1 recovered for the velocity!*
*Q1-Q1 - Pressure error and convergence*
Refinement level L2Error-Pressure Ratio
1 0.171975
2 0.0964683 1.33
3 0.0301739 1.78
4 0.00844618 1.89
5 0.0023758 1.885
6 0.000698421 1.84
7 0.000216927 1.79
8 7.07505e-05 1.75
*Conclusion : Order = ??? The error does not reach the right asymptotic
convergence rate. This is with a very low non-linear solver residual.
Clearly, There is a significant loss of accuracy in the pressure. The error
on the pressure is also orders of magnitudes higher than the grad-div
solver...*
*Q2-Q1 - Velocity error and convergence*
cells error
64 1.4729e-02 -
256 1.9368e-03 2.93
1024 2.4521e-04 2.98
4096 3.0749e-05 3.00
16384 3.8466e-06 3.00
65536 4.8237e-07 3.00
*Conclusion : Order = p+1 recovered for the velocity!*
*Q2-Q1 - pressure error and convergence*
Refinement level L2Error-Pressure Ratio
1 0.0306665
2 0.00399144 2.77183454825005
3 0.00084095 2.1786111158165
4 0.000206301 2.01899117611792
5 5.1481e-05 2.00182993535217
6 1.28777e-05 1.99942139684848
*Conclusion : Order = p+1 recovered for the velocity!*
If you want, I would be more than glad to discuss this issue even more. We
could organize maybe a skype call to discuss this? Clearly we are having
very very similar issues :)!
Looking forward to more of this discussion.
Best
Bruno
On Tuesday, 10 March 2020 12:05:26 UTC-4, Richard Schussnig wrote:
>
> Hi everyone,
> I am back with good & bad news!
> - Good news first, I implemented a parallel direct solver (mumps via
> petsc) and going from a 2x2 blocked system to
> a regular one (actually 1x1 still blocked system)
> one just needs to adapt the setup phase, not reordering per component &
> not use the "get_view"
> but rather use all the same vector and matrix (types) with 1 block instead
> of 2. Once you realize that, it takes ~15 mins for a
> non-experienced deal.ii user like me.
>
> - Bad news is, that I still do not obtain the optimal convergence rate in
> the pressure, i.e., 2 for Q1Q1 PSPG-stabilized Stokes.
>
> The obtained rates in the Poisuille benchmark (quadratic velocity profile
> & linear pressure) using the direct solver are:
>
> L2-errors in v & p: Q1Q1 + PSPG
>
> --------------------------------------------------------------------------------------------------
> lvl 0: | e_v = 1 | eoc_v = 0 | e_p = 1
> | eoc_p = 0 |
>
> --------------------------------------------------------------------------------------------------
> lvl 1: | e_v = 0.531463 | eoc_v = 0.91196 | e_p = 0.479585
> | eoc_p = 1.06014 |
>
> --------------------------------------------------------------------------------------------------
> lvl 2: | e_v = 0.218196 | eoc_v = 1.28434 | e_p = 0.208925
> | eoc_p = 1.1988 |
>
> --------------------------------------------------------------------------------------------------
> lvl 3: | e_v = 0.0661014 | eoc_v = 1.72287 | e_p = 0.0674415
> | eoc_p = 1.63128 |
>
> --------------------------------------------------------------------------------------------------
> lvl 4: | e_v = 0.0176683 | eoc_v = 1.90352 | e_p = 0.0194532
> | eoc_p = 1.79363 |
>
> --------------------------------------------------------------------------------------------------
> lvl 5: | e_v = 0.00452588 | eoc_v = 1.96489 | e_p = 0.00549756
> | eoc_p = 1.82315 |
>
> --------------------------------------------------------------------------------------------------
> lvl 6: | e_v = 0.00114238 | eoc_v = 1.98616 | e_p = 0.00158448
> | eoc_p = 1.79478 |
>
> --------------------------------------------------------------------------------------------------
> lvl 7: | e_v = 0.000286765 | eoc_v = 1.9941 | e_p = 0.000475072
> | eoc_p = 1.73779 |
>
> --------------------------------------------------------------------------------------------------
> lvl 8: | e_v = 7.1824e-05 | eoc_v = 1.99733 | e_p = 0.000149148
> | eoc_p = 1.67141 |
>
> --------------------------------------------------------------------------------------------------
> lvl 9: | e_v = 1.79717e-05 | eoc_v = 1.99874 | e_p = 4.88044e-05
> | eoc_p = 1.61166 |
>
> --------------------------------------------------------------------------------------------------
> lvl 10: | e_v = 4.49483e-06 | eoc_v = 1.99939 | e_p = 1.64706e-05
> | eoc_p = 1.56712 |
>
> --------------------------------------------------------------------------------------------------
> where e_v & e_p are relative errors, i.e., int(p-ph)dx/int(p)dx
>
> The obtained rates using stable Taylor-Hood Q2Q1 are:
> L2-errors in v & p: Q2Q1
>
> --------------------------------------------------------------------------------------------------
> lvl 0: | e_v = 2.20205e-16 | eoc_v = 0 | e_p = 3.58011e-16
> | eoc_p = 0 |
>
> --------------------------------------------------------------------------------------------------
> lvl 1: | e_v = 2.1471e-16 | eoc_v = 0.0364598 | e_p = 1.02455e-15
> | eoc_p = -1.51692 |
>
> --------------------------------------------------------------------------------------------------
> lvl 2: | e_v = 4.23849e-16 | eoc_v = -0.98116 | e_p = 4.1472e-16
> | eoc_p = 1.30479 |
>
> --------------------------------------------------------------------------------------------------
> lvl 3: | e_v = 6.9843e-16 | eoc_v = -0.720565 | e_p = 4.49311e-15
> | eoc_p = -3.4375 |
>
> --------------------------------------------------------------------------------------------------
> lvl 4: | e_v = 1.53993e-15 | eoc_v = -1.14068 | e_p = 9.82242e-15
> | eoc_p = -1.12837 |
>
> --------------------------------------------------------------------------------------------------
> lvl 5: | e_v = 7.60182e-15 | eoc_v = -2.30348 | e_p = 5.44953e-14
> | eoc_p = -2.47198 |
>
> --------------------------------------------------------------------------------------------------
> lvl 6: | e_v = 1.60118e-14 | eoc_v = -1.07472 | e_p = 1.86328e-13
> | eoc_p = -1.77364 |
>
> --------------------------------------------------------------------------------------------------
> ... which shows, that I get the exact solution, which is expected.
>
> Using Q2Q2 + PSPG i get:
> L2-errors in v & p: Q2Q2 + PSPG
>
> --------------------------------------------------------------------------------------------------
> lvl 0: | e_v = 3.5608e-16 | eoc_v = 0 | e_p = 6.5321e-16
> | eoc_p = 0 |
>
> --------------------------------------------------------------------------------------------------
> lvl 1: | e_v = 7.51502e-16 | eoc_v = -1.07758 | e_p = 1.59993e-15
> | eoc_p = -1.29239 |
>
> --------------------------------------------------------------------------------------------------
> lvl 2: | e_v = 7.47141e-16 | eoc_v = 0.00839586 | e_p = 7.92418e-16
> | eoc_p = 1.01367 |
>
> --------------------------------------------------------------------------------------------------
> lvl 3: | e_v = 1.36773e-15 | eoc_v = -0.872327 | e_p = 4.87009e-15
> | eoc_p = -2.61961 |
>
> --------------------------------------------------------------------------------------------------
> lvl 4: | e_v = 2.16714e-15 | eoc_v = -0.664013 | e_p = 2.01012e-14
> | eoc_p = -2.04526 |
>
> --------------------------------------------------------------------------------------------------
> lvl 5: | e_v = 1.14488e-14 | eoc_v = -2.40134 | e_p = 2.04823e-14
> | eoc_p = -0.0270959 |
>
> --------------------------------------------------------------------------------------------------
> lvl 6: | e_v = 1.15543e-14 | eoc_v = -0.0132335 | e_p = 1.7867e-13
> | eoc_p = -3.12484 |
>
> --------------------------------------------------------------------------------------------------
> also (!) giving me the exact solution for all meshes used, but actually
> one has to use the PSPG here.
>
> The only suspicious thing I found, is that when I use in the preparation
> of the testfunctions:
>
> if (get_laplace_residual)
> {
> Tensor<3,dim> phi_hessian_v;
> phi_hessian_v =
> fe_hessians[velocities].hessian(k,q);
>
> for (int d = 0; d < dim; ++d)
> {
> phi_laplacians_v[k][d] =
> trace(phi_hessian_v[d]);
> }
> if (phi_hessian_v.norm() > 1e-14)
> std::cout << "|H| = " << std::scientific <<
> phi_hessian_v.norm() << std::endl;
> }
>
> There is actually a non-zero norm printed during execution. Linear shape
> functions on a rectangular(!) grid should
> despite being subject to a bilinear mapping still give zero second
> derivatives, right?
> Since my understanding of the hessian() function is somewhat limited, I
> also tried with
>
> if (get_laplace_residual)
> {
> phi_hessian_v = 0;
> phi_hessian_v =
> fe_hessians[velocities].hessian(k,q);
>
> for (int d = 0; d < dim; ++d)
> {
> phi_laplacians_v[k][d] =
> trace(phi_hessian_v[d]); // ###
>
> if (fe.system_to_component_index(k).first <
> dim)
> phi_laplacians_v[k][d] =
> trace(fe_hessians.shape_hessian_component(k,q,d));
> else
> phi_laplacians_v[k][d] = 0.0;
>
> Tensor<1,dim> diff;
> diff = 0;
> diff = trace(phi_hessian_v[d]) -
> trace(fe_hessians.shape_hessian_component(k,q,d));
>
> if (diff.norm () > 1e-14)
> std::cout << "|d| = " << std::scientific
> << diff.norm() << std::endl;
>
> }
> if (phi_hessian_v.norm() > 1e-14)
> std::cout << "|H| = " << std::scientific <<
> phi_hessian_v.norm() << std::endl;
> }
>
> Showing, that I get same (correct?) results using any function to compute
> the laplacian, since nothing is printed to screen when executing.
>
> Can someone confirm, that the use of these 2 functions is correct?
> I (almost) use the same call as the proposed code 'lethe'.
>
> Does the Q2Q2+PSPG not show, that the implementation is correct? Is the
> (internal) treatment of linear shape functions different?
> The non-zero second derivatives are still somewhat suspicious to me.
>
> If all else fails, I am gonna strip off all the extra stuff I have in my
> flow solver and post a simple Q(p)xQ(q) + pspg.
>
> Kind regards & 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/bb1a58f6-99fc-4fdf-a678-9e9c0716fd9c%40googlegroups.com.
