Pressure-correction scheme: Behaviour of the H1 norm on the Taylor-Green vortex

[jose.a...@gmail.com]

Hello everyone,

I am working on a parallel solver based upon the incremental rotational pressure-scheme  as seen in the overview paper of Guermond and implemented in step-35.  To verify the code, the following problems were replicated

• Numerical test proposed in section 3.7.2 of Guermond's paper. It is based on the method of manufactured solutions. The domain is the unit square and a circle of radius 0.5. Dirichlet boundary conditions are imposed on the whole boundary for the velocity. The pressure field is constrained by setting its mean value to zero.
  

• The Taylor-Green vortex. The domain is the square (0,1)^2. The velocity and pressure field are constrained by periodic boundary conditions. Furthermore, the pressure field is constrained by setting its mean value to zero.

In order to compare the numerical pressure solution to its analytical counterpart, the mean value of the numerical pressure field is made to match that of the analytical one at the end of each time step.

The results for the Guermond's numerical tests for the square domain

                               Velocity convergence table
==============================================================================================
level    dt         cells     dofs       hmax           L2                             H1                             Linfty      
    7 1.00e-01  8572 132098 1.10e-02 8.557567e-03     -    6.456461e-02      -   1.103922e-02     -
    7 2.50e-02  8572 132098 1.10e-02 5.354964e-04 -2.00 5.090698e-03 -1.83 6.822255e-04 -2.01
    7 6.25e-03  8572 132098 1.10e-02 3.358435e-05 -2.00 4.221552e-04 -1.80 4.245914e-05 -2.00
    7 1.56e-03  8572 132098 1.10e-02 2.105764e-06 -2.00 3.632448e-05 -1.77 2.656204e-06 -2.00

                               Pressure convergence table
==============================================================================================
level    dt         cells   dofs       hmax           L2                             H1                             Linfty      
    7 1.00e-01  8572 16641 1.10e-02 5.543160e-03     -    2.687030e-02     -    3.649990e-02     -
    7 2.50e-02  8572 16641 1.10e-02 3.394874e-04 -2.01 3.684615e-03 -1.43 4.209359e-03 -1.56
    7 6.25e-03  8572 16641 1.10e-02 2.126561e-05 -2.00 2.738788e-03 -0.21 3.958393e-04 -1.71
    7 1.56e-03  8572 16641 1.10e-02 3.049911e-06 -1.40 2.727803e-03 -0.00 2.751900e-05 -1.92

and for the circle domain

                               Velocity convergence table
==============================================================================================
level    dt         cells      dofs      hmax              L2                              H1                        Linfty      
    6 1.00e-01 10979 164354 1.34e-02 7.593191e-03     -    5.404824e-02     -    9.465067e-03     -
    6 2.50e-02 10979 164354 1.34e-02 4.879026e-04 -1.98 4.298417e-03 -1.83 6.072627e-04 -1.98
    6 6.25e-03 10979 164354 1.34e-02 3.069737e-05 -2.00 3.539538e-04 -1.80 3.808221e-05 -2.00
    6 1.56e-03 10979 164354 1.34e-02 1.926069e-06 -2.00 3.267549e-05 -1.72 2.500277e-06 -1.96

                               Pressure convergence table
==============================================================================================
level    dt         cells      dofs      hmax              L2                              H1                        Linfty      
    6 1.00e-01 10979 20609 1.34e-02 2.855089e-03     -    1.235457e-02     -    7.053574e-03     -
    6 2.50e-02 10979 20609 1.34e-02 1.983995e-04 -1.92 2.829075e-03 -1.06 4.839667e-04 -1.93
    6 6.25e-03 10979 20609 1.34e-02 1.308016e-05 -1.96 2.705119e-03 -0.03 3.305446e-05 -1.94
    6 1.56e-03 10979 20609 1.34e-02 3.290324e-06 -1.00 2.704622e-03 -0.00 2.323934e-05 -0.25

