Newton 1

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

/*
Incompressible Navier Stokes
with Taylor-Hood Finite element
No linearity : Newton method
continuation on Reynols Number

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

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

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