solving stabilized Stokes

[Richard Schussnig]

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

[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

[Wolfgang Bangerth]

On 2/14/20 5:34 AM, Richard Schussnig wrote:
>
> 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?

We have all spent much time trying to figure out these sorts of issues
:-) My usual approach is to switch to a direct solver to eliminate the
possibility that the problem lies with the solver. Just speaking about,
say, step-22, there is also the issue that the iterative solver actually
just ignores the pressure-pressure block of the matrix (I think that we
put a mass matrix in there as that is used in the preconditioner -- or
at least that's what we used to do). So it's useful to just remove
everything that could be problematic, and the iterative solvers are
definitely a point in case.

The other thing I usually check is what happens if I increase the order
of quadrature.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Richard Schussnig]

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

[Richard Schussnig]

Dear Wolfgang,

Thanks also to you for taking the time, I really appreciate it!

I already tested the integration, increasing up by order+1, +2, +3 ... gave identical results.

Now it seems, that I should introduce some switches to a direct solver! This might give me a hint, if I need some rescaling of some rows!

I will report back, if that helped!

Kind regards,

Richard

[Bruno Blais]

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

[Richard Schussnig]

Dear Bruno,

thanks again for bothering!

I will try to do that, but it is a bit involved, since I was setting up block-matrices, which i need to change!

For pure Dirichlet, enclosed flow problems I did the same as for step-55f, which is basically nothing
(the iterative solver handles the pressure scaling) but after solving rescale the pressure to mean value of zero,
which is working like a charm for q2q1 for example!

I will report back, when I get some more insight using a direct solver!

Kind regards,
Richard

[Richard Schussnig]

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

[Bruno Blais]

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

[Wolfgang Bangerth]

On 3/11/20 8:54 AM, Bruno Blais wrote:
> 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./

If I understand the chain of emails correctly, this is the question you
have: Using a Q1/Q1 element with GLS/PSPG, what convergence rates does
one expect.

You seem to suggest that it should be 2 for both velocity and pressure.
But is that true? What does the theory say? I would actually be quite
surprised if that's what you get for these kinds of problems. My gut
feeling is that you should expect 2=p+1 for the velocity, but something
between p and p+1 for the pressure. That would be consistent with other
stabilized approaches.

[Bruno Blais]

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

[Bruno Blais]

Dear Wolfgang,

I spent some time re-reading the theory and you are right, nothing shows that the convergence rate should be conserved when we have stabilized equal order elements. However it is interesting to note that when we use stabilization and we revert to LBB stable elements (Q2-Q1 for instance) the right order is recovered.

It seems that Richard was expecting that the order would be conserved however when PSPG stabilization is used. In our case, at least,I do not manage to reproduce that (and that is across the usage of two codes using GLS implementations).

[Wolfgang Bangerth]

On 3/11/20 11:25 AM, Bruno Blais wrote:
> I spent some time re-reading the theory and you are right, nothing shows
> that the convergence rate should be conserved when we have stabilized
> equal order elements. However it is interesting to note that when we use
> stabilization and we revert to LBB stable elements (Q2-Q1 for instance)
> the right order is recovered.

Yes, that is correct.

But I still think that the best you can hope for with equal-order
(Qk-Qk) discretizations is optimal order for 'u', but not for 'p'. My
gut reaction is that that is going to be true for any stabilization
method. In other words, your numerical experiments do not surprise me --
they should what I would have naively expected.

[Richard Schussnig]

Dear Wolfgang, dear Bruno,

Thank you very much again for your time & effort!

I think one should leave grad-div stabilization out of the discussion, since it is only used as a counter-measure. It has nothing to do with inf-sup stability,

but is rather used to additionally penalize violations of continuity of mass in the case of stabilized formulations through a penalty-term in the balance of linear momentum.

I also checked up on the literature in V. John s book from 2016: FEM for incompr flow problems
and what I found is, that in example 4.100 for stabilized Stokes shows nicely, that the velocity is converging quadratically

& the order of the pressure convergence depends on viscosity and is in [1,2] for q1q1 elements (nu = 1e-10 giving ~2).

So summing up, I was simply expecting too much it seems!

Best wishes & kind regards from Austria,
Richard