# IncNavierStokesSolver: Pressure field incorrect despite stable velocity field using Taylor-Hood Approximation on Polynomial order

**URL:** https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17
**Category:** Usage Support
**Created:** [July 16, 2026, 7:19pm UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17 "2026-07-16T19:19:46Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![yaso](https://forum.nektar.info/letter_avatar_proxy/v4/letter/y/8e7dd6/32.png) [@yaso](https://forum.nektar.info/u/yaso)
#### Post date: [July 16, 2026, 7:19pm UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/1 "2026-07-16T19:19:46Z")

</div>

Hello Dear Nektar forum, hope my message finds everyone well.

I have been working on a film-cooling hole nektar++ set-up case. I mainly followed the guidelines from Imperial front wing study and the user guide. I’ve run into a puzzling Taylor-Hood pressure issue that I’d appreciate some insight on.

**Case set-u p**

I started with a coarse linear mesh, has 200K nodes and 1million DOF at P=1. This mesh has one macro prism layer and uses NekMesh prism layer splitting to create the prism layers .

Channel U∞ = 10 m/s, with coolant supplied from a plenum at 0.1309 m/s to match a 10 m/s jet velocity at the hole. The geometry has two sharp features that make meshing difficult; since the bottom corner also has high gradients, I added a fillet there to avoid solver instabilities from skewed cel ls.

I use a RANS field to initialise start-up and start the solution from P=1. This calculation continued for 10 through-flows based on channel length (60D/Uinf; where D=0.00775 m and Uinf=10ms), just to make sure I have a stable starting point when increasing the polynomial order. (P= 1, CFL= 0.04, for testing purposes)

**The issue**

I then used part of the IFW condition xml file as a template to increase the polynomial order to P=2(NUMMODES=3) for velocity and leave the polynomial order at P=1(NUMMODES=2) for pressure, using a Taylor-Hood polynomial expansion. I tested 2 different variants of this.  
Explicit form: (as shown in the my example following IFW below)

```auto
<EXPANSIONS>
<!-- C[0] = Prism C[1] = Tetra -->

<E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMMODES="3,3,3" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="5,5,4" FIELDS="u,v,w" />

<E COMPOSITE="C[1]" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMMODES="3,3,3" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0" NUMPOINTS="5,4,4" FIELDS="u,v,w" />

<E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMMODES="2,2,2" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="5,5,4" FIELDS="p" />

<E COMPOSITE="C[1]" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMMODES="2,2,2" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0" NUMPOINTS="5,4,4" FIELDS="p" />

</EXPANSIONS>

```

Secondly, I tested a compact form following the user guide as below;

```auto
    <E COMPOSITE="C[0,1]" NUMMODES="3" TYPE="MODIFIED" FIELDS="u,v,w" />
    <E COMPOSITE="C[0,1]" NUMMODES="2" TYPE="MODIFIEDQUADPLUS1" FIELDS="p" />

```

I attached the full version of my condition xml file.

The solver can run with a stable CFL throughout the calculation. However unfortunately, I cannot obtain a meaningful pressure field (unphysical values) although the velocity field is showing correct time evolution.

 ![Velocity field (left), pressure field (right)](https://forum.nektar.info/uploads/default/original/1X/c0047d727891f08c2e6508aea6256f9f4b294e40.png)

dt=10-7, CFL=0.04 (for testing purposes)

Example log output:  
 Writing: “cad3p\_p12\_dt107\_6599.chk” (0.105439s, XML)  
CPU time for filter Checkpoint: 0.106312s  
Percentage of time integration: 44.765 %  
Total filters CPU Time: 0.10631 2s  
CG iterations made = 13 using tolerance of 1e-09 (error = 3.52241e-10, rhs\_mag = 65.06 82)  
CG iterations made = 51 using tolerance of 1e-09 (error = 9.66572e-10, rhs\_mag = 2.40845e +08)  
CG iterations made = 45 using tolerance of 1e-09 (error = 9.95594e-10, rhs\_mag = 2.28485 e+08)  
CG iterations made = 49 using tolerance of 1e-09 (error = 9.5522e-10, rhs\_mag = 2.1676 1e+08)  
Steps: 65995 Time: 0.0065995 CPU Time: 0. 230819s  
CFL: 0.0412955 (in elm t 77680)  
CG iterations made = 13 using tolerance of 1e-09 (error = 3.52414e-10, rhs\_mag = 65.0671)  
CG iterations made = 51 using tolerance of 1e-09 (error = 9.66698e-10, rhs\_mag = 2. 40845e+08)  
CG iterations made = 45 using tolerance of 1e-09 (error = 9.9577e-10, rhs\_mag = 2 .28485e+08)  
CG iterations made = 49 using tolerance of 1e-09 (error = 9.54988e-10, rhs\_mag = 2.16761e+08)  
Steps: 65997 Time: 0.0065997 CPU Ti me: 0.233592s

In the case I don’t use Taylor Hood Approximation and set the polynomial order to P=2 for both velocity and pressure, I am able to capture the accurate evolution for both velocity and pressure fields. (P= 2, dt=10-7, CFL= 0.04)

Another detail is the DOF calculated using FieldConvert utility, which is 3million for both cases that using Taylor-Hood appx. and P=2 for both pressure and velocity. Since the pressure polynomial order is lower with Taylor-Hood appx., I would have expected the total DOF to be lower accordingly.

I had several tests to overcome the issue such as initialising pressure to 0 at the start-up(s), testing different tolerances for pressure iterative solver or testing Taylor Hood appx with simplified case; a channel flow domain where the jet is defined as a BC at the wall to isolate the effect of the mesh difficulties, however these also resulted s imilarly.

- I would be very happy to have a bit more understanding of the Taylor-Hood approximation, and how I should be calculating number of points for different element types in case there are any mistakes in my de finitions.

- Are there any specific pressure boundary conditions should I be using at the walls for this f ormulation?

- I was also wondering whether the pressure is obtained through an iterative corrector loop within each timestep (as in PISO-type schemes) or reconstructed directly from the predicted velocity field via a single Poisson solve per step (fractional-step) to get a better understanding on velocity pressure coupling and how sensitive that reconstruction is to initialising pressure to 0 at start-up.

Any pointers or advice would be highly appreciated and thank you very much for your times and kind support.

> **Summary**
>
> This text will be hidden

\<?xml version="1.0" encoding="utf-8" ?\> 

```auto
<E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMMODES="3,3,3" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="5,5,4" FIELDS="u,v,w" />

<E COMPOSITE="C[1]" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMMODES="3,3,3" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0" NUMPOINTS="5,4,4" FIELDS="u,v,w" />

<E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMMODES="2,2,2" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="5,5,4" FIELDS="p" />

<E COMPOSITE="C[1]" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMMODES="2,2,2" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0" NUMPOINTS="5,4,4" FIELDS="p" />

</EXPANSIONS>

<CONDITIONS>

    <SOLVERINFO>
        <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
        <I PROPERTY="EqType" VALUE="UnsteadyNavierStokes" />
        <I PROPERTY="AdvectionForm" VALUE="Convective" />
        <I PROPERTY="Projection" VALUE="Galerkin" />
        <I PROPERTY="SpectralVanishingViscosity" VALUE="True" />
        <I PROPERTY="SPECTRALHPDEALIASING" VALUE="True" />
    </SOLVERINFO>

    <GLOBALSYSSOLNINFO>

        <V VAR="u,v,w">
            <I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
            <I PROPERTY="Preconditioner" VALUE="LowEnergyBlock" />
        <!-- <I PROPERTY="IterativeSolverTolerance" VALUE="1e-9" />
            <I PROPERTY="NekLinSysMaxIterations" VALUE="20000" /> -->
        </V>
        <V VAR="p">
            <I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
            <I PROPERTY="Preconditioner" VALUE="FullLinearSpaceWithLowEnergyBlock" />
        <!-- <I PROPERTY="IterativeSolverTolerance" VALUE="1e-9" />
            <I PROPERTY="NekLinSysMaxIterations" VALUE="20000" /> -->
        </V>
    </GLOBALSYSSOLNINFO>

    <TIMEINTEGRATIONSCHEME>
        <METHOD> IMEX </METHOD>
        <ORDER> 1 </ORDER>
    </TIMEINTEGRATIONSCHEME>

    <PARAMETERS>
        <P> TimeStep = 0.0000001 </P>
        <P> NumSteps = 150000 </P>
        <P> IO_InfoSteps = 1 </P>
        <P> IO_CFLSteps = 1 </P>
        <P> Re = 5100 </P>
        <P> Kinvis = 0.0000149 </P>
        <P> SVVCutoffRatio = 0.3 </P>
        <P> SVVDiffCoeff = 1.0 </P>
    </PARAMETERS>

    <VARIABLES>
        <V ID="0"> u </V>
        <V ID="1"> v </V>
        <V ID="2"> w </V>
        <V ID="3"> p </V>
    </VARIABLES>

    <BOUNDARYREGIONS>
        <B ID="0"> C[11] </B> <!-- coolant inlet -->
        <B ID="1"> C[3] </B> <!-- main inlet -->
        <B ID="2"> C[5] </B> <!-- main outlet -->
        <B ID="3"> C[7] </B> <!-- periodic1 -->
        <B ID="5"> C[9] </B> <!-- periodic2 -->
        <B ID="7"> C[2,12,14,15] </B> <!-- walls -->
    </BOUNDARYREGIONS>

    <BOUNDARYCONDITIONS>
        <REGION REF="0"> <!-- coolant inlet: T_c = 196 K -->
            <D VAR="u" VALUE="0" />
            <D VAR="v" VALUE="0" />
            <D VAR="w" VALUE="0.1309" />
            <N VAR="p" VALUE="0" />
        </REGION>
        <REGION REF="1"> <!-- main inlet: T_inf = 295 K -->
            <D VAR="u" VALUE="10" />
            <D VAR="v" VALUE="0" />
            <D VAR="w" VALUE="0" />
            <N VAR="p" VALUE="0" />
        </REGION>
        <REGION REF="2"> <!-- main outlet -->
            <N VAR="u" VALUE="0" />
            <N VAR="v" VALUE="0" />
            <N VAR="w" VALUE="0" />
            <D VAR="p" VALUE="0" />
        </REGION>
        <REGION REF="3"> <!-- periodic1 -->
            <P VAR="u" VALUE="[5]" />
            <P VAR="v" VALUE="[5]" />
            <P VAR="w" VALUE="[5]" />
            <P VAR="p" VALUE="[5]" />
        </REGION>
        <REGION REF="5"> <!-- periodic2 -->
            <P VAR="u" VALUE="[3]" />
            <P VAR="v" VALUE="[3]" />
            <P VAR="w" VALUE="[3]" />
            <P VAR="p" VALUE="[3]" />
        </REGION>
        <REGION REF="7"> <!-- adiabatic walls: no-slip, zero flux -->
            <D VAR="u" VALUE="0" />
            <D VAR="v" VALUE="0" />
            <D VAR="w" VALUE="0" />
            <N VAR="p" VALUE="0" />
        </REGION>
    </BOUNDARYCONDITIONS>

    <FUNCTION NAME="InitialConditions">
        <F VAR="u,v,w,p" FILE="cad3p_dt3_106_2_77.chk" />
            </FUNCTION>
</CONDITIONS>

<FILTERS>
    <FILTER TYPE="Checkpoint">
        <PARAM NAME="OutputFile">cad3p_p12_dt107</PARAM>
        <PARAM NAME="OutputFrequency">10</PARAM>
    </FILTER>
</FILTERS>

```

---

<div class="post-metadata">

### Author: ![hwustenb](https://forum.nektar.info/letter_avatar_proxy/v4/letter/h/71c47a/32.png) [@hwustenb](https://forum.nektar.info/u/hwustenb)
#### Post date: [July 24, 2026, 10:39am UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/2 "2026-07-24T10:39:23Z")

</div>

Hi yaso,

many thanks for the detailed setup.

I have a few ideas what might be the issue:

- The velocity correction scheme requires a consistent pressure boundary condition for higher-order time accuracy, see section 3.4.6.2 of the User-guide. Set all pressure Neumann condition to use ‘USERDEFINEDTYPE=“H”’. This might already solve the issue.
- The preconditioner `FullLinearSpaceWithLowEnergyBlock` is not heavily used and hence tested. Consider using a diagonal preconditioner for testing.
- The Taylor-Hood and equal-order approximations can become unstable at very small time step sizes consider using a time step aiming for a CFL estimate of CFL ~= 0.3

A few more ideas that are not directly related to your issue:

- Did you try a higher polynomial order that is NUMMODES \> 3?
- If you want to try larger time step sizes, you could also use the implicit velocity correction scheme by replacing `VelocityCorrectionScheme` with `VCSImpliict` and changing the iterative solver to `GMRES`.

Let me know if this helps.

Best wishes,  
Henrik

---

<div class="post-metadata">

### Author: ![alexander.schukmann](https://forum.nektar.info/letter_avatar_proxy/v4/letter/a/c4cdca/32.png) [@alexander.schukmann](https://forum.nektar.info/u/alexander.schukmann)
#### Post date: [October 2, 2026, 7:58am UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/3 "2026-10-02T07:58:35Z")

</div>

I had a similar issue and found it to be caused by the checkpoint filter when used in combination with the T-H approx. Switching to IO\_CheckSteps solved it for me with version 5.7.0.

Edit: There is a more recent post mentioning this issue  
[https://forum.nektar.info/t/checkpoint-filter-fails-with-mixed-expansion-order/21?](https://forum.nektar.info/t/checkpoint-filter-fails-with-mixed-expansion-order/21)

---

<div class="post-metadata">

### Author: ![gaoak](https://forum.nektar.info/letter_avatar_proxy/v4/letter/g/94ad74/32.png) [@gaoak](https://forum.nektar.info/u/gaoak)
#### Post date: [October 3, 2026, 4:58am UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/4 "2026-10-03T04:58:06Z")

</div>

Yes, the solver is all right. Only the output pressure field is broken. I asked gpt, and it said the checkpoint filter does not set the expansion information correctly.

---

<div class="post-metadata">

### Author: ![yaso](https://forum.nektar.info/letter_avatar_proxy/v4/letter/y/8e7dd6/32.png) [@yaso](https://forum.nektar.info/u/yaso)
#### Post date: [October 6, 2026, 1:36pm UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/5 "2026-10-06T13:36:28Z")

</div>

Dear All,

This is very helpful, thank you very much for your time and replies on the issue.

@hwustenb thank you very much also for the recommendations, these were also very helpful for improving the set-up!

I was wondering if the issue was fixed on the new version, nektar++ v5.10?

Thank you very much for your time and support once again

---

<div class="post-metadata">

### Author: ![hwustenb](https://forum.nektar.info/letter_avatar_proxy/v4/letter/h/71c47a/32.png) [@hwustenb](https://forum.nektar.info/u/hwustenb)
#### Post date: [October 7, 2026, 7:03am UTC](https://forum.nektar.info/t/incnavierstokessolver-pressure-field-incorrect-despite-stable-velocity-field-using-taylor-hood-approximation-on-polynomial-order/17/6 "2026-10-07T07:03:19Z")

</div>

Hi,

we discovered an issue with the checkpointing for Taylor-Hood expansions. There is a merge request open for fixing this on the [gitlab](https://gitlab.nektar.info/nektar/nektar/-/merge_requests/2792). Potentially this is the issue both of you encountered.

Best wishes,  
Henrik