Pressure does not converge for Stokes problem with RT elements and inhomogeneous Dirichlet boundary conditions

[HC Zhang]

Hello everyone,

I am solving the Stokes equations using Raviart–Thomas elements in deal.II. When the boundary conditions are homogeneous Dirichlet, everything works fine and I observe the expected convergence behavior.

However, when I switch to inhomogeneous Dirichlet boundary conditions, I encounter problems: although I apply ``VectorTools::project_boundary_values_div_conforming(...)`` to set the boundary values, the pressure error does not converge.

Has anyone experienced a similar issue, or could point me to the right way of handling inhomogeneous Dirichlet conditions for the Stokes problem with RT elements? Am I missing an additional step to make the pressure converge?

Any advice or references would be very much appreciated.

Thanks in advance!

[Wolfgang Bangerth]

HC: I think there isn't enough information here. For example, I'm not
even sure what variable you're providing boundary values.

Regardless, what have you tried already? When you compare the computed
and the exact solution, which variable differs? The velocity or the
pressure? If you look at a visualization of the solution, in which
specific ways does the computed solution differ from the expected one?
Is it offset by a constant? Are the boundary values completely wrong?
For example, do you expect the boundary values to be g(x) but they are
zero? When you say "the pressure error does not converge", does that
mean that it doesn't converge to the right value, or that it in fact
*diverges*?

I often find it useful to be specific in which ways the solution does
not match my expectations when trying to find errors.

Best
W.

[blais...@gmail.com]

In addition to what Wolfgang mentioned, I was wondering how you were managing the pressure constant in this case. This may severely hinder your "apparent" convergence, especially for pressure.

[HC Zhang]

Deal all, 
Thank you very much for the quick replies — you are right that my previous description lacked important details. I’ll summarize the setup more precisely and then show the numerical behaviour.

Problem / setup:

• I follow the paper “Multigrid methods for H(div)-conforming discontinuous Galerkin methods for the Stokes equations”.

• Finite elements: RT_k-DGQ_k (so an H(div)-conforming RT velocity, DG pressure).

• Solver: FGMRES for the global linear system; multigrid used as a preconditioner; the multigrid smoother is a multilevel Schwarz method.

• PDE form: I use the stress-tensor form of the Stokes equations (same form as in deal.II step-22). But I think under the divergence-free condition these two forms should be equivalent.

• Boundary conditions: I implement inhomogeneous Dirichlet boundary conditions. Nitsche’s method is applied in the weak form for enforcing the velocity BCs. In setup_system() I also apply

  `` VectorTools::project_boundary_values_div_conforming(...) `` to the velocity block.

• Pressure: I did not apply any extra constraint or projection in the weak form — the pressure block is assembled directly from the weak formulation (no extra stabilization there).

About the pressure constant / mean value:
Before computing the error I tried to remove the mean pressure via

const double mean_pressure = VectorTools::compute_mean_value( dof_handler, quadrature_over_integration, solution, dim); 

solution.block(1).add(-mean_pressure);

However ``mean_pressure`` turned out to be very small (at least much smaller than the pressure error I observe), so this subtraction did not practically change the results.

Observed behaviour: (The penalty parameters follow the choice in step-74)

• With homogeneous Dirichlet BCs everything behaves as expected (convergence rates ok).

• With inhomogeneous Dirichlet BCs the pressure L2 error does not converge (velocity behaves better only in RT1)
  ======================= RT1-DGQ1 ========================
  cells           uL2                          uH1                      pL2     
      4      2.067e-01    -         3.128e+00    -        2.055e+00    -
     16     5.324e-02  1.96    1.596e+00  0.97   4.447e-01 2.21
     64     1. 307e-02 2.03    7.968e-01  1.00    1.170e-01 1.93
    256    3.256e-03  2.01    3.983e-01  1.00    4.811e-02 1.28
   1024   8.195e-04 1.99    1.997e-01  1.00     2.965e-02 0.70
   4096   2.086e-04 1.97    1.005e-01  0.99     2.049e-02 0.53
  16384  5.409e-05 1.95    5.090e-02  0.98     1.444e-02 0.51
  
  ======================= RT2-DGQ2 ==========================
  cells            uL2                     uH1                     pL2      
      4        3.405e-02    -       6.559e-01    -       4.237e-01     -
     16       4.760e-03  2.84  1.666e-01  1.98  6.974e-02  2.60
     64       7.869e-04  2.60  5.200e-02  1.68  4.539e-02  0.62
    256      2.285e-04  1.78  2.778e-02  0.90  4.671e-02 -0.04
   1024    8.213e-05  1.48  1.923e-02  0.53  3.686e-02  0.34
   4096    2.966e-05  1.47  1.371e-02  0.49  2.707e-02  0.45
  16384   1.061e-05  1.48  9.751e-03  0.49  1.944e-02  0.48
  

If helpful I can post minimal code snippets (the weak form assembly, the boundary projection call, and the lines where I subtract mean_pressure) and the convergence tables/logs. Thank you again for any guidance.

Best regards,
HC

[Wolfgang Bangerth]

On 10/1/25 05:39, HC Zhang wrote:
> *

>
> ======================= RT1-DGQ1 ========================
> cells           uL2                          uH1
>   pL2
>     4      2.067e-01    -         3.128e+00    -        2.055e+00    -
>    16     5.324e-02  1.96    1.596e+00  0.97   4.447e-01 2.21
>    64     1. 307e-02 2.03    7.968e-01  1.00    1.170e-01 1.93
>   256    3.256e-03  2.01    3.983e-01  1.00    4.811e-02 1.28
>  1024   8.195e-04 1.99    1.997e-01  1.00     2.965e-02 0.70
>  4096   2.086e-04 1.97    1.005e-01  0.99     2.049e-02 0.53
> 16384  5.409e-05 1.95    5.090e-02  0.98     1.444e-02 0.51
>
> ======================= RT2-DGQ2 ==========================
> cells            uL2                     uH1                     pL2
>     4        3.405e-02    -       6.559e-01    -       4.237e-01     -
>    16       4.760e-03  2.84  1.666e-01  1.98  6.974e-02  2.60
>    64       7.869e-04  2.60  5.200e-02  1.68  4.539e-02  0.62
>   256      2.285e-04  1.78  2.778e-02  0.90  4.671e-02 -0.04
>  1024    8.213e-05  1.48  1.923e-02  0.53  3.686e-02  0.34
>  4096    2.966e-05  1.47  1.371e-02  0.49  2.707e-02  0.45
> 16384   1.061e-05  1.48  9.751e-03  0.49  1.944e-02  0.48
>
> If helpful I can post minimal code snippets (the weak form assembly, the
> boundary projection call, and the lines where I subtract mean_pressure)
> and the convergence tables/logs. Thank you again for any guidance.

So if I see this right, the solution actually *does* converge, just not
at the expected rates. That's a different situation than "solution does
not converge" because it means that the function that projects the
boundary values must be doing "something" right at least -- it's not
clear whether it's doing everything right, but it's also at least not
completely wrong :-)

So start to poke some more. Can you pick boundary values that are in the
discrete space? In that case, the discrete boundary values are exact (or
should be) and if you then get the correct convergence rates, you know
that the difference in behavior is due to having exact vs projected
boundary values. If you still get wrong convergence rates, you can
suspect that the problem lies elsewhere.

Best
W.