In addition to the above results. I ran the same MMS with the GLS
stabilized solver but using Q2-Q1 and Q2-Q2.
What i obtain for Q2-Q1 is exactly what you would expect:
*Velocity*
Cell error -
64 1.4972e-02 -
256 1.9438e-03 2.95
1024 2.4542e-04 2.99
4096 3.0755e-05 3.00
*Pressure*
cell error -
64 0.0232558 -
256 0.00355176 2.55
1024 0.000831086 2.06
4096 0.000206113 2.01
And for Q2-Q2 the results become similar to Q1-Q1 in that the pressure
convergence is lost:
Velocity
cells error
64 1.5715e-02 -
256 1.9762e-03 2.99
1024 2.4655e-04 3.00
4096 3.0792e-05 3.00
16384 3.8480e-06 3.00
Pressure
cells error
64 0.0796909 2.03
256 0.0192026 2.01
1024 0.00476045 2.00
4096 0.00118731 2.00
16384 0.000296642 2.00
On Wednesday, 11 March 2020 10:54:21 UTC-4, Bruno Blais wrote:
>
> 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/be4aa671-7937-4be9-8b7c-1b637039b46c%40googlegroups.com.
