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