Newton 2

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\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} & =& \vec{f},\quad \Omega\\ \vec{u} & =& 0,\quad \Gamma_D\\ \frac{\partial{\vec{u}}}{\partial\vec{n}} & =& 0,\quad \Gamma_N \end{array}

include Dirichlet boundary conditions on $U$ and take

$$F(\vec{u}):=\alpha\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} -\vec{f}$$

a simple computation gives

$$<dF(\vec{u}),\vec{h}>=\alpha\vec{h}-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}$$

and Newton method is

Want $\vec{u}_\star\in U$ such that $F(\vec{u}_\star)=0$. Given $\vec{u} \approx \vec{u}_\star$ improve $\vec{u}$ by taking $\vec{u}:=\vec{u}+\vec{h}$ where $\vec{h}$ satisfies the linear problem \begin{array}{rcl} \alpha\vec{h}-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u} & =& -(\alpha\vec{u}-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} -\vec{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 \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u} +\vec{\nabla p}& =& \vec{f},\quad \Omega\\ div(\vec{u})& =& 0,\quad \Omega\\ \vec{u} & =& 0,\quad \Gamma_D\\ \frac{\partial{\vec{u}}}{\partial\vec{n}} & =& 0,\quad \Gamma_N \end{array}

include Dirichlet boundary conditions on $U$ space for the pair $(\vec{u},p)$ and maintain the restriction $div( \vec{u})=0$. Then take

$$F(\vec{u},p):=(-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} -\vec{f},\; div(\vec{u}))$$

again, a simple computation gives

$$<dF((\vec{u},p)),(\vec{h},q)>=(-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q},\; div(\vec{h}))$$

and Newton method is

Want $(\vec{u}_\star,p_\star)\in U$ such that $F(\vec{u}_\star,p_\star)=0$. Given $(\vec{u},p) \approx (\vec{u}_\star,p_\star)$ improve $(\vec{u},p)$ by taking $\vec{u}:=\vec{u}+\vec{h}$ and $p:=p+q$ where $(\vec{h},q)$ satisfies the linear problem \begin{array}{rcl} (-\mu \Delta \vec{h}+ (\vec{u}\cdot \vec{\nabla})\vec{h} +(\vec{h}\cdot \vec{\nabla})\vec{u}+\vec{\nabla q}, \; div(\vec{h}) ) & =& (-\mu \Delta \vec{u}+ (\vec{u}\cdot \vec{\nabla})\vec{u}+\vec{\nabla p} -\vec{f},\; div(\vec{u})). \end{array} Solved by a variational formulation using a pair of test function $(\vec{v} ,s)$.

Incompressible Navier Stokes around a 3d Cylinder. Continuation on viscosity

Adapted on April 2014 from code by Author: F. Hecht, jan 2012.

 Triangulation Viscosity 0.02 Viscosity 0.002
 Close up viscosity 0.02 Close up viscosity 0.002

//
// adapted on April 2014 from code by Author: F. Hecht
// jan 2012
//
//Stationnary imcompressible Navier Stokes Equation with Newton method.
//  a round a 3d Cylinder
// build the  Mesh
//

real R = 5,L=15;
border cc(t=0,2*pi){ x=cos(t)/2;y=sin(t)/2;label=1;}
border ce(t=pi/2,3*pi/2) { x=cos(t)*R;y=sin(t)*R;label=1;}
//cluster triangulation points on boundary y=-R near x=0
border beb(tt=0,1) { real t=tt^1.2; x= t*L; y= -R; label = 1;}
//cluster triangulation points on boundary y=R near x=0
border beu(tt=1,0) { real t=tt^1.2; x= t*L; y= R; label = 1;}
border beo(t=-R,R) {  x= L; y= t; label = 0;}
/*
//
// to see the boundary curves one by one
//
plot(cc(-50),wait=1);
plot(cc(-50)+ce(30),wait=1);
plot(cc(-50),ce(30),beb(20),wait=1);
plot(cc(-50),ce(30),beb(20),beu(20),wait=1);
plot(cc(-50),ce(30),beb(20),beu(20),beo(10),wait=1);
*/
int mm=30;
mesh Th=buildmesh(cc(-2*mm)+ce(3*mm)+beb(3*mm)+beu(3*mm)+beo(2*mm));
plot(cmm=" This is the mesh...Wait for computations...",Th);
//
// for a close view
//
func bb=[[-1,-2],[4,2]];
//
// Macro operators
//
macro Grad(u1,u2) [ dx(u1),dy(u1), dx(u2),dy(u2)]//
macro UgradV(u1,u2,v1,v2) [ [u1,u2]'*[dx(v1),dy(v1)] , [u1,u2]'*[dx(v2),dy(v2)] ]//
macro div(u1,u2)  (dx(u1)+dy(u2))//

//  FE Space
fespace Xh(Th,P2);// Use P1b for a less expensive computation
fespace Mh(Th,P1);

Xh u1,u2;// solutions for velocity variables
Xh du1,du2;// increments for velocity variables
Xh uc1,uc2;// computed solutions for velocity variables at previous viscosity
Mh p;// solution for the pressure variable
Mh dp;// increment for the pressure variable
Xh v1,v2; // test functions for the velocity variables
Mh q;// test function for pressure variable
Mh pc;//  computed solutions for pressure variable at previous viscosity