The Fig. 4 of the paper presents the log-log plot of the infinity norm of the pressure vs time step size using 8 time step sizes. On the square domain it shows convergence order of 1.6 on the square domain. Whereas on the circle domain it shows a convergence order of 2.0 up to a time step size of 2e-3. The last time step shown in Fig. 4 hints the reaching of the asymptotic value of the error.

The results for the Taylor-Green vortex problem

                               Velocity convergence table
==============================================================================================
level    dt         cells      dofs       hmax           L2                              H1                       Linfty      
    8 1.00e-01 17980 526338 5.52e-03 3.057096e-02     -    4.354711e-01     -    4.998748e-02     -
    8 5.00e-02 17980 526338 5.52e-03 9.313187e-03 -1.71 1.289422e-01 -1.76 1.550793e-02 -1.69
    8 2.50e-02 17980 526338 5.52e-03 2.493096e-03 -1.90 3.395720e-02 -1.92 4.112861e-03 -1.91
    8 1.25e-02 17980 526338 5.52e-03 6.411600e-04 -1.96 8.668379e-03 -1.97 1.052791e-03 -1.97
    8 6.25e-03 17980 526338 5.52e-03 1.623733e-04 -1.98 2.190470e-03 -1.98 2.661297e-04 -1.98
    8 3.12e-03 17980 526338 5.52e-03 4.084716e-05 -1.99 5.601806e-04 -1.97 6.694580e-05 -1.99

                               Pressure convergence table
==============================================================================================
level    dt         cells      dofs       hmax           L2                              H1                       Linfty      
    8 1.00e-01 17980 66049 5.52e-03 1.582475e-02     -    2.265999e-01     -    4.902018e-02     -
    8 5.00e-02 17980 66049 5.52e-03 4.698715e-03 -1.75 6.899043e-02 -1.72 1.415426e-02 -1.79
    8 2.50e-02 17980 66049 5.52e-03 1.249606e-03 -1.91 2.664252e-02 -1.37 3.746148e-03 -1.92
    8 1.25e-02 17980 66049 5.52e-03 3.211439e-04 -1.96 2.054385e-02 -0.38 9.905812e-04 -1.92
    8 6.25e-03 17980 66049 5.52e-03 8.191443e-05 -1.97 2.004909e-02 -0.04 2.840482e-04 -1.80
    8 3.12e-03 17980 66049 5.52e-03 2.280548e-05 -1.84 2.000950e-02 -0.00 1.051814e-04 -1.43

While the velocity looks alright the H1 norm of the pressure locks pretty fast and to a high value, while the L2 and infinity norm indicate a decline in the convergence rate.

I performed then a spatial convergence test on both problems, which delivered for the Guermond problem

                               Velocity convergence table
==============================================================================================
level    dt         cells    dofs      hmax           L2                               H1                          Linfty      
    5 2.00e-02   361   8450     4.42e-02 3.390456e-04     -    4.005312e-03     -    4.264292e-04     -
    6 2.00e-02  1228  33282   2.21e-02 3.379001e-04 -0.00 3.928219e-03 -0.03 4.266117e-04  0.00
    7 2.00e-02  4495  132098 1.10e-02 3.375465e-04 -0.00 3.881829e-03 -0.02 4.266296e-04  0.00
    8 2.00e-02 17170 526338 5.52e-03 3.374516e-04 -0.00 3.866594e-03 -0.01 4.265951e-04 -0.00

                               Pressure convergence table
==============================================================================================
level    dt         cells    dofs      hmax           L2                               H1                          Linfty      
    5 2.00e-02   361   1089   4.42e-02 3.021009e-04     -    1.117284e-02     -    1.817554e-03    -
    6 2.00e-02  1228  4225   2.21e-02 2.991258e-04 -0.01 6.139539e-03 -0.86 2.426530e-03 0.42
    7 2.00e-02  4495  16641 1.10e-02 2.988485e-04 -0.00 4.104065e-03 -0.58 2.943609e-03 0.28
    8 2.00e-02 17170 66049 5.52e-03 2.988060e-04 -0.00 3.494146e-03 -0.23 3.440934e-03 0.23

and for the Taylor-green vortex

                               Velocity convergence table
