1/****************************************************************************/
2/* This file is part of FreeFEM.                                            */
3/*                                                                          */
4/* FreeFEM is free software: you can redistribute it and/or modify          */
5/* it under the terms of the GNU Lesser General Public License as           */
6/* published by the Free Software Foundation, either version 3 of           */
7/* the License, or (at your option) any later version.                      */
8/*                                                                          */
9/* FreeFEM is distributed in the hope that it will be useful,               */
10/* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
11/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
12/* GNU Lesser General Public License for more details.                      */
13/*                                                                          */
14/* You should have received a copy of the GNU Lesser General Public License */
15/* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.          */
16/****************************************************************************/
17// Author: F. Hecht
18// Stationnary imcompressible Navier Stokes Equation with Newton method (jan. 2012)
19
20verbosity=0;
21
22// Parameters
23func bb = [[-1, -2], [4, 2]]; // bounding box for the plot
24real R = 5, L = 15;
25// physical parameters
26real nu = 1./50, nufinal = 1/200., cnu = 0.5;
27// stop test for Newton
28real eps = 1e-4;
29
30// Mesh
31border cc(t=0, 2*pi) {x=cos(t)/2; y=sin(t)/2; label=1;}
32border ce(t=pi/2, 3*pi/2) {x=cos(t)*R; y=sin(t)*R; label=1;}
33border beb(tt=0, 1) {real t=tt^1.2; x=t*L; y=-R; label=1;}
34border beu(tt=1, 0) {real t=tt^1.2; x=t*L; y=R; label=1;}
35border beo(t=-R, R) {x=L; y=t; label=0;}
36border bei(t=-R/4, R/4) {x=L/2; y=t; label=0;}
37mesh Th = buildmesh(cc(-50) + ce(30) + beb(20) + beu(20) + beo(10) + bei(10));
38plot(Th);
39
40// Macro
41macro Grad(u1, u2) [dx(u1), dy(u1), dx(u2), dy(u2)]//
42macro UgradV(u1, u2, v1, v2) [[u1, u2]'*[dx(v1), dy(v1)], [u1, u2]'*[dx(v2), dy(v2)] ]//
43macro div(u1, u2)  (dx(u1) + dy(u2))//
44
45//  FE Space
46fespace Xh(Th, P2);
47Xh u1, u2;
48Xh v1, v2;
49Xh du1, du2;
50Xh u1p, u2p;
51
52fespace Mh(Th, P1);
53Mh p, q;
54Mh dp, pp;
55
56// Intial guess with B.C.
57u1 = (x^2+y^2) > 2;
58u2 = 0;
59
60// Loop on vicosity
61while(1) {
62	int n;
63	real err = 0;
64	// Newton Loop
65	for (n = 0; n < 15; n++) {
66		solve Oseen([du1, du2, dp], [v1, v2, q])
67			= int2d(Th)(
68					nu*(Grad(du1,du2)'*Grad(v1,v2))
69				+ UgradV(du1,du2, u1, u2)'*[v1,v2]
70				+ UgradV( u1, u2,du1,du2)'*[v1,v2]
71				- div(du1,du2)*q - div(v1,v2)*dp
72				- 1e-8*dp*q // stabilization term
73			)
74			- int2d(Th)(
75				  nu*(Grad(u1,u2)'*Grad(v1,v2))
76				+ UgradV(u1,u2, u1, u2)'*[v1,v2]
77				- div(u1,u2)*q - div(v1,v2)*p
78				- 1e-8*p*q
79			)
80			+ on(1,du1=0,du2=0)
81			;
82
83		u1[] -= du1[];
84		u2[] -= du2[];
85		p[] -= dp[];
86
87		real Lu1 = u1[].linfty, Lu2 = u2[].linfty, Lp = p[].linfty;
88		err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
89
90		cout << n << " err = " << err << " " << eps << " rey  = " << 1./nu << endl;
91		if(err < eps) break; // converge
92		if( n > 3 && err > 10.) break; // Blowup ?
93	}
94	if(err < eps) {// if converge  decrease nu (more difficult)
95		plot([u1, u2], p, wait=1, cmm=" rey = " + 1./nu, coef=0.3, bb=bb);
96		if(nu == nufinal) break;
97		if(n < 4) cnu = cnu^1.5; // fast converge => change faster
98		nu = max(nufinal, nu*cnu); // new vicosity
99		u1p = u1;
100		u2p = u2;
101		pp = p; // save correct solution
102	}
103	else {  // if blowup, increase nu (more simple)
104		assert(cnu < 0.95); // final blowup
105		nu = nu/cnu; //  get previous value of viscosity
106		cnu = cnu^(1./1.5); // no conv. => change lower
107		nu = nu* cnu;  // new vicosity
108		cout << " restart nu = " << nu << " Rey = " << 1./nu << "  (cnu = " << cnu << " ) \n";
109		// restore correct solution
110		u1 = u1p;
111		u2 = u2p;
112		p = pp;
113	}
114}
115