Hello all,
I am new to FreeFem++ and I hope anyone here could give some help.
I am trying to export the data from FreeFem++ however I am not able to replicate the results in Octave. Probably due to the way I am trying to export the results.
Taking this example: https://www.um.es/freefem/ff++/pmwiki.php?n=Main.POD
Considering only the High Definition Model (HDM) of the Lid-Driven Cavity:
load “lapack”
int n=20;
int nev=4;
int nsnsh=10;
mesh Th=square(n,n);
//Th=movemesh(Th,[(1-cos(pix))/2.,(1-cos(piy))/2.]); // refine near the boundaries - refinement removed
fespace Vh(Th,[P2,P2,P1]);
Vh [u,v,p],[uu,vv,pp],[up,vp,q];
// snapshot space:
Vh[int] usnsh,vsnsh,psnsh;
// average values:
Vh [uavg,vavg,pavg];
uavg=0.; up=0.;
macro div(u,v)(dx(u)+dy(v))//
macro grad(u)[dx(u),dy(u)]//
int i,j,k;
real Re=1,nu=1/Re,dt=.1;
real[int] viso=[0.33,0.21,0.17,0.12,0.07,0.03,-0.017,-0.0065,-0.112,-0.16,-0.21,-0.25,-0.3,-0.35,-0.4,-0.44,-0.49,-0.54,-0.585,-0.66];
problem NSunst([u,v,p],[uu,vv,pp],solver=sparsesolver)=
int2d(Th)(nu*(grad(u)‘*grad(uu)+grad(v)’*grad(vv)))
+int2d(Th)([up,vp]'grad(u)uu+[up,vp]'grad(v)vv)
-int2d(Th)(pdiv(uu,vv)+ppdiv(u,v))
+int2d(Th)(1.e-10ppp)
+on(1,2,4,u=0,v=0)
+on(3,u=1,v=0)
;
for (j=0;j<nsnsh;j++){
for (i=0;i<10;i++){ // here starts the fixed-point iteration for the nonlinearity
NSunst;
up=u;
}
usnsh[j]=u;
uavg+=usnsh[j];
plot([usnsh[j],vsnsh[j]],cmm=“Reynolds=”+Re,ps=“Re”+j+Re+“.eps”,wait=1);
Re += 50;nu=1/Re;
}
uavg/=nsnsh;
plot([uavg,vavg],value=1,cmm=“uavg”,ps=“avg.eps”,wait=1);
for (int i=0;i<nsnsh;i++)
usnsh[i]-=uavg;
real[int,int] C(nsnsh,nsnsh);
for (i=0;i<nsnsh;i++)
for (j=i;j<nsnsh;j++)
C(j,i) = int2d(Th)(usnsh[i]*usnsh[j]+vsnsh[i]*vsnsh[j]);
for (j=0;j<nsnsh;j++)
for (i=j;i<nsnsh;i++)
C(j,i)=C(i,j);
cout << "correlation matrix: "<< endl;
cout << C << endl;
// identity matrx
real[int,int] II(nsnsh,nsnsh);
II=0;
for (i=0;i<nsnsh;i++) II(i,i)=1.;
real[int] ev(nsnsh);
real[int,int] eVec(nsnsh,nsnsh);
int k1=dsygvd(C,II,ev,eVec);
cout << "eigenvalues of the correlation matrix: " << ev << endl;
real[int,int] matVec(nsnsh,Vh.ndof);
for (i=0;i<nsnsh;i++)
matVec(i,:)=usnsh[i];
real normsn=1;
// recover the modes and store them in the usnsh matrix
for (i=0;i<nev;i++){
usnsh[i]=0.;
for (j=0;j<nsnsh;j++){
usnsh[i]+=eVec(j,nsnsh-1-i)*matVec(j,:); // reordering: most energetic modes first
}
normsn = sqrt(int2d(Th)(square(usnsh[i])+square(vsnsh[i])));
usnsh[i]/=normsn;
}
for(i=0;i<nev;i++) {
plot([usnsh[i],vsnsh[i]],cmm=“mode”+i,ps=“129PODmode”+i+“.eps”,wait=1);
}
NSunst;
I am trying to export by ofstream file:
ofstream file(“datafile.dat”);
for (int i = 0; i < Th.nv; i++) {
real xCoord = Th(i).x;
real yCoord = Th(i).y;
real val01 = usnsh[0](xCoord, yCoord); //[#] to choose the # of snsh so it retrieves the wanted Re response
real val02 = vsnsh[0](xCoord, yCoord); //[#] to choose the # of snsh so it retrieves the wanted Re response
file << xCoord << " " << yCoord << " " << val01 << " " << val02 << endl;
However, the exported file doesn’t match the result given by FreeFem++ and the papers. For instance, at the lid, it should be 1 as u=1 v=0 was imposed at that boundary.
Many thanks in advance for your help.