# Defining 2D mesh for jet stability problem

**URL:** <https://community.freefem.org/t/defining-2d-mesh-for-jet-stability-problem/4266>\
**Category:** General Discussion\
**Created:** [May 3, 2026, 7:33am UTC](https://community.freefem.org/t/defining-2d-mesh-for-jet-stability-problem/4266 "2026-05-03T07:33:10Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Amitkaush](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/amitkaush/32/1978_2.png) [@Amitkaush](https://community.freefem.org/u/Amitkaush)\
**Post date:** [May 3, 2026, 7:33am UTC](https://community.freefem.org/t/defining-2d-mesh-for-jet-stability-problem/4266/1 "2026-05-03T07:33:10Z")

</div>

I am trying to validate my 2D local spatial stability code written in Cartesian coordinates system. I am using the results of Spatial stability of single jet to validate my code. Here is the code snippet where I am defining the mesh:load “PETSc-complex”  
include “ffmatlib.idp”  
macro dimension()2//  
include “macro\_ddm.idp”

real cpu = clock();

// =============================================================  
// PARAMETERS AND FE DEFINITIONS  
// =============================================================  
func Pb2 = [P2, P2];  
func Pc = [P2, P2, P2, P2, P2];

// =============================================================  
// DOMAIN  
// =============================================================  
real x0 = -5, x1 = 5;  
real y0 = -3, y1 = 3;

int n = 60, m = 60;

border C1(t = 0, 2_pi) { x = 8_cos(t); y = 8\*sin(t); }  
mesh Th = buildmesh(C1(15000));

// =============================================================  
// INITIAL MESH ADAPTATION  
// =============================================================  
fespace Kh(Th, P1);  
Kh r = sqrt(x^2 + y^2);

real Rj = 1.2;  
real deltam = 1.0;

Kh indicator = 2 \* exp(-((r - Rj)^2) / (deltam^2));

Th = adaptmesh(Th, indicator,  
err = 0.01,  
hmin = 0.002,  
hmax = 0.3,  
iso = false,  
nbvx = 200000);

// =============================================================  
// BASE FLOW PROFILE  
// =============================================================  
real zpos = 1;

real a = 0.59 + 0.09 \* tanh(sqrt(zpos) - 2.9);  
real delta = (39 + 24_zpos + 0.11_zpos^4) / (1000 + zpos^3);

func real Vz(real rr)  
{  
return 0.5 \* (tanh((rr + a) / delta) - tanh((rr - a) / delta));  
}

func real vcoflow(real rr)  
{  
return 0.005 \* (tanh((rr + a) / delta) - tanh((rr - a) / delta));  
}

real a0 = 0.59 + 0.09 \* tanh(sqrt(0.0) - 2.9);  
real delta0 = (39 + 24_0.0 + 0.11_0.0^4) / (1000 + 0.0^3);

func real Vzref(real rr)  
{  
return 0.5 \* (tanh((rr + a0) / delta0) - tanh((rr - a0) / delta0));  
}

func real Vzsim(real rr)  
{  
return (Vz(rr) + vcoflow(rr)) / (Vzref(0.0) + vcoflow(rr));  
}

// =============================================================  
// BASE FIELD ON MESH  
// =============================================================  
fespace Vh(Th, P2);

Vh W, Wx, Wy, Wxx, Wyy;  
W = Vzsim(sqrt(x_x + y_\*y);

 ![Screenshot from 2026-05-03 13-01-04](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/e0fa2388b451906063f923cc3df5e4d923690027.png)

Above is the base velocity profile.

But I am not getting the good convergence for the eigen values with this mesh definition. What changes should I do?