Newton method for \(F:\mathbb{R}\to \mathbb{R}\) is given by:
Want $ x_\star $ such that \(F(x_\star)=0\). Given \(x \approx x_\star\), by Taylor expansion we have \(F(x_\star)=0\approx F(x)+(x_\star-x)F'(x)\) . Then, we improve \(x\) by taking \(x:=x+h\) where \(h=x_\star-x\) satisfies \begin{array}{rcl} F'(x)h & =& -F(x). \end{array}
Take now \(F:U\to V\) where \(U\) and \(V\) are Banach spaces. The Gateaux derivative of \(F\) at a point \(x_\star\in U\) is the linear mapping \(dF(x_\star):U\to V\) that satisfies (denote the value of \(dF(x_\star) \) at \(h\in U\) by \(<dF(x_\star),h>\in V\) )
$$ \begin{array}{rcl} <dF(x_\star),h> & :=&\lim_{s\to 0}\frac{1}{s}(F(x_\star+sh)-F(x_\star)),\quad \forall h\in U. \end{array} $$
Newton method for \(F:U\to V\) is given by:
Want \(x_\star\in U\) such that \(F(x_\star)=0\). Given \(x \approx x_\star\) improve \(x\) by taking \(x:=x+h\) where \(h\) satisfies \begin{array}{rcl} <dF(x),h> & =& -F(x). \end{array}
Example 1. For the nonlinear advection equation on \(\Omega\) with \(\partial\Omega=\Gamma_D\cup\Gamma_N\)
\begin{array}{rcl} \alpha\,u-\mu \Delta u+ (u\cdot \nabla)u & =& f,\quad \Omega\\ u & =& 0,\quad \Gamma_D\\ \frac{\partial{u}}{\partial n} & =& 0,\quad \Gamma_N \end{array}
include Dirichlet boundary conditions on \(U\) and take
$$F(u):=\alpha u-\mu \Delta u+ (u\cdot \nabla)u -f$$
a simple computation gives
$$<dF(u),h>=\alpha h-\mu \Delta h+ (u\cdot \nabla)h +(h\cdot \nabla)u $$
and Newton method is
Want \(u_\star\in U\) such that \(F(u_\star)=0\). Given \(u \approx u_\star\) improve \(u\) by taking \(u:=u+h\) where \(h\) satisfies the linear problem \begin{array}{rcl} \alpha h-\mu \Delta h+ (u\cdot \nabla)h +(h\cdot \nabla)u & =& -(\alpha u-\mu \Delta u+ (u\cdot \nabla) u -f). \end{array}
Example 2. For the nonlinear Navier-Stokes equation on \(\Omega\) with \(\partial\Omega=\Gamma_D\cup\Gamma_N\)
\begin{array}{rcl} -\mu \Delta u+ (u\cdot \nabla)u +\nabla p& =&f,\quad \Omega\\ div(u)& =& 0,\quad \Omega\\ u& =& 0,\quad \Gamma_D\\ \frac{\partial{u}}{\partial n} & =& 0,\quad \Gamma_N \end{array}
include Dirichlet boundary conditions on \(U\) space for the pair \((u,p)\) and maintain the restriction \(div( u)=0\). Then take
$$F(u,p):=(-\mu \Delta u+ (u\cdot \nabla)u+\nabla p -f,\; div(u))$$
again, a simple computation gives
$$<dF((u,p)),(h,q)>=(-\mu \Delta h+ (u\cdot \nabla)h +(h\cdot \nabla)u+\nabla q,\; div(h)) $$
and Newton method is
Want \((u_\star,p_\star)\in U\) such that \(F(u_\star,p_\star)=0\). Given \((u,p) \approx (u_\star,p_\star)\) improve \((u,p)\) by taking \(u:=u+h\) and \(p:=p+q\) where \((h,q)\) satisfies the linear problem \begin{array}{rcl} (-\mu \Delta h+ (u\cdot \nabla)h +(h\cdot \nabla)u+\nabla q, \; div(h) ) & =& (-\mu \Delta u+ (u\cdot \nabla)u+\nabla p -f,\; div(u)). \end{array} Solved by a variational formulation using a pair of test function \((v ,s)\).
Incompressible Navier Stokes with Taylor-Hood Finite element No linearity : Newton method continuation on Reynols Number and Mesh adaptation
Test case Laminar Flow over a Backward Facing Step Gamm Workshop
K.Morgan, J.Périaux and F.Thomasset, Analysis of laminar flow over a backward facing step, Vol9 of Notes on Num. Fluid Mech., Vieweg, 1984.
![]() | ![]() |
Reynolds 300 | Reynolds 500 |
download example: NS-BackwardStep.edp or return to main page
/* Incompressible Navier Stokes with Taylor-Hood Finite element No linearity : Newton method continuation on Reynols Number Mesh adaptation Test case Laminar Flow over a Backward Facing Step Gamm Workshop K.Morgan, J.Périaux and F.Thomasset, Analysis of laminar flow over a backward facing step, Vol9 of Notes on Num. Fluid Mech., Vieweg, 1984. */ real[int] Reynold=[50,150,300,400,500]; real[int] HH=[1.5,1]; real[int,int] reattachP=[ [ 2.8, 2 ], [ 5.16, 3.7 ]] ; // reattachP[irey,cas] int nerr=0; bool adapt=1; // do adap or not bool dplot=0; // debug plot for(int cas=0;cas<2;++cas) { real h=HH[cas]-0.5,H=HH[cas],l=3,L=22; int[int] ll=[3,2,5,1]; // label 1 in // 2 out // 3 down wall // 4 step // 5 up wall func zoom=[[2.5,0],[10,H]]; // // A way to construct the geometry of the step not using // a parametric description of the boundary // mesh Th=square(6,22,[x*22,y*H],label=ll); plot(Th,cmm=" ORIGINAL MESH ",wait=1); Th=trunc(Th, (x>l) | (y >0.5),label=4); plot(Th,cmm=" TRUNCATED MESH ",wait=1); func meshsize= 2*max(0.05,max(max(x-l,0.0)/19./5.,max(l-x,0.0)/3./8. )); func uin=(H-y)*(y-0.5)/square((H-0.5)/2.); Th=adaptmesh(Th,meshsize,IsMetric=1); plot(Th,cmm=" ADAPTED MESH FIRST TIME...",wait=1); Th=adaptmesh(Th,meshsize,IsMetric=1); plot(Th,cmm=" ADAPTED MESH SECOND TIME...",wait=1); // // Finite element spaces // fespace Xh(Th,P2); fespace Mh(Th,P1); fespace XXMh(Th,[P2,P2,P1]); XXMh [u1,u2,p]; XXMh [v1,v2,q]; // // macros for the variational formulation // macro div(u1,u2) (dx(u1)+dy(u2))// macro grad(u1,u2) [dx(u1),dy(u2)]// macro ugrad(u1,u2,v) (u1*dx(v)+u2*dy(v)) // macro Ugrad(u1,u2,v1,v2) [ugrad(u1,u2,v1),ugrad(u1,u2,v2)]// // // Stokes problem // solve Stokes ([u1,u2,p],[v1,v2,q],solver=UMFPACK) = int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) ) + p*q*(0.000001) - p*div(v1,v2)-q*div(u1,u2) ) + on(1,u1=uin,u2=0) + on(3,4,5,u1=0,u2=0); Xh uu1=u1,uu2=u2; plot(coef=0.2,cmm="Stokes [u1,u2] et p ",p,[uu1,uu2],wait=0); plot(coef=0.2,cmm="Stokes p ",p,wait=0); Mh psi,phi; // // Stream function // solve streamlines(psi,phi) = int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi)) + int2d(Th)( -phi*(dy(u1)-dx(u2))) + on(3,4,psi=0)+ on(5,psi=-2./3.*(H-0.5)) ; // // To see predefined contour lines // real[int] psiviso(31); {int k=0; for(int i=-20;i<0;i++) psiviso[k++] = i*2./3.*(H-0.5)/20; for(int i=0;i<=10;i++) psiviso[k++] = i*2./3.*(H-0.5)/100/(H*H*H); } plot(psi,cmm="CONTOUR LINES FOR THE STREAM FUNCTION ",wait=1,viso=psiviso); int i=0; real nu=1./100.; real dt=0.1; real alpha=1/dt; XXMh [up1,up2,pp]; // // To build the matrix for the linear system in Newton method // // this is the matrix that multiply the increment, dF // varf vDNS ([u1,u2,p],[v1,v2,q]) = int2d(Th)( nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) +dx(u2)*dx(v2) + dy(u2)*dy(v2) ) + p*q*1e-8// stabilization term - p*(dx(v1)+dy(v2)) - (dx(u1)+dy(u2))*q + Ugrad(u1,u2,up1,up2)'*[v1,v2] + Ugrad(up1,up2,u1,u2)'*[v1,v2] ) + on(1,3,4,5,u1=0,u2=0) ; // // To build right hand side for the linear system in Newton method // varf vNS ([u1,u2,p],[v1,v2,q]) = int2d(Th)( -nu * ( dx(up1)*dx(v1) + dy(up1)*dy(v1) +dx(up2)*dx(v2) + dy(up2)*dy(v2) ) + pp*q*1e-8// stabilization term + pp*(dx(v1)+ dy(v2)) + (dx(up1)+ dy(up2))*q - Ugrad(up1,up2,up1,up2)'*[v1,v2] ) + on(1,3,4,5,u1=0,u2=0) ; // // Continuation on Reynolds number // // Compute with one Rey. // If convergence, use solution as initial data and increase Rey // for(int krey=0;krey<Reynold.n;++krey) { real re=Reynold[krey]; real lerr=0.01; for(int step=0;step<(adapt?2:1) ;step++) { if(adapt) { Th=adaptmesh(Th,[u1,u2],p,abserror=1,cutoff=1e-5,err=lerr,nbvx=100000,hmin=0.01); if(dplot) plot(Th,wait=0,bb=zoom); } [u1,u2,p]=[u1,u2,p]; [up1,up2,pp]=[up1,up2,pp]; // Newton iterations until convergence for a given Rey // for (i=0;i<=20;i++) { nu = (H-h)/re; up1[]=u1[]; real[int] b = vNS(0,XXMh); // Right hand side for the linear system matrix Ans=vDNS(XXMh,XXMh);// Matrix for the linear system set(Ans,solver=UMFPACK);// Set solver to matrix real[int] w = Ans^-1*b; // Solve the system u1[] += w; // Perform the update of the increment in both variables at the same time cout << " iter = "<< i << " " << w.l2 << " rey = " << re << endl; if(w.l2<1e-6) break; // uu1=u1;uu2=u2; if(dplot) plot(coef=0.2,cmm="H="+H+" re "+re+ " [u1,u2] et p ",p,[uu1,uu2],bb=zoom); } ; } uu1=u1;uu2=u2; streamlines; real rp1=1./(H-h)*int1d(Th,3)( real( (x>=l & x < (l+0.5)) | (x>(l+0.4)) & (x<10)& (dy(psi) >= 1e-5)) ) ; real rp2=1./(H-h)*int1d(Th,3)( real( (x>=l & x < (l+0.5)) | (x>(l+0.4)) & (x<10)& (dy(psi) >= -1e-5)) ) ; real rp3=1./(H-h)*int1d(Th,3)( real( (x>=l & x < (l+0.5)) | (x>(l+0.4)) & (x<10)& (dy(u1)<=0) ) ) ; cout << " Reattach point " << rp2 << " " << rp2 << " " << rp3 << endl; real rp = (rp1+rp2)/2; real rppaper = krey < 2 ? reattachP(krey,cas) : rp; real err= abs(rppaper - rp)/rp; if( err>0.3 ) nerr++;// cout << "\n\n\n"; cout << "H= " << H << " Re " << re << " Reattach point " << rp << " paper=" << rppaper << " err "<< err << " psi max = " << psi[].max <<endl; cout << "\n\n\n"; plot(coef=0.2,cmm="H="+H+", rey="+re+" [u1,u2] et p ",p,[uu1,uu2],wait=0,nbiso=20,bb=zoom);//,ps="Upstep-"+H+"-"+re plot(coef=0.2,cmm="H="+H+", rey="+re+" [u1,u2] et p ",p,[uu1,uu2],wait=0,nbiso=20,bb=zoom);//,ps="Upstep-"+H+"-"+re+".ps"); plot(coef=0.2,cmm="H="+H+", rey="+re+" [u1,u2] et p ",psi,bb=zoom,viso=psiviso);//,ps="psi-step-"+H+"-"+re+".ps"); } } assert(nerr==0);