# Spatial stability analysis for compressible jet flow
**URL:** https://community.freefem.org/t/spatial-stability-analysis-for-compressible-jet-flow/4151
**Category:** General Discussion
**Created:** [December 11, 2025, 1:57pm UTC](https://community.freefem.org/t/spatial-stability-analysis-for-compressible-jet-flow/4151 "2025-12-11T13:57:01Z")
**Posts on this page:** 1
**Page:** 1
### Author:  [@Kumar](https://community.freefem.org/u/Kumar)
#### Post date: [December 11, 2025, 1:57pm UTC](https://community.freefem.org/t/spatial-stability-analysis-for-compressible-jet-flow/4151/1 "2025-12-11T13:57:01Z")
Hi all,
I am working on spatial stability analysis of 1D jet flow in a polar coordinate system.
I am trying to validate the growth rate plot from this Linear and nonlinear mechanisms of sound
radiation by instability waves in subsonic jets VICTORIA SUPONITSKY, research paper.
Link for paper: [Linear and nonlinear mechanisms of sound radiation by instability waves in subsonic jets | Journal of Fluid Mechanics | Cambridge Core](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/linear-and-nonlinear-mechanisms-of-sound-radiation-by-instability-waves-in-subsonic-jets/408DB2A13CDB2C74F58C6AB6034BFCCC)
I wasn’t able to match the result.
I have attached the code and the current result plot for reference.
Thanks in advance.
```auto
//x=z, r=y, z=theta
// Run with MPI: ff-mpirun -np 4 jetstab.edp
load "PETSc-complex"
load "msh3"
macro dimension()3L//
include "macro_ddm.idp"
func Pb2 = [P2, P2];
func Pc = [P2, P2, P2, P2, P2]; // [rho, Vr, Vtheta, Vz, Tt]
real rmin = 1e-6; // avoid singularity at r=0
//real rmin = 0;
real rmax = 10;
border rline(t=rmin,rmax){x=0; y=t;}
int np = 45000;
meshL thL=buildmeshL(rline(np));
meshL thLzero = thL;
// Define separate complex FE functions
fespace VhC(thL,Pc);
VhC [rho, vr, vtheta, vz, Tt], [X1, X2, X3, X4, X5];
fespace VhBzero(thLzero,Pb2);
fespace Vh(thL,P1);
Vh U,Ur,Urr;
// Define r as the coordinate variable
Vh r = y;
// *********for the profile********
/*
U=1/(1+(r*r))^2;
ofstream fout("U_output.txt"); // open file for writing
fout << U[] << endl;
*/
real zpos = 0; // renamed from z
real vcoflow = 0.01; // fixed incomplete line
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 r) {
return 0.5 * (tanh((r + a) / delta) - tanh((r - a) / delta));
}
real a0 = 0.59 + 0.09 * tanh(sqrt(0) - 2.9);
real delta0 = (39 + 24*0 + 0.11*0^4) / (1000 + 0^3);
func real Vzref(real r) {
return 0.5 * (tanh((r + a0) / delta0) - tanh((r - a0) / delta0));
}
func real Vzsim(real r) {
return (Vz(r) + vcoflow) / (Vzref(0) + vcoflow);
}
U = Vzsim(r);
ofstream fout("U_output.txt"); // open file for writing
fout << U[] << endl;
// Define dr as a macro for dx
macro dr(u) dy(u) // EOM
Ur=dr(U);
Urr=dr(Ur);
int[int] n2othL2;
macro thL2N2O()n2othL2//
fespace Bh(thL,P2);
real Pr=0.72;
real ga=1.4;
real M=0.9;
real M2=M*M;
Vh T,Tr,Trr;
T=1+M2*(ga-1)*(U-U*U);
Tr=M2*(ga-1)*(dr(U)-2*U*dr(U));
Trr=M2*(ga-1)*(dr(Ur)-2*dr(U)*dr(U)-2*U*dr(Ur));
real C=0.5;
Vh Mu,Mut,K,Kt,Mutt,Ktt;
Mu=(T^1.5)*((1+C)/(T+C));
K=(T^1.5)*((1+C)/(T+C));
Mut= 1.5*sqrt(T)*(1+C)*(1/(T+C)) - (T^(1.5))*((1+C)/((T+C)^(2)));
Mutt=(3/4)*(1/sqrt(T))*((1+C)/(T+C)) -3*sqrt(T)*( (1+C)/((T+C)*(T+C)) ) +2*(T^1.5)*((1+C)/(T+C)^3);
Kt=1.5*sqrt(T)*(1+C)*(1/(T+C)) - (T^(1.5))*((1+C)/((T+C)^(2)));
Ktt=(3/4)*(1/sqrt(T))*((1+C)/(T+C)) -3*sqrt(T)*( (1+C)/((T+C)*(T+C)) ) +2*(T^1.5)*((1+C)/((T+C)^3));
// Viscosity derivatives with respect to r
Vh Mur = Mut * Tr; // ∂μ/∂r = (∂μ/∂T)(∂T/∂r)
Vh Murr = (Mutt * (Tr * Tr)) + (Mut * Trr); // ∂²μ/∂r² = (∂²μ/∂T²)(∂T/∂r)² + (∂μ/∂T)(∂²T/∂r²)
Vh Kr = Kt * Tr; // First derivative ∂κ/∂r
Vh Krr =( Ktt * (Tr * Tr)) + (Kt * Trr); // Second derivative ∂²κ/∂r²
// Temperature and density profiles
Vh Rho,Rhor,Rhorr;
Rho=1/T;
Rhor=dr(Rho);
Rhorr=dr(Rhor);
/* ################################################################################## */
/* ########################## Solution of the Eigenproblem ########################## */
/* ################################################################################## */
real Re = getARGV("-Re", 3600);
real om= getARGV("-om", 9);
real shiftr = getARGV("-shiftr", 15);
real shifti = getARGV("-shifti", -1.4);
complex shift = shiftr + 1i*shifti;
int nEV = getARGV("-nEV", 15);
int nKryl = getARGV("-nKryl", 200);
real m = getARGV("-m", 0); // Azimuthal wavenumber // m = 0 - axisymmetric;m = 1-first helical mode; m = 2 - higher modes
// Define test functions
VhC [a1, a2, a3, a4, a5];
// Weak form definitions based on the PDF formulation
varf A0mat([rho, vr, vtheta, vz, Tt],[a1,a2,a3,a4,a5]) =
int1d(thL) (
// Row 1: Continuity Equation (ν1)
( a1 * (( vr *r*r* Rhor) + (Rho *r*r* dr(vr)) + (Rho * r* vr) + (1i * m *r * Rho * vtheta )) ) - (a1* 1i * om*r*r* rho) // contribution from -i*omega*K
// Row 2: r-Momentum Equation (ν2)
+(a2 * ( ((1.0 / (ga * M * M)) *r*r * Tr * rho)
+ ((1.0 / (ga * M * M)) *r*r *T * dr(rho))
- (1.0 / Re) * ( (-(7.0/3.0) * (1i * m * Mu * vtheta))
- ((2.0/3.0) * (1i * m* r* Mur * vtheta))
+ ((1.0/3.0) * (1i * m * r* Mu * dr(vtheta))) )
+ ((1.0 / (ga * M * M)) *r*r * Rhor * Tt )
+ ((1.0 / (ga * M * M)) *r*r* Rho * dr(Tt))
- (1.0 / Re) *( ((4.0/3.0) * (Mu* r * dr(vr)) )
+ ((4.0/3.0) *r*r* Mur * dr(vr))
- (m*m * Mu * vr)
- ((4.0/3.0) * (Mu * vr))
- ((2.0/3.0) * Mur* r* vr ) ) ) )
+ (1.0 / Re) * (4.0/3.0) * ( (r*r* dr(a2) * Mu * dr(vr))+ ( r*r* a2 * Mur * dr(vr)) + (2* a2*r* Mu* dr(vr)) )
- (a2* 1i * r* r*om * Rho * vr)
// Row 3: θ-Momentum Equation (ν3)
+ (a3 * ( ((1i * m / (ga * M * M))*r * T * rho)
+ ((1i * m / (ga * M * M ))*r *Rho * Tt)
- (1.0 / Re) * ( ((7.0/3.0) * (1i * m * Mu * vr))
+ (1i * m *r * Mur * vr)
+ (1i * m* r*(1.0/3.0) * Mu * dr(vr)) )
- (1.0 / Re) *( (Mu *r * dr(vtheta))
+ (r*r* Mur * dr(vtheta))
- ((4.0/3.0) * (m*m * Mu * vtheta))
- (Mu * vtheta)
- (r* Mur * vtheta) ) ) )
+ (1.0 / Re) * ( dr(a3) *r*r* Mu * dr(vtheta) + a3 * r*r*Mur * dr(vtheta) + 2*a3*Mu*r* dr(vtheta) )
-( a3* 1i *r* r* om * Rho * vtheta)
// Row 4: z-Momentum Equation (ν4)
+ (a4* ( ( Rho * Ur *r*r* vr)
- (1.0 / Re) * ( (Mu *r * dr(vz))
+ (r*r * Mur * dr(vz))
- (m*m * Mu * vz))
- (1.0 / Re) * ((r*r* Urr * Mut * Tt)
+ (r* Mut * Ur*Tt)
+ (r*r* Mut*Ur * dr(Tt)) ) ) )
+ (1.0 / Re) * ( dr(a4) *r*r* Mu * dr(vz) + a4 *r*r* Mur * dr(vz) + 2*a4*r* Mu* dr(vz))
- (a4* 1i *r* r*om * Rho * vz)
// Row 5: Energy Equation (ν5)
+ (a5* ( (Rho * Tr *r*r* vr)
+ ((ga - 1.0) * r * Rho * T* vr)
+ ((ga - 1.0) * Rho * T *r*r* dr(vr))
+ ((ga - 1.0) * (1i * m *r * Rho * T * vtheta))
- ( ((M*M * ga * (ga - 1.0)) / Re ) *r*r* ( 2.0 * Mu * Ur * dr(vz) ))
- (((M*M * ga * (ga - 1.0)) / Re )*r*r* ((Ur * Ur) * Mut * Tt ))
- ( ga / (Pr * Re) ) * (( r*r* Trr * Kt * Tt)
+ (r * Kr * Tt)
- (m*m* K * Tt)
+ (r*r*Kr * dr(Tt))
+ (r*r*Kt*Tr * dr(Tt))
+ (K*r* dr(Tt)) ) ) )
+ ( ga / (Pr * Re) ) * ( dr(a5) *r*r* K * dr(Tt) + a5 *r*r* Kr * dr(Tt) + 2* a5*r*K* dr(Tt) )
- (a5* 1i * r* r*om * Rho * Tt)
)+ on((abs(m) == 0)*1, vr=0 , vtheta=0) // |m| = 0 case
+ on((abs(m) == 1)*1, rho=0, vz=0, Tt=0) // |m| = 1 case
+ on((abs(m) >= 2)*1, rho=0, vr=0, vtheta=0, vz=0, Tt=0 ) // |m| ≥ 2
+ on(2, rho=0, vr=0, vtheta=0, vz=0, Tt=0); // r = Rmax
varf A1mat([rho, vr, vtheta, vz, Tt], [a1, a2, a3, a4, a5]) =
int1d(thL)(
(
a1*(1i*U*r*r*rho + 1i*Rho*r*r*vz)
+ (a2*1i*Rho*U*r*r*vr)
+ (a3* 1i*Rho*U*r*r*vtheta)
+ a4*(((1i* T)/(ga*M*M))*r*r*rho + 1i*Rho*U*r*r*vz + ((1i* Rho)/(ga*M*M)*r*r*Tt ))
+a5* (((ga-1)*1i*r*r*Rho*T*vz)+ 1i*Rho*r*r*U*Tt)
- (1/Re)*(
a2 * (((1.0 / 3.0)* 1i *Mu *r*r* dr(vz))
- ((2.0 / 3.0) * 1i *r*r *Mur * vz)
+ (1i *r*r* Ur * Mut * Tt) )
- a3 * ((1.0 / 3.0) *m *r*r * Mur * vz)
+ a4 * ( (1i * Mu*(1.0/3.0) *r*r * dr(vr))
+ ((1.0 / 3.0)*1i* r* Mu*vr)
+ (1i * Mur*r*r* vr)
- ((m* Mu*r* (1.0/3.0) * dr(vtheta))) )
+ a5*(2.0*M*M*ga*(ga-1.0)*1i*Mu*r*r*Ur*vr)
)
)
)+ on((abs(m) == 0)*1, vr=0 , vtheta=0) // |m| = 0 case
+ on((abs(m) == 1)*1, rho=0, vz=0, Tt=0) // |m| = 1 case
+ on((abs(m) >= 2)*1, rho=0, vr=0, vtheta=0, vz=0, Tt=0 ) // |m| ≥ 2
+ on(2, rho=0, vr=0, vtheta=0, vz=0, Tt=0); // r = Rmax
varf A2mat([rho, vr, vtheta, vz, Tt],[a1,a2,a3,a4,a5]) =
int1d(thL) (
(
-(1.0/Re)*( ( Mu*a2*vr*r*r)
+ (Mu*a3*vtheta*r*r)
+ ((4.0/3.0)*Mu*a4*vz*r*r)
+ ((ga/(Re*Pr))*a5*r*r*K*Tt) )
)
) + on((abs(m) == 0)*1, vr=0 , vtheta=0) // |m| = 0 case
+ on((abs(m) == 1)*1, rho=0, vz=0, Tt=0) // |m| = 1 case
+ on((abs(m) >= 2)*1, rho=0, vr=0, vtheta=0, vz=0, Tt=0) // |m| ≥ 2
+ on(2, rho=0, vr=0, vtheta=0, vz=0, Tt=0); // r = Rmax
Mat[int] Apep(3);
macro def(u)[u, u#B, u#C,u#D,u#E]//EOM
macro init(i)[i, i, i,i,i]//EOM
{
MatCreate(thL, Apep[0], Pc); // Create the matrices array first
MatCreate(thL, Apep[1], Pc);
MatCreate(thL, Apep[2], Pc);
}
Apep[2] = A2mat(VhC,VhC,tgv=-1);
Apep[1] = A1mat(VhC,VhC,tgv=-1);
Apep[0] = A0mat(VhC,VhC,tgv=-1);
VhC[int] def(eigenvecVec)(nEV);
complex[int] eigvalVec(nEV);
real[int] errestVec(nEV);
string PEPParamsbase =
" -pep_basis monomial " +
" -pep_general " +
" -st_type sinvert " +
" -st_pc_type lu " +
" -pep_tol 1e-8 ";
string PEPspar = PEPParamsbase + " -pep_target " + shift +
" -pep_nev " + nEV +
" -pep_ncv " + nKryl + " ";
int nEvalConv = PEPSolve(Apep, vectors = eigenvecVec, values = eigvalVec, sparams = PEPspar, errorestimate = errestVec);
if(mpirank == 0) {
cout << "Same output manually" << endl;
cout << " PEP nconv=" << nEvalConv << " Values (Errors)";
for(int i = 0; i < nEvalConv; i++) {
cout << " " << real(eigvalVec[i]) << "+" << imag(eigvalVec[i]) << "i (" << errestVec[i] << ")";
}
cout << endl;
}
if(mpirank == 0) {
ofstream ofile("Blasius_EV.dat");
ofile.precision(16);
ofile << eigvalVec << endl;
}
/*---------------------------------------------------------------------
Extract the components of the first eigenvector
from the 5-component FE space: [ρ, v_r, v_θ, v_z, T_t]
---------------------------------------------------------------------*/
// Define complex finite element functions for each component
Bh Rh, Vr, Vtheta,Vzz ,Tth;
Rh=eigenvecVec[0]; // density perturbation (rho) component of the first eigenmode
Vr=eigenvecVecB[0]; //x-velocity perturbation component
Vtheta=eigenvecVecC[0]; // y-velocity perturbation component
Vzz=eigenvecVecD[0]; // z-velocity perturbation component
Tth=eigenvecVecE[0]; //temperature/energy perturbation component
//---------------------------------------------
// Write results to files for visualization
//---------------------------------------------
ofstream fileRho("rho_eig.txt");
fileRho.precision(12);
fileRho << Rh[] << endl;
ofstream fileVr("vr_eig.txt");
fileVr.precision(12);
fileVr << Vr[] << endl;
ofstream fileVtheta("vtheta_eig.txt");
fileVtheta.precision(12);
fileVtheta << Vtheta[] << endl;
ofstream fileVzz("vzz_eig.txt");
fileVzz.precision(12);
fileVzz << Vzz[] << endl;
ofstream fileTt("Tt_eig.txt");
fileTt.precision(12);
fileTt << Tth[] << endl;
```