/usr/share/doc/freefem++/examples/examples++-bug/aaa.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 | // remark: the sign of p is correct
bool classique=0;
real s0=clock();
mesh Th=square(10,10);
fespace Vh2(Th,P2);
fespace Vh(Th,P1);
fespace Wh(Th,[P2,P2,P1]);
Vh2 u2,v2,up1,up2;
Vh2 u1,v1;
Vh u1x=0,u1y,u2x,u2y, vv;
real reylnods=1000;
//cout << " Enter the reynolds number :"; cin >> reylnods;
assert(reylnods>1 && reylnods < 100000);
up1=0;
up2=0;
func g=(x)*(1-x)*4;
Vh p=0,q;
real alpha=0;
real nu=1;
int i=0,iter=0;
real dt=0;
real sig = 2*classique-1;
varf vNS([u1,u2,p],[v1,v2,q]) =
int2d(Th)(
alpha*( u1*v1 + u2*v2)
+ nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
+ p*q*(0.000001)
- p*dx(v1) - p*dy(v2)
- dx(u1)*q - dy(u2)*q
)
+ int2d(Th) ( sig*(-alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 ) )
+ on(3,u1=g,u2=0)
+ on(1,2,4,u1=0,u2=0) ;
solve NS ([u1,u2,p],[v1,v2,q],solver=UMFPACK,init=i,save="toto") = vNS;
cout << "-- n " << p[].n << " stokes " << endl;
cout << "-- u1 : " << u1[].min << " " << u1[].max << endl;
cout << "-- u2 : " << u2[].min << " " << u2[].max << endl;
cout << "-- p : " << p[].min << " " << p[].max << endl;
//plot(coef=0.2,cmm=" [u1,u2] et p ",p,[u1,u2],ps="StokesP2P1.eps",value=1,wait=1);
dt = 0.1;
int nbiter = 2;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter) , cut=0.01;
real tol=0.5,coeftol = 0.25^(1./nbiter);
nu=1./reylnods;
Wh [uu1,uu2,pp];
Wh [vv1,vv2,qq];
Wh [f1,f2,fp];
for (iter=1;iter<=nbiter;iter++)
{
cout << " dt = " << dt << " ------------------------ sig ="<< sig << endl;
alpha=1/dt;
for (i=0;i<=1;i++)
{
up1=u1;
up2=u2;
matrix A=vNS(Wh,Wh,solver=UMFPACK);
// set(A,solver=UMFPACK); // set a solver
verbosity=3;
if(classique) {
//NS;
solve NS1 ([uu1,uu2,pp],[vv1,vv2,qq],solver=UMFPACK,init=i,save="toto") = vNS;
}
else {
f1[] = vNS(0,Wh);
{
ofstream tt("tt.matrix");
tt << A << endl;
}
{
ofstream tt("tt.b");
tt << f1[] << endl;
}
uu1[] = A^-1*f1[];
}
verbosity=1;
u1=uu1;
u2=uu2;
p = pp;
cout << "-- n " << p[].n << endl;
cout << "-- u1 : " << u1[].min << " " << u1[].max << endl;
cout << "-- u2 : " << u2[].min << " " << u2[].max << endl;
cout << "-- p : " << p[].min << " " << p[].max << endl;
if ( !(i % 10))
plot(coef=0.2,cmm=" [u1,u2] et p ",p,[u1,u2],ps="plotNS_"+iter+"_"+i+".eps");
cout << "CPU " << clock()-s0 << "s " << endl;
}
if (iter>= nbiter) break;
Th=adaptmesh(Th,[dx(u1),dy(u1),dx(u1),dy(u2)],
abserror=0,cutoff=cut,err=tol, inquire=0,ratio=1.5,hmin=1./1000);
plot(Th,ps="ThNS.eps");
dt = dt*coefdt;
tol = tol *coeftol;
cut = cut *coefcut;
}
cout << "CPU " << clock()-s0 << "s " << endl;
|