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);