// intial guess with B.C.

u1 = ( x^2+y^2) > 2;
u2=0;
p=0;
//  Physical parameter
real nu= 1./50; // starting viscosity
real nufinal=1/500.; // wanted final viscosity
//
// cnu reducing factor for viscosity
// We reduce viscosity by dividing by cnu
//
real cnu=2;
//
// Number of reductions on viscosity
//
int countermax=10;
int counter=0;
// stop test for Newton
real eps=1e-6;

//
verbosity=0;
real err=0;
real nuc=nu;// current viscosity
cout <<  "  ----------------------------  " << endl;
cout <<  "  --- START COMPUTATIONS ---  " << endl;
cout <<  "  ----------------------------  " << endl;
//
//  Loop on vicosity until reach nufinal, countermax computations or stop reducing
//
while((nuc>nufinal) && (counter<countermax)&&(cnu>1.05))
{ // start loop on viscosity
counter =counter+1;
int n; // need value of n out of the for loop
cout <<  " Compute with viscosity  " << nuc <<" (counter "<<counter<<")"<< endl;
for( n=0;n< 15;n++) // Newton Loop for current value of cnu
{ // --->star for loop
//   compute increments for the three variables u1,u2,p
solve Oseen([du1,du2,dp],[v1,v2,q],init=n) =
+ UgradV(du1,du2, u1, u2)'*[v1,v2]
+ UgradV( u1, u2,du1,du2)'*[v1,v2]
- div(du1,du2)*q - div(v1,v2)*dp
- 1e-8*dp*q // stabilization term
)
+ UgradV(u1,u2, u1, u2)'*[v1,v2]
- div(u1,u2)*q - div(v1,v2)*p
- 1e-8*p*q
)
+ on(1,du1=0,du2=0)
;

//
// update values u1, u2,p to compute again
//
u1[] += du1[]; //u1[]=u1[]+du1[]
u2[] += du2[]; //u2[]=u2[]+du2[]
p[]  += dp[];  //p[]=p[]+dp[]

//
// computation of norms for the increments
//
real Lu1=u1[].linfty,  Lu2 = u2[].linfty , Lp = p[].linfty;
//
// normalization of the norms
//
err= du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
cout<<" it = "+ n+ ", err = "+err <<endl;
//
// plot iterates
//
plot([u1,u2],p,wait=0,cmm="Newton Method for viscosity = " + nuc +", it = "+ n+ ", err = "+err, coef=0.3);
//
//--> In case of no convergence...
//      (n>3 AND err>2) is blowup
//      (n>12 AND err>eps) is lack of convergence or very slow
//
if( (n>3 && err > 2.) || (n>12 && err > 1e4*eps) )  // (n>3 AND err>2) or (n>12 AND err>eps)
{// start ---> if
cout<<" NO convergence for nu = "+nuc+ ", it = "+ n+ ", err = "+err <<endl;
nuc=nuc*cnu; // get back to previous nu
cout << " LAST viscosity TO CONVERGE = " << nuc<< endl;
//
// take computed values for previous nuc to start again
//
u1=uc1; u2=uc2; p=pc;
//
// New reduction for cnu
//
cnu= cnu^(0.6); // no conv. => change lower (cnu^(0.6) makes cnu smaller)
nuc = nuc/cnu;  // new vicosity
cout << " RESTART with NEW viscosity  = " << nuc << endl;
cout << " Using as Inital Data the solution for PREVIOUS viscosity "  << endl;
break; //  Blowup ???? --->out of  for loop  go to while loop
}// end ---> if( n>3 && err > 2.)
//
//-->In case of Convergence then...
//
if(err < eps)
{// start ---> if(err < eps)
cout << " CONVERGENCE  for viscosity = " << nuc << " on iters = "<< n <<", err = "+err << endl;
uc1=u1; uc2=u2; pc=p; // save computed values of this convergent solution
plot([u1,u2],p,wait=1,cmm="Solution for viscosity = " + nuc+". Press Enter for a close up ",coef=0.3);
plot([u1,u2],p,wait=1,coef=0.3,bb=bb);
if( n < 4)  cnu=cnu^1.5; // fast converge => change faster ( cnu^ 1.5 makes bigger cnu)
nuc = max(nufinal, nuc/cnu); // reduce to a new vicosity
cout << " Continuation to new viscosity  = " << nuc <<" (dividing factor "<<cnu<<")"<< endl;
cout << " Initial data IS the FINAL computation  for previous viscosity " << endl;
break; // converge ---> out of for loop. Start again with  u1,u2,p and new nuc
}// end ---> if(err < eps)
}//end ---> for loop
} // end  while loop on viscosity
plot([uc1,uc2],pc,cmm="Last convergence of Newton Method for nu = " + nuc+". Press Enter for a close up ",coef=0.3,wait=1);
plot([uc1,uc2],pc,cmm="Close up nu = " + nuc,coef=0.3,bb=bb);
cout << " LAST viscosity TO CONVERGE nu = " << nuc << endl;

Page last modified on April 15, 2014, at 06:04 PM