==============================================================================================
level    dt         cells    dofs      hmax           L2                               H1                          Linfty      
    5 1.00e-02   472   8450     4.42e-02 4.176710e-04     -    2.838026e-02     -    7.339885e-04     -
    6 1.00e-02  1444  33282   2.21e-02 4.135497e-04 -0.01 6.707143e-03 -2.08 6.833520e-04 -0.10
    7 1.00e-02  4912  132098 1.10e-02 4.127579e-04 -0.00 5.601192e-03 -0.26 6.776613e-04 -0.01
    8 1.00e-02 17980 526338 5.52e-03 4.125038e-04 -0.00 5.569924e-03 -0.01 6.767872e-04 -0.00

                               Pressure convergence table
==============================================================================================
level    dt         cells    dofs      hmax           L2                               H1                          Linfty      
    5 1.00e-02   472   1089   4.42e-02 6.878170e-04     -     1.609285e-01     -     3.492529e-03     -
    6 1.00e-02  1444  4225   2.21e-02 2.630181e-04 -1.39 8.027648e-02 -1.00 1.329719e-03 -1.39
    7 1.00e-02  4912  16641 1.10e-02 2.105590e-04 -0.32 4.018553e-02 -1.00 7.874550e-04 -0.76
    8 1.00e-02 17980 66049 5.52e-03 2.067516e-04 -0.03 2.024105e-02 -0.99 6.528251e-04 -0.27

I was expecting the error to lock everywhere to the asymptotic value of the error set by the time step but was surprised to see that on the Taylor-Green vortex the refinement level sets a convergence rate of 1 on the pressure H1-norm and has a weaker influence on the L2 and infinity pressure norms and the H1-norm of the velocity. The Guermond problem shows a similar behaviour but milder. My question to anyone with experience on projection schemes and the Taylor-Green vortex, is whether this behaviour is problem dependent, is intrinsic of the pressure-correction scheme, both or am I overseeing an implementation/coding error and this behaviour should not occurr at all?

Extra info: A Taylor-Hood element was implemented. The convergence rates shown in the tables are the experimental order of convergence using the maximal circumradius of the triangulation (hmax) as reference column. The level column indicates the number of global refinements. Please ignore the cells column as on some tables the number of cells were obtained using the wrong deal.ii method (n_active_cells instead of n_global_active_cells).

Cheers,

Jose

[Wolfgang Bangerth]

On 10/18/20 10:21 AM, jose.a...@gmail.com wrote:
>
> I was expecting the error to lock everywhere to the asymptotic value of the
> error set by the time step but was surprised to see that on the Taylor-Green
> vortex the refinement level sets a convergence rate of 1 on the pressure
> H1-norm and has a weaker influence on the L2 and infinity pressure norms and
> the H1-norm of the velocity. The Guermond problem shows a similar behaviour
> but milder. My question to anyone with experience on projection schemes and
> the Taylor-Green vortex, is whether this behaviour is problem dependent, is
> intrinsic of the pressure-correction scheme, both or am I overseeing an
> implementation/coding error and this behaviour should not occurr at all?

Jose -- your question is very specific, and I'm in doubt whether there is
anyone around with the exact piece of information you're looking for.

I'm not entirely clear about what your question is. Are you seeing convergence
rates that are too low or too large? It is not uncommon to have cases where a
scheme converges too fast (the convergence rate is too large); this is
typically the case because the solution has a symmetry.

Best
W.

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

[Martin Kronbichler]

Dear Jose,

While I have not experimented in detail with the step-35 program, we have done extensive studies on similar problems in https://doi.org/10.1016/j.jcp.2017.09.031 (or https://arxiv.org/abs/1706.09252 for an earlier preprint version of the same manuscript) including the pressure correction scheme. While the spatial discretization is DG where some of the issues are , the experience from our experiments suggests that the pressure correction scheem should not behave too differently from other time discretization schemes with similar ingredients (say BDF-2 on the fully coupled scheme). In other words, I would expect second order convergence in the pressure for Taylor-Hood elements. I should note that there are some subtle issues with boundary conditions in projection schemes, so I cannot exclude some hidden problems with the step-35 implementation or the way you set up the experiments.

