/usr/share/doc/freefem++/examples/examples++-tutorial/NSUzawaCahouetChabart.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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 | // correct bug in CahuetChabard routine
// sign of pressure correct
// and change Bx, By in Bx' By' (back and forth) in version>1.18
// correction in CahouetChabard routine in version 1.26
assert(version>1.18); // for big bug is non symetric matrix see HISTORY
real s0=clock();
mesh Th=square(10,10);
fespace Xh(Th,P2),Mh(Th,P1);
Xh u1,u2,v1,v2;
Mh p,q,ppp;
real[int] pwork(p.n);
int i=0;
real one =1;
varf bx(u1,q) = int2d(Th,qforder=4)( dx(u1)*q );
varf by(u1,q) = int2d(Th,qforder=4)( dy(u1)*q );
varf a(u1,u2)= int2d(Th,qforder=4)( dx(u1)*dx(u2) + dy(u1)*dy(u2) )
+ on(3,u1=1) + on(1,2,4,u1=0) ;
varf vfMass(p,q) = int2d(Th)(p*q);
matrix MassMh=vfMass(Mh,Mh,solver=CG);
p[]=1;
ppp[]= MassMh*p[];
real aire = ppp[].sum;
cout << " area of Omega = " << aire << endl;
Xh bc1; bc1[] = a(0,Xh);
Xh b;
matrix A= a(Xh,Xh,solver=CG);
matrix Bx= bx(Xh,Mh); // Mh corresponding to line and Xh to column
matrix By= by(Xh,Mh);
Xh bcx=(1-x)*x*4,bcy=0;
Xh b1=0,b2=0; // right hand side (0) for Stokes
func real[int] divup(real[int] & pp)
{
pwork= MassMh*pp;
cout << " pp moy = " << pwork.sum/aire << " " ;
pp -= pwork.sum/aire;;
b = b1; b[] += Bx'*pp; b[] += bc1[] .*bcx[];
u1[] = A^-1*b[];
b = b2; b[] += By'*pp; b[] += bc1[] .*bcy[];
u2[] = A^-1*b[];
ppp[] = Bx*u1[];
ppp[] += By*u2[];
pwork= MassMh*ppp[];
cout << " moy = " << pwork.sum/aire << endl;
// ppp[] -= pwork.sum/aire;
return ppp[] ;
};
p=0;q=0;u1=0;v1=0;
//cout << " -------- A = " << A << endl;
LinearCG(divup,p[],q[],eps=1.e-6,nbiter=50);
divup(p[]);
plot([u1,u2],p,wait=1,value=true,coef=0.1);
real dt=0.05, alpha=1/dt;
if ( alpha > 1e30) exit(1);
real xnu=1./400.;
cout << " alpha = " << alpha << " nu = " << xnu << endl;
cout << "------------------------------------------ " << endl;
varf at(u1,u2)= int2d(Th)( xnu*dx(u1)*dx(u2) + xnu*dy(u1)*dy(u2) + u1*u2*alpha )
+ on(3,u1=1) + on(1,2,4,u1=0) ;
A = at(Xh,Xh,solver=CG);
//cout << " -------- AA = " << AA << endl;
varf vfconv1(uu,vv) = int2d(Th,qforder=5) (convect([u1,u2],-dt,u1)*vv*alpha);
varf vfconv2(v2,v1) = int2d(Th,qforder=5) (convect([u1,u2],-dt,u2)*v1*alpha);
int idt;
real temps=0;
Mh pprec,prhs;
varf vfLap(p,q) = int2d(Th)(dx(p)*dx(q) + dy(p)*dy(q) + p*q*1e-3);
matrix LapMh= vfLap(Mh,Mh,solver=Cholesky,factorize=true);
real[int] unMh(pprec.n);
pprec[]=1;
unMh = MassMh*pprec[];
func real[int] CahouetChabart(real[int] & xx)
{ // xx = \int (div u) w_i
// $ \alpha \Delta^{-1} + \nu I $
// $ \alpha \LapMh ^{-1} + \nu MassMh^-1 $
real mxx= unMh'*xx;
xx -= mxx*Th.area;
pprec[]= LapMh^-1* xx;
prhs[] = MassMh^-1*xx;
pprec[] = alpha*pprec[]+xnu* prhs[];
// xx = LapMh*pprec[];
// pprec[] -= xx.sum;
return pprec[];
};
Xh up1,up2;
for (idt = 1; idt < 50; idt++)
{
up1=u1;
up2=u2;
temps += dt;
cout << " --------- temps " << temps << " \n ";
b1[] = vfconv1(0,Xh);
b2[] = vfconv2(0,Xh);
cout << " min b1 b2 " << b1[].min << " " << b2[].min
<< " max b1 b2 " << b1[].max << " " << b2[].max << endl;
LinearCG(divup,p[],q[],eps=-1.e-6,nbiter=50,precon=CahouetChabart);
divup(p[]);
real errl2 = sqrt(int2d(Th)( (u1-up1)^2 + (u2-up2)^2 ) );
cout << " errl2 " << errl2 << endl;
plot([u1,u2],p,wait=!(idt%10),value= 1,coef=0.1,cmm="[u1,u2],p || u^n+1 - u^n ]]_L2 ="+errl2);
if (errl2 < 1e-4) break;
}
|