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