Best,
Martin

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f10a98c1-9ec7-4944-a850-97cd6b43f148n%40googlegroups.com.

[blais...@gmail.com]

Hello Jose,

I wish I could help, but I second Wolfgang's question.

Is your code available somewhere? I would be glad to take a look at it and compare the solutions for the same problems using different formulations. I would expect that if you fix the issue with boundary conditions (those described in the Guermond paper, that is the  "pressure boundary layer") then you would recover exactly what you should get with traditional schemes using Taylor-Hood element (as Martin discussed).

On another note, I remember having a discussion about this with Timo Heister at the deal.II workshop in 2019. Maybe Timo has ideas on this? I know he is quite the expert on algorithms to solve the Stokes / Navier-Stokes equations (e.g. his paper on the grad-div scheme, etc.)

Sorry for not being to help more.

Best

Bruno

[richard....@gmail.com]

Hi Jose,

just a thought from my side:

if the velocity error you have is low enough, the somewhat time-independent PPE you solve given that velocity,

you might get high enough rates up to a certain point -> and that point might lie below the error you see on those 3 levels.

So, try to go for smaller timesteps (keeping the Re the same) and use more spatial refinement levels.

In general I would also recommend changing the Reynold's number a bit around and see what happens - maybe it is an effect that is limited to low Re?

Anyways, having a bigger convergence rate than expected is a nice problem to have, isn't it? ; )
-I would not think it is caused by a bug given the other rates looking just as expected!

All the best,
Richard

[Timo Heister]

> On another note, I remember having a discussion about this with Timo Heister at the deal.II workshop in 2019. Maybe Timo has ideas on this? I know he is quite the expert on algorithms to solve the Stokes / Navier-Stokes equations (e.g. his paper on the grad-div scheme, etc.)

I am happy to give further comments, but -- like Wolfgang -- I don't
quite understand what the precise question is. That said:

1. With projection schemes you will need to be careful about pressure
boundary layers. A good starting point might be the Elman, Silvester,
Wathen book.
2. Specific numerical test setups can be more or less sensitive to
this fact (size of pressure error vs velocity error, smoothness of
solutions in time, specific behavior on the boundary, ...)


--
Timo Heister
http://www.math.clemson.edu/~heister/

[Jose Lara]

Hello Wolfgang, Marthin, Bruno, Richard and Timo,

 

'm not entirely clear about what your question is. Are you seeing convergence
rates that are too low or too large? It is not uncommon to have cases where a
scheme converges too fast (the convergence rate is too large); this is
typically the case because the solution has a symmetry.

Best
  W.

Apologies for not explaining myself clearer. I will try it again:

From my understanding, the lowest bound of the error on each norm is set either by the spatial or the temporal discretization. I was kind of expecting that the L2- and H1-Norms share a similar spatial and time dependence, i.e. that each field reaches its lowest bound simultaneously, and that they do so with a similar convergence rate evolution. Stated differently, they start with an order of convergence which remains constant for a given time step range. After reaching a small time step size, the convergence order tends to zero as the lowest bound of the error is reached.

From the tests' results I can see that H1-Norm of the pressure has a considerably stronger spatial dependence than the velocity, as it reaches its lowest bound while the velocity still has a constant convergence order. This behaviour is also seen in the L2- and Linfty-Norm but in a much more milder scale, as seen in the spatial convergence test. My question is if this stronger spatial dependency of the pressure is problem dependent or if it is intrinsic to the pressure-correction scheme.

While I have not experimented in detail with the step-35 program, we have done extensive studies on similar problems in https://doi.org/10.1016/j.jcp.2017.09.031 (or https://arxiv.org/abs/1706.09252 for an earlier preprint version of the same manuscript) including the pressure correction scheme. While the spatial discretization is DG where some of the issues are , the experience from our experiments suggests that the pressure correction scheem should not behave too differently from other time discretization schemes with similar ingredients (say BDF-2 on the fully coupled scheme). In other words, I would expect second order convergence in the pressure for Taylor-Hood elements. I should note that there are some subtle issues with boundary conditions in projection schemes, so I cannot exclude some hidden problems with the step-35 implementation or the way you set up the experiments.

