# Navier-Stokes Equation

**URL:** <https://getfem.discourse.group/t/navier-stokes-equation/97>\
**Category:** General\
**Created:** [July 28, 2025, 12:57pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97 "2025-07-28T12:57:30Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Giovanni.Gelmetti](https://avatars.discourse-cdn.com/v4/letter/g/a698b9/32.png) [@Giovanni.Gelmetti](https://getfem.discourse.group/u/Giovanni.Gelmetti)\
**Post date:** [July 28, 2025, 12:57pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/1 "2025-07-28T12:57:30Z")

</div>

Hello, my goal is to solve the full Navier-Stokes equations for a beam immersed in a fluid (air). The inlet velocity profile is parabolic, with a ramping factor applied to ensure a smoother transition at the start of the simulation. However, after approximately 100 iterations, the solution stops converging. Could this be caused by how I’ve configured the time integration method or the boundary conditions?

```python
import getfem as gf 
from regions_check import verify_regions
import numpy as np
import os

gf.util_trace_level(1)
gf.util_warning_level(1)
π = np.pi
################################
# Beam in a Fluid:
#################################
'''
The goal is to solve the full Navier-Stokes equations for a beam immeresed in a fluid.
'''

##################
## PROBLEM DATA ##
##################
"""
The fluid it will be the air
"""

scale_factor = 1 #m -> m (it can be used if we want to pass from meter to anything else cm, dm etc...) 

# Fluid properties (air):
ν_fluid = 0.015e-4 * scale_factor**2 # m²/s 
rho_fluid = 1.23 /(scale_factor**3 ) # kg/m³ 

# Boundaries values: 
# Inlet velocity parameters
U_mean = 1.0 # Mean inlet velocity
H = 0.41 # Channel height

# Transient paramters: 
T = 10.0 # Total simulation time
dt = 1e-3 # Time step
theta = 0.5 # Theta parameter (0.5 = Crank-Nicolson)

#############
## MESH ##
#############

Mesh_fluid= gf.Mesh('Import', 'gmsh','MESH_GMSH/Fluid.msh')

#############
## REGIONS ##
#############
"""
regions flagging in gmsh does not work thus it's necessary to give a region for each line or physical surface. 
 
"""

Beam_left = 1 
Beam_bottom = 2 
Beam_right = 3 
Beam_top = 4 
INLET = 5 
Wall_bottom = 6
OUTLET = 7  
Wall_top = 8 

BEAM_INTERFACE = 100

Mesh_fluid.region_merge(BEAM_INTERFACE, Beam_left)
Mesh_fluid.region_merge(BEAM_INTERFACE, Beam_right)
Mesh_fluid.region_merge(BEAM_INTERFACE, Beam_bottom)
Mesh_fluid.region_merge(BEAM_INTERFACE, Beam_top)

# physical surface: 
fluid1 = 9
fluid2 = 11
fluid3 = 10
fluid4 = 13

FLUID = 101
Mesh_fluid.region_merge(FLUID,fluid1)
Mesh_fluid.region_merge(FLUID,fluid2)
Mesh_fluid.region_merge(FLUID,fluid3)
Mesh_fluid.region_merge(FLUID,fluid4)

WALLS = 102
Mesh_fluid.region_merge(WALLS,Wall_top)
Mesh_fluid.region_merge(WALLS,Wall_bottom)

verify_regions(Mesh_fluid, 'meshfluid') #it cna be used verify regions againg to see if the regions are set up correctly.

#verify_regions(Mesh_solid, 'meshsolid')
########################
## INTEGRATION METHOD ##
########################
"""
The integration method is quadrature with 17 points, which is suitable for QK elements.  
"""
mim_fluid = gf.MeshIm(Mesh_fluid, gf.Integ('IM_QUAD(17)'))

#########################
## FEM ELEMENTS METHOD ##
########################
"""
THe meshes are quadrilater, hence QK elements are used. 
"""

## FEM ELEMENTS:

mfv_fluid_ = gf.MeshFem(Mesh_fluid, 2)
mfv_fluid_.set_fem(gf.Fem('FEM_QK(2,2)'))

mfp_fluid = gf.MeshFem(Mesh_fluid, 1)
mfp_fluid.set_fem(gf.Fem('FEM_QK(2,1)'))

######################
## REMOVE FIXED DOF ##
######################
"""
Instead of applying dirichelt boundary coniditons, the degree of freedom are removed directly here.
"""

kept_dofs_v_f = list(
                set(range(mfv_fluid_.nbdof()))
                -set(mfv_fluid_.basic_dof_on_region(WALLS))
                -set(mfv_fluid_.basic_dof_on_region(BEAM_INTERFACE))
                )
                

mfv_fluid = gf.MeshFem('partial', mfv_fluid_, kept_dofs_v_f)

###########
## MODEL ##
###########

md = gf.Model("real")

###################
## FEM VARIABLES ##
###################
"""
Fem variable are defined on the MeshFEM with removed degree of freedom. 

"""

md.add_fem_variable("v_f", mfv_fluid) # Fluid velocity  
md.add_fem_variable("p", mfp_fluid) # Pressure

# Previous time step data
md.add_fem_data("Previous_v_f", mfv_fluid)

###########################
## INITIALIZED CONSTANTS ##
###########################
"""
The material properties and the times variable are initiualized in the model 
"""
md.add_initialized_data("rho_f", rho_fluid)
md.add_initialized_data("nu_f", ν_fluid)
md.add_initialized_data("dt", dt)
md.add_initialized_data("theta", theta)
md.add_initialized_data("H", H)
md.add_initialized_data("U_mean", U_mean)
md.add_initialized_data("p_out", 0)
# Time variable
md.add_initialized_data("t", 0.0)
# Define inlet profile as a macro (parabolic profile)
md.add_macro("inlet_profile(y)", "1.5*U_mean*4*y*(H-y)/(H*H)")

########################
## INITIAL CONDITIONS ##
########################
"""
Initialize all to zero (at rest), since we're using filtered fem variable, it's not possibble to 
use md.interpolation. 
"""

# Initialize with zero vectors of the correct size
md.set_variable("v_f", np.zeros(mfv_fluid.nbdof()))
md.set_variable("p", np.zeros(mfp_fluid.nbdof()))

# Also set previous values
md.set_variable("Previous_v_f", np.zeros(mfv_fluid.nbdof()))

######################
## WEAK FORMULATION ##
# ######################
"""
The problem is formulated in weak form, and for the time the theta method is used.
"""

# FLUID : 
# stress tensors: 
md.add_macro('sigma_f_vu(v)', "2*rho_f*nu_f*(Grad_v + (Grad_v)) ")
md.add_macro('sigma_f_p(p)', '-p*Id(2)')

md.add_macro("Convection(v)", "(rho_f*v.Grad_v)")
md.add_macro('Incompressibility(v)', "Trace(Grad_v)")
md.add_macro("Stress_vu(v)", 'sigma_f_vu(v)')
md.add_macro("Stress_p(p)", "sigma_f_p(p)")

Transient_fluid = '(rho_f*theta*(v_f - Previous_v_f)). Test_v_f' \
            '+ (rho_f*(1-theta)*(v_f - Previous_v_f)). Test_v_f '

md.add_nonlinear_term(mim_fluid, Transient_fluid, FLUID)

md.add_linear_term(mim_fluid, ' dt*Incompressibility(v_f)*Test_p', FLUID) 

md.add_nonlinear_term(mim_fluid,"dt*theta*(Convection(v_f)).Test_v_f +" \
                                "dt*(1-theta)*(Convection(Previous_v_f)).Test_v_f", FLUID)

md.add_linear_term(mim_fluid,"dt*theta*Stress_vu(v_f):Grad_Test_v_f +" \
                                "dt*(1-theta)*(Stress_vu(Previous_v_f):Grad_Test_v_f)", FLUID)
md.add_linear_term(mim_fluid, 'dt*Stress_p(p):Grad_Test_v_f', FLUID)

#####################
# BOUNDARY CONDITIONS
# Inlet velocity profile (parabolic, time-dependent with smooth ramp)
ramp_factor = 0.0
V_inlet_expr = f"{ramp_factor}*[inlet_profile(X(2)), 0]"
V_inlet_full = md.interpolation(V_inlet_expr, mfv_fluid_)
V_inlet = V_inlet_full[kept_dofs_v_f]

md.add_initialized_fem_data('V_inlet', mfv_fluid, V_inlet) 

# Replace penalty terms with:
md.add_Dirichlet_condition_with_penalization(mim_fluid, "v_f", 10**3, INLET, "V_inlet")
md.add_Dirichlet_condition_with_penalization(mim_fluid, 'p', 10**3, OUTLET)

####################
## MODEL SOLOTION ##
####################

# Create output directory
output_dir = "Results"
os.makedirs(output_dir, exist_ok=True)

# TIME STEPPING LOOP

print(f"Time step: {dt}, Total time: {T}")
print(f"Theta: {theta} (Crank-Nicolson)" if theta == 0.5 else f"Theta: {theta}")

t = 0.0
step = 0
export_every = int(0.01/dt) # Export every 0.01s

# Main time loop
while t < T:
    
    md.set_variable("t", t)
    if t < 2.0:
        ramp_factor = 0.5 * (1 - np.cos(np.pi * t / 2.0))
    else:
        ramp_factor = 1.0
    V_inlet_expr = f"{ramp_factor}*[inlet_profile(X(2)), 0]"
    V_inlet_full = md.interpolation(V_inlet_expr, mfv_fluid_)
    V_inlet = V_inlet_full[kept_dofs_v_f]
    md.set_variable("V_inlet", V_inlet)
    if step > 0: # After first solve
        
        print(f'time is: {t} and ramp_factor is {ramp_factor} and velocity is {V_inlet}')
    
    
    # More robust solver parameters
    md.solve("very noisy", 
         "max_iter", 100,
         "max_res", 1e-6,  
         "lsolver", "superlu",
         "lsearch", "basic",  
         "alpha min", 1e-4,  
         "alpha mult", 0.5)    

    # Extract current solution

    v_f = md.variable("v_f")
    p = md.variable("p")
   
    

    mfv_fluid.export_to_vtu(f"{output_dir}/fluid_{step:05d}.vtu",
                            mfv_fluid, v_f, "Velocity", 
                            mfp_fluid, p, "Pressure")
    
    
    md.set_variable("Previous_v_f", v_f)
    
    
    # Advance time
    t += dt
    step += 1
    

```

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [July 29, 2025, 11:54am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/2 "2025-07-29T11:54:01Z")

</div>

Hello, nice project, something weird that I have spotted is this: `IM_QUAD(17)`

 ![image](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/522b70c8c807e8a3dc1933c475007e4ffbaf0b88.png)

A 17th order integration looks a bit extreme.

Otherwise I lack the mesh file for running your code.

---

<div class="post-metadata">

**Author:** ![Giovanni.Gelmetti](https://avatars.discourse-cdn.com/v4/letter/g/a698b9/32.png) [@Giovanni.Gelmetti](https://getfem.discourse.group/u/Giovanni.Gelmetti)\
**Post date:** [July 30, 2025, 8:35am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/3 "2025-07-30T08:35:13Z")

</div>

Hello,

thank you for the answer.

I know that the integration order is high I was trying to see if something would change if I increased it and I forgot it to change it back, but it doesn’t seem to affect the solution.

I’m not sure how to share the .msh file with you, but this is the .geo file. Could I also ask you which is, in your opinion, the best way to impose boundary conditions and which solver to set in case of a non-linear problem?"

```c

////////////////////////////////////////////////////////////
// Constants
h_ref = 0.5;
L = 10;
thickness = 0.005*L;

// Number of elements along each direction
nx1 = 20; // Left of beam
nx2 = 40; // Along beam
nx3 = 20; // Right of beam
ny1 = 20; // Below/Above beam (symmetric)
ny2 = 5; // Along beam thickness

// Fluid domain dimensions (made symmetric)
domain_x_min = -2*L;     
domain_x_max = 3*L;      
domain_y_min = -2*L; // Made symmetric with domain_y_max
domain_y_max = 2*L;      

////////////////////////////////////////////////////////////
// Points Definition
////////////////////////////////////////////////////////////
// Outer corners
Point(1) = {domain_x_min, domain_y_min, 0, h_ref}; // Bottom left
Point(2) = {domain_x_max, domain_y_min, 0, h_ref}; // Bottom right
Point(3) = {domain_x_max, domain_y_max, 0, h_ref}; // Top right
Point(4) = {domain_x_min, domain_y_max, 0, h_ref}; // Top left

// Points for beam
Point(5) = {0, -thickness/2, 0, h_ref}; // Beam bottom left
Point(6) = {L, -thickness/2, 0, h_ref}; // Beam bottom right
Point(7) = {L, thickness/2, 0, h_ref}; // Beam top right
Point(8) = {0, thickness/2, 0, h_ref}; // Beam top left

//+
Line(1) = {8, 5};
//+
Line(2) = {5, 6};
//+
Line(3) = {6, 7};
//+
Line(4) = {7, 8};

//+
Line(5) = {4, 1};
//+
Line(6) = {1, 2};
//+
Line(7) = {2, 3};
//+
Line(8) = {3, 4};
//+

// Meshing lines
//+
Line(9) = {8, 4};
//+
Line(10) = {5, 1};
//+
Line(11) = {7, 3};
//+
Line(12) = {6, 2};

Curve Loop(1) = {5, 6, 7, 8};
//+
Curve Loop(2) = {2, 3, 4, 1};

//+
Physical Curve("beam_left", 9) = {1};
//+
Physical Curve("beam_bottom", 10) = {2};
//+
Physical Curve("beam_right", 11) = {3};
//+
Physical Curve("beam_top", 12) = {4};
//+
Physical Curve("Inlet", 14) = {5};
//+
Physical Curve("Outlet", 15) = {7};
//+
Physical Curve("Walls", 16) = {8, 6};
//+

//+
Curve Loop(3) = {5, -10, -1, 9};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {4, 9, -8, -11};
//+
Plane Surface(4) = {4};
//+
Curve Loop(5) = {2, 12, -6, -10};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {12, 7, -11, -3};
//+
Plane Surface(6) = {6};
//+

// Physical Surfaces
Physical Surface("fluid", 18) = {3, 4, 5, 6};

//+
Transfinite Surface {3} = {1, 5, 8, 4};
//+
Transfinite Surface {5} = {1, 2, 6, 5};
//+
Transfinite Surface {6} = {6, 2, 3, 7};
//+
Transfinite Surface {4} = {3, 4, 8, 7};
//+
Transfinite Surface {2} = {8, 5, 6, 7};

//+
N_v = 20;
N_o = 60;
N= 80;
p = 1.05;
Transfinite Curve {5, 1} = N_v Using Progression 1;
//+
Transfinite Curve {9, 10} = N Using Progression p;

//+
Transfinite Curve {6, 2} = N_o Using Progression 1;
//+
Transfinite Curve {10, 12} = N Using Progression p;

//+
Transfinite Curve {8, 4} = N_o Using Progression 1;
//+
Transfinite Curve {11, 9} = N Using Progression p;

//+
Transfinite Curve {7, 3} = N_v Using Progression 1;
//+
Transfinite Curve {12, 11} = N Using Progression p;

Recombine Surface{3,4,5,6}; 
Mesh.Algorithm = 6;

Mesh 2;
Save "Fluid.msh";

```

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [July 30, 2025, 9:24am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/4 "2025-07-30T09:24:00Z")

</div>

already in the iteration that converge you see a checkerboard pattern in the velocity field

 ![image](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/69f2a384886559a4e77bedf76f07fbbb7c18bb32.jpeg)

the discretization of your v-p fields does not seem to be stable

I have also noticed this

```auto
(Grad_v + (Grad_v))

```

do you mean

```auto
(Grad(v) + Grad(v)')

```

?

One can also us the sym operator

```auto
2*Sym(Grad(v))

```

In `md.solve`, do not use “very noisy“ it is just annoying. Just use “noisy“. I also prefer to use mumps instead of superlu. If your finite element spaces and equations are ok, then you will not really need linesearch. In general, linesearch for mixed problems is ugly. Unless I model contact, or some other semismooth problem, I prefer to write a standard undamped Newton loop.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [July 30, 2025, 9:50am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/5 "2025-07-30T09:50:38Z")

</div>

This here I do not understand either

> [@Giovanni.Gelmetti](#):
>
> ```auto
> Transient_fluid = '(rho_f*theta*(v_f - Previous_v_f)). Test_v_f' \
> '+ (rho_f*(1-theta)*(v_f - Previous_v_f)). Test_v_f '
> 
> ```

theta will cancel out in the sum of these two subexpressions

---

<div class="post-metadata">

**Author:** ![Giovanni.Gelmetti](https://avatars.discourse-cdn.com/v4/letter/g/a698b9/32.png) [@Giovanni.Gelmetti](https://getfem.discourse.group/u/Giovanni.Gelmetti)\
**Post date:** [October 20, 2025, 7:39am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/6 "2025-10-20T07:39:38Z")

</div>

Hello, Thanks for the advice! I’ve fixed the errors and modified the problem to simulate a cylinder immersed in fluid at a lower Reynolds number. However, I’m running into some issues: 1. The solution appears to be laminar with no vortex shedding occurring 2. The velocity field should be divergence-free, but I’m consistently getting divergence values around 0.05. Here’s my current implementation:

```python
import getfem as gf 
import numpy as np
import os

############################{####
# Cylinder in a Fluid:
#################################
'''
The goal is to solve the full Navier-Stokes equations for a cylinder immeresed in a fluid.
'''

##################
## PROBLEM DATA ##
##################

# Fluid properties (artificial values for a Reynolds number of 60 ):
ν_fluid = 1/60 
rho_fluid = 1.204 # kg/m³ 
mu_fluid = rho_fluid * ν_fluid # Dynamic viscosity kg/(m·s)

# Boundaries values: 
# Inlet velocity parameters
U_mean = 1 # Mean inlet velocity m/s

# Reference values for coefficients
D = 1 # Cylinder diameter (you need to set this based on your mesh)
A_ref = D # Reference area (diameter × unit depth for 2D)
q_inf = 0.5 * rho_fluid * U_mean**2 # Dynamic pressure

# Transient paramters: 
T = 5.0 # Total simulation time
dt = 0.002 # Time step
theta = 0.5 # Theta parameter (0.5 = Crank-Nicolson)

print(f"Reynolds number is: {U_mean*D/ν_fluid} and number of time steps: {int(T/dt)}")
print(f"Strouhal number is approximately: {0.2}, thus a peiod is T = {1/0.2*D/U_mean}")

############
## MESH ##
#############

Mesh_fluid= gf.Mesh('Import', 'gmsh','Mesh/cylinder.msh')
h = min(Mesh_fluid.convex_radius())
print( f"Minimum mesh size h ={h}, and CFL = "f"{1}, thus dt = {dt} should be less than {1*h/U_mean}" )

#############
## REGIONS ##
#############
"""
regions flagging in gmsh does not work thus it's necessary to give a region for each line or physical surface. 

"""

Bottom_left = 1 
Bottom_right = 2
Top_left = 3
Top_right = 4 
INLET = 5 
OUTLET = 7  
Cylinder_1 = 8
Cylinder_2 = 9
Cylinder_3 = 10
Cylinder_4 = 11
print("Regions in the mesh are:", Mesh_fluid.regions())

CYLINDER_INTERFACE = 2001

Mesh_fluid.region_merge(CYLINDER_INTERFACE, Cylinder_1)
Mesh_fluid.region_merge(CYLINDER_INTERFACE, Cylinder_2)
Mesh_fluid.region_merge(CYLINDER_INTERFACE, Cylinder_3)
Mesh_fluid.region_merge(CYLINDER_INTERFACE, Cylinder_4)

# physical surface: 
fluid1 = 6 
fluid2 = 12
fluid3 = 13
fluid4 = 14
fluid5 = 15

FLUID = 2002
Mesh_fluid.region_merge(FLUID,fluid1)
Mesh_fluid.region_merge(FLUID,fluid2)
Mesh_fluid.region_merge(FLUID,fluid3)
Mesh_fluid.region_merge(FLUID,fluid4)
Mesh_fluid.region_merge(FLUID,fluid5)

WALLS = 2003
Mesh_fluid.region_merge(WALLS,Bottom_left)
Mesh_fluid.region_merge(WALLS,Bottom_right)
Mesh_fluid.region_merge(WALLS,Top_left)
Mesh_fluid.region_merge(WALLS,Top_right)

########################
## INTEGRATION METHOD ##
########################
"""
The integration method is quadrature with 5 points, which is suitable for QK elements.  
"""
mim_fluid = gf.MeshIm(Mesh_fluid, gf.Integ('IM_QUAD(5)'))

#########################
## FEM ELEMENTS METHOD ##
########################1
"""
THe meshes are quadrilater, hence QK elements are used. 
"""

## FEM ELEMENTS:

mfv_fluid = gf.MeshFem(Mesh_fluid, 2)
mfv_fluid.set_fem(gf.Fem('FEM_QK(2,2)'))

mfp_fluid = gf.MeshFem(Mesh_fluid, 1)
mfp_fluid.set_fem(gf.Fem('FEM_QK(2,1)'))

###########
## MODEL ##
###########

md = gf.Model("real")

###################
## FEM VARIABLES ##
###################
"""
Fem variable are defined on the MeshFEM with removed degree of freedom. 

"""

md.add_fem_variable("v_f", mfv_fluid) # Fluid velocity  
md.add_fem_variable("p", mfp_fluid) # Pressure

# Previous time step data
md.add_fem_data("Previous_v_f", mfv_fluid)

###########################
## INITIALIZED CONSTANTS ##
###########################
"""
The material properties and the times variable are initiualized in the model 
"""
md.add_initialized_data("rho_f", rho_fluid)
md.add_initialized_data("nu_f", ν_fluid)
md.add_initialized_data("dt", dt)
md.add_initialized_data("theta", theta)
md.add_initialized_data("U_mean", U_mean)

# Time variable
md.add_initialized_data("t", 0.0)
# Define inlet profile as a macro (parabolic profile)

########################
## INITIAL CONDITIONS ##
########################
"""
Initialize all to zero (at rest), since we're using filtered fem variable, it's not possibble to 
use md.interpolation. 
"""

# Initialize with zero vectors of the correct size
md.set_variable("v_f", np.zeros(mfv_fluid.nbdof()))
md.set_variable("p", np.zeros(mfp_fluid.nbdof()))

# Also set previous values
md.set_variable("Previous_v_f", np.zeros(mfv_fluid.nbdof()))

######################

## WEAK FORMULATION ##
# ######################
"""
The problem is formulated in weak form and for the time the theta method is used.
"""

# FLUID : 
# stress tensors: 
md.add_macro('Stress_vu(v)', "rho_f*nu_f*2*Sym(Grad_v) ")
md.add_macro('Stress_p(p)', '-p*Id(2)')
md.add_macro("Convection(v)", "(rho_f*v.Grad_v)")
md.add_macro('Incompressibility(v)', "Trace(Grad_v)")

Transient_fluid = '(rho_f/dt)*(v_f - Previous_v_f). Test_v_f' 
            
md.add_nonlinear_term(mim_fluid, Transient_fluid, FLUID)

md.add_linear_term(mim_fluid, ' Incompressibility(v_f)*Test_p', FLUID) 

md.add_nonlinear_term(mim_fluid,"theta*(Convection(v_f)).Test_v_f +" \
 "(1-theta)*(Convection(Previous_v_f)).Test_v_f", FLUID)

md.add_linear_term(mim_fluid,"theta*Stress_vu(v_f):Grad_Test_v_f +" \
                                "(1-theta)*(Stress_vu(Previous_v_f):Grad_Test_v_f)", FLUID)
md.add_linear_term(mim_fluid, 'Stress_p(p):Grad_Test_v_f', FLUID)

#####################
# BOUNDARY CONDITIONS
#####################

V_inlet = md.interpolation("[U_mean, 0]", mfv_fluid)
md.add_initialized_fem_data('V_inlet', mfv_fluid, V_inlet) 

V_walls = md.interpolation("[U_mean, 0]", mfv_fluid)
md.add_initialized_fem_data('V_walls', mfv_fluid, V_walls) 

V_cylinder = md.interpolation("[0, 0]", mfv_fluid)
md.add_initialized_fem_data('V_cylinder', mfv_fluid, V_cylinder) 

# Apply Dirichlet conditions with multipliers
md.add_Dirichlet_condition_with_multipliers(mim_fluid, "v_f", mfv_fluid, INLET, "V_inlet")
md.add_Dirichlet_condition_with_multipliers(mim_fluid, "v_f", mfv_fluid, WALLS, "V_inlet")
md.add_Dirichlet_condition_with_multipliers(mim_fluid, "v_f", mfv_fluid, CYLINDER_INTERFACE, "V_cylinder")

####################
## MODEL SOLOTION ##
####################

# Create output directory
output_dir = "Results"
os.makedirs(output_dir, exist_ok=True)

# TIME STEPPING LOOP

print(f"Time step: {dt}, Total time: {T}, number of steps: {int(T/dt)}")

t = 0.0
step = 0

time_history = []
cd_history = []
cl_history = []

# Main time loop
while t < T:
    
  md.set_variable("t", t)
  md.solve("noisy", 
         "max_iter", 100,
         "max_res", 1e-8,  
         "lsolver", "superlu",  
         "alpha min", 1e-4,  
         "alpha mult", 0.5)    

  # Extract current solution
  v_f = md.variable("v_f")
  p = md.variable("p")

  print(f'time is: {t} the velocity is {v_f}')

  
  time_ms = int(t * 1000) # Convert to milliseconds
  mfv_fluid.export_to_vtu(f"{output_dir}/fluid_{time_ms:06d}.vtu",
                        mfv_fluid, v_f, "Velocity", 
                        mfp_fluid, p, "Pressure")
    

  md.set_variable("Previous_v_f", v_f)

  ###### CL and CD computation ######
  traction = gf.asm_generic(mim_fluid, 0, "Stress_vu(v_f)*Normal+ Stress_p(p)*Normal", CYLINDER_INTERFACE, md)
  Fx, Fy = -traction[0], -traction[1]
  Cd = Fx / (q_inf * A_ref)
  Cl = Fy / (q_inf * A_ref)
  print(f"Time: {t:.4f}, Cd: {Cd:.6f}, Cl: {Cl:.6f}")

  time_history.append(t)
  cd_history.append(Cd)
  cl_history.append(Fy)
  np.savetxt(f"force_coefficients.txt", 
        np.column_stack([time_history, cd_history, cl_history]),
        header="Time, Cd, Cl")
  
  ## some checks to see if the solution is resonable ##
  div_norm = np.sqrt(gf.asm_generic(mim_fluid, 0, 'pow((Trace(Grad_v_f)),2)', FLUID, md))
  print("‖div(v_f)‖ₗ₂ =", div_norm)

  # Advance time
  t += dt
  step += 1
  
print("simulation completed.")

```

and here there is the .geo file for the mesh:

```c
// Cylinder Mesh

Nx1 = 61; Rx1 = 1.00;
Nx2 = 41; Rx2 = 1.035;
Ny = 91; Ry = 4.50;
Nb = 121; Rb= 0.96;
Nc = 61; Rc= 1.00;
// Domain points
//+
Point(1) = {-10, -10, 0, 1.0};
//+
Point(2) = {10, -10, 0, 1.0};
//+
Point(3) = {40, -10, 0, 1.0};
//+
Point(4) = {-10, 10, 0, 1.0};
//+
Point(5) = {10, 10, 0, 1.0};
//+
Point(6) = {40, 10, 0, 1.0};
// Cylinder points
// radius r=0.5 so 0.5cos(45) = 0.35355339
//+
Point(7) = {-0.35355339, -0.35355339, 0, 1.0};
//+
Point(8) = {0.35355339, -0.35355339, 0, 1.0};
//+
Point(9) = {-0.35355339, 0.35355339, 0, 1.0};
//+
Point(10) = {0.35355339, 0.35355339, 0, 1.0};
//+
Point(11) = {0, -0, 0, 1.0};
//+

// Domain Lines
Line(1) = {1, 2}; Transfinite Curve {1} = Nx1 Using Progression Rx1;
//+
Line(2) = {2, 3}; Transfinite Curve {2} = Nx2 Using Progression Rx2;
//+
Line(3) = {4, 5}; Transfinite Curve {3} = Nx1 Using Progression Rx1;
//+
Line(4) = {5, 6}; Transfinite Curve {4} = Nx2 Using Progression Rx2;
//+
Line(5) = {1, 4}; Transfinite Curve {5} = Ny Using Bump Ry;
//+
Line(6) = {2, 5}; Transfinite Curve {6} = Ny Using Bump Ry;
//+
Line(7) = {3, 6}; Transfinite Curve {7} = Ny Using Bump Ry;

// Cylinder arc
//+
Circle(8) = {7, 11, 8}; Transfinite Curve {8} = Nc Using Progression Rc;
//+
Circle(9) = {8, 11, 10}; Transfinite Curve {9} = Ny Using Progression Rc;
//+
Circle(10) = {10, 11, 9}; Transfinite Curve {10} = Nc Using Progression Rc;
//+
Circle(11) = {9, 11, 7}; Transfinite Curve {11} = Ny Using Progression Rc;
//+

// diagonals block lines 
Line(12) = {1, 7}; Transfinite Curve {12} = Nb Using Progression Rb;
//+
Line(13) = {2, 8}; Transfinite Curve {13} = Nb Using Progression Rb;
//+
Line(14) = {5, 10}; Transfinite Curve {14} = Nb Using Progression Rb;
//+
Line(15) = {4, 9}; Transfinite Curve {15} = Nb Using Progression Rb;

// Surfaces
//+
Curve Loop(1) = {12, 8, -13, -1};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {13, 9, -14, -6};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {14, 10, -15, 3};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {15, 11, -12, 5};
//+
Plane Surface(4) = {4};
//+
Curve Loop(5) = {7, -4, -6, 2};
//+
Plane Surface(5) = {5};
//+
Transfinite Surface {1};
//+
Transfinite Surface {2};
//+
Transfinite Surface {3};
//+
Transfinite Surface {4};
//+
Transfinite Surface {5};
//+
Recombine Surface {1};
//+
Recombine Surface {2};
//+
Recombine Surface {3};
//+
Recombine Surface {4};
//+
Recombine Surface {5};
//+
Physical Surface("Fluid", 101) = {4, 1, 2, 3, 5};
//+
Physical Curve("Cylinder", 1005) = {8, 9, 10, 11};
//+
Physical Curve("Inlet", 1001) = {5};
//+
Physical Curve("Outlet", 1002) = {7};
//+
Physical Curve("Wall", 1003) = {4, 3, 1, 2};
//+
Save "cylinder.msh";

```

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [October 20, 2025, 9:14am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/7 "2025-10-20T09:14:04Z")

</div>

before trying your code, I would like to ask you to try a modification yourself. Use a well established element formulation for your benchmarks, such as the Taylor-Hood element.  
From gmsh export a mesh with linear triangles. In getfem define a quadratic meshfem for velocities and a linear one for pressure.

```auto
mfv_fluid.set_classical_fem(2) # equivalent to "FEM_PK(2,2)" on triangles
mfp_fluid.set_classical_fem(1) # equivalent to "FEM_PK(2,1)" on triangles

```

---

<div class="post-metadata">

**Author:** ![Giovanni.Gelmetti](https://avatars.discourse-cdn.com/v4/letter/g/a698b9/32.png) [@Giovanni.Gelmetti](https://getfem.discourse.group/u/Giovanni.Gelmetti)\
**Post date:** [October 22, 2025, 1:26pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/8 "2025-10-22T13:26:25Z")

</div>

Thank you! It seems to work better like this. I have a couple more questions:

1. For quadrilateral elements, which combination of elements is considered stable?

2. When I define a macro like this:

3. Finally, if I stop a simulation before it’s finished and want to restart it, is it possible to retrieve the velocity and pressure fields from a `.vtu` file?

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [October 22, 2025, 4:24pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/9 "2025-10-22T16:24:55Z")

</div>

1. Quadrilateral elements for Stokes flow are not as easy. There is an element known as “mini” but I am not sure if GetFEM implements it, because as far as I know we only have bubble shape functions for triangles/tetrahedrals but not for quads. @yves.renard is that correct?

See also:

> **[Matlab source code for the Stokes equation example — GetFEM](https://getfem.org/screenshots/stokes-source.html)**

1. yes, exactly

2. no, vtu is not an input/output format, vtu is suitable only for output. You can just dump the state of the model to a numpy file and reload it afterwards, with

```python
 np.save("model_state.npy", md.from_variables())

```

and

```python
 md.to_variables(np.load("model_state.npy"))

```

where `md` is yout model object.

---

<div class="post-metadata">

**Author:** ![yves.renard](https://avatars.discourse-cdn.com/v4/letter/y/df705f/32.png) [@yves.renard](https://getfem.discourse.group/u/yves.renard)\
**Post date:** [October 22, 2025, 5:00pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/10 "2025-10-22T17:00:30Z")

</div>

Yes, I confirm. The mini element is not available, there is not bubble enriched quadrialteral element in getfem. Adding this element would be rather easy, but it has not been done yet.

---

<div class="post-metadata">

**Author:** ![yves.renard](https://avatars.discourse-cdn.com/v4/letter/y/df705f/32.png) [@yves.renard](https://getfem.discourse.group/u/yves.renard)\
**Post date:** [October 22, 2025, 5:03pm UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/11 "2025-10-22T17:03:11Z")

</div>

However, Q2 x Q1 is also stable I think, and available.

---

<div class="post-metadata">

**Author:** ![Giovanni.Gelmetti](https://avatars.discourse-cdn.com/v4/letter/g/a698b9/32.png) [@Giovanni.Gelmetti](https://getfem.discourse.group/u/Giovanni.Gelmetti)\
**Post date:** [December 3, 2025, 10:27am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/12 "2025-12-03T10:27:02Z")

</div>

hello, i was trying to use the FEM\_PK\_WITH\_CUBIC\_BUBBLE but i get the following error: “Only Lagrange fems can be used in interpolation”, what should i do instead ?

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [December 3, 2025, 11:07am UTC](https://getfem.discourse.group/t/navier-stokes-equation/97/13 "2025-12-03T11:07:38Z")

</div>

well what you actually miss is the syntax for projecting some field between different FE spaces

let’s say mfL is a Lagrange mesh\_fem and mfNL is a non Lagrange one. You can still interpolate on mfL and the project the result onto mfNL with

```auto
V0 = md.interpolation("[U_mean, 0]", mfL)
mass_mat = gf.asm_mass_matrix(mim, mfNL)
Vasm = gf.asm_generic(mim, 1, "v0*Test_t", -1,
                       "v0", False, mfL, V0,
                       "t", True, mfNL, np.zeros(mfNL.nbdof())
                       )
V = gf.linsolve_mumps(mass_mat, Vasm).T

```

by the way, you can also transfer interpolation results between reduced and unreduced mesh\_fems using the methods

`mf.reduce_vector(V)`

and

`mf.extend_vector(V)`

see documentation at

> **[MeshFem — GetFEM](https://getfem.org/python/cmdref_MeshFem.html#getfem.MeshFem.MeshFem)**