Best,
Martin

Thanks for the manuscript, there I notice that the results shown are those of the behaviour of the L2-Norm. My finite element implementation behaves similarly in the L2-Norm to the convergence rates in your paper (BFD2 leads to a 2nd order convergence on both velocity and, for the most part, on pressure). Did you also analyse the H1-Norm by any chance? There is where I see the stronger spatial dependency of the pressure.

On another note, it caught my eye that you split the Neumann boundary conditions. I have not done tests with them yet, but what is the benefit of doing this or why is it necessary? Furthermore, my next step would be the DFG 2D-2 benchmark. There, you computed the traction force using the symmetric gradient instead of the normal gradient. While your formulation would be for me the correct formulation, as the stress tensor is so defined, on the benchmark papers they use the normal gradient. This has caused me some confusion, as to which formulation should I implement for the benchmark.

Hello Jose,

I wish I could help, but I second Wolfgang's question.

Is your code available somewhere? I would be glad to take a look at it and compare the solutions for the same problems using different formulations. I would expect that if you fix the issue with boundary conditions (those described in the Guermond paper, that is the  "pressure boundary layer") then you would recover exactly what you should get with traditional schemes using Taylor-Hood element (as Martin discussed).

I tried reformulating the question above. Hopefully it is clearer now. I will clean the code up and get back to you. The pressure-correction scheme I am using is the incremental rotational, which has the smallest error caused due to the boundary layer. Could the boundary layer in this case still cause such a strong influence?

if the velocity error you have is low enough, the somewhat time-independent PPE you solve given that velocity,

you might get high enough rates up to a certain point -> and that point might lie below the error you see on those 3 levels.

So, try to go for smaller timesteps (keeping the Re the same) and use more spatial refinement levels.

In general I would also recommend changing the Reynold's number a bit around and see what happens - maybe it is an effect that is limited to low Re?

Anyways, having a bigger convergence rate than expected is a nice problem to have, isn't it? ; )
-I would not think it is caused by a bug given the other rates looking just as expected!

All the best,
Richard

Yes, this is on my todo list. Until now I have been working on my laptop but now I got access to a cluster so I can expand the scope of the tests. For the Guermond problem I had been using Re = 100 and for the Taylor-Green vortex Re = 5 and as you propose I should investigate of this effect is driven by the Reynolds number. Thanks for the suggestion!

I am happy to give further comments, but -- like Wolfgang -- I don't
quite understand what the precise question is. That said:

1. With projection schemes you will need to be careful about pressure
boundary layers. A good starting point might be the Elman, Silvester,
Wathen book.
2. Specific numerical test setups can be more or less sensitive to
this fact (size of pressure error vs velocity error, smoothness of
solutions in time, specific behavior on the boundary, ...)


--
Timo Heister

Sorry for not being clearer. I reformulated the question at the beginning of the E-Mail. Thanks for the book suggestion, I will get on reading.

The pressure-correction scheme I am using is the incremental rotational, for which the intrinsic boundary layer error is the smallest of all the schemes. Furthermore, for the Taylor-Green vortex I am using periodic boundary conditions on all boundaries which rules out the corner singularities which plague the pressure-correction scheme. The solution in itself is smooth, so I had not thought the error of the boundary layer could have had such an effect on the H1-Norm. Interestingly, while the Guermond problem has Dirichlet boundary conditions and a non-smooth boundary, it also shows a better behaviour on the H1-Norm.

Cheers,

Jose

[Timo Heister]

> From my understanding, the lowest bound of the error on each norm is set either by the spatial or the temporal discretization. I was kind of expecting that the L2- and H1-Norms share a similar spatial and time dependence, i.e. that each field reaches its lowest bound simultaneously, and that they do so with a similar convergence rate evolution. Stated differently, they start with an order of convergence which remains constant for a given time step range. After reaching a small time step size, the convergence order tends to zero as the lowest bound of the error is reached.

I don't support this idea. If you take the Taylor-Green vortex for
example, the velocity decays with exp(-2 nu t), while the pressure
decays with the square of that term. Why do you expect error from your
pressure-correction scheme and your error from the time discretization
to converge in the same way? Note that they are not completely
independent of course.

> Furthermore, for the Taylor-Green vortex I am using periodic boundary conditions on all boundaries which rules out the corner singularities which plague the pressure-correction scheme. The solution in itself is smooth, so I had not thought the error of the boundary layer could have had such an effect on the H1-Norm.

I am not convinced that this completely eliminates all influence of
the pressure-correction scheme. I assume that it still gives an O(dt)
additional error (or maybe something higher order depending on your
scheme).

[Jose Lara]

I don't support this idea. If you take the Taylor-Green vortex for
example, the velocity decays with exp(-2 nu t), while the pressure
decays with the square of that term. Why do you expect error from your
pressure-correction scheme and your error from the time discretization
to converge in the same way? Note that they are not completely
independent of course.

The assumption was naive on my part. I had not considered that the pressure-correction scheme could introduce such a strong spatial dependency on the error saturation of the pressure field alone.

I am not convinced that this completely eliminates all influence of
the pressure-correction scheme. I assume that it still gives an O(dt)
additional error (or maybe something higher order depending on your
scheme).

Yes, you are right. If those factors would reduce/eliminate the influence of the pressure-correction scheme, the Guermond problem would then showcase a lower convergence order and/or faster saturation than the Taylor-Green vortex, which is not the case.

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/0MkM-7k6tjo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAMRj59GKEjXp0t_5EED8myCmi0pcdwjA7GgZ2mb5vL8DpDQCCA%40mail.gmail.com.

[blais...@gmail.com]

Hello Jose,

I wish I could help, but I second Wolfgang's question.

Is your code available somewhere? I would be glad to take a look at it and compare the solutions for the same problems using different formulations. I would expect that if you fix the issue with boundary conditions (those described in the Guermond paper, that is the  "pressure boundary layer") then you would recover exactly what you should get with traditional schemes using Taylor-Hood element (as Martin discussed).

I tried reformulating the question above. Hopefully it is clearer now. I will clean the code up and get back to you. The pressure-correction scheme I am using is the incremental rotational, which has the smallest error caused due to the boundary layer. Could the boundary layer in this case still cause such a strong influence? 

The TGV case you are simulating only has periodic boundary conditions. Artificial pressure boundary layer only appear close to walls where you apply a boundary condition on the pressure correction (this is discussed in the Guermond paper). If you do not have a wall (or any Dirichlet BC on the velocity), you should not get any artificial pressure boundary layer. If you are seeing the introduction of an additionnal error in time, it means it is introduce by the corrector step, but it is not related to the pressure boundary layer.

Otherwise, I second Timo's answer on the convergence of the error.

From my experience, the best way to analyze the TGV is to try to decouple the errors in time and space as much as possible by analyzing the evolution of the error in time with a large time step and a very fine mesh and analyzing the error in space with a very fine time-step and coarser meshes. Otherwise the two errors overlap and it is very hard to draw any sort of meaningful conclusions. This is especially true if you use high-order in time (BDF2, BDF3, SDIRK22, etc.).

 

[blais...@gmail.com]

Hello Jose,

I hope you are well.

I was wondering if you had managed to get your code on github or any sharing service? I would be really interested in comparing how it behaves with monolothic FEM approaches. This is something that is very relevant to our research (and on which I would be glad to collaborate :) )

[richard....@gmail.com]

Hi Jose,
I am not too sure if you will actually read this, but I was somehow unable to extract your email from this thread?!
Thing is, me and my colleague are organizing a Minisymposium at the YIC 2021 in Valencia, Spain, 7th-9th of July and your work seems to be fitting pretty well into our scope!

See, e.g., these links:

https://yic2021.upv.es/

https://yic2021.upv.es/wp-content/uploads/2020/11/15-1.pdf
In case you are interested (or in fact anyone reading this is up for a nice trip to Valencia, meeting some new people and in general having a blast), please just let me know and I would be glad to invite you!

Kind regards,
Richard