algo.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

//  cleanning version 07/2008  FH in Sevilla.
real umax=0;
{
  // minimisation of $J(u) = \frac12\sum (i+1) u_i^2 - b_i $	
  // work array 
  real[int] b(10),u(10); 
  
  func real J(real[int] & u)
    {
      real s=0;
      for (int i=0;i<u.n;i++)
	s +=(i+1)*u[i]*u[i]*0.5 - b[i]*u[i];
      cout << "J ="<< s << " u =" <<  u[0] << " " << u[1] << "...\n" ;
      return s;
    }

//  the grad of J (this is a affine version (the RHS is in  )
  func real[int] DJ(real[int] &u)
    { 
      for (int i=0;i<u.n;i++)
	u[i]=(i+1)*u[i]-b[i];
      return u;  // return of global variable ok 
    };

// the grad of the bilinear part of J (the RHS in remove)
  func real[int] DJ0(real[int] &u)
    { 
      for (int i=0;i<u.n;i++)
	u[i]=(i+1)*u[i];
      return u;  // return of global variable ok 
    };


  func real error(real[int] & u,real[int] & b)
   {
   real s=0;
     for (int i=0;i<u.n;i++)
	s += abs((i+1)*u[i] - b[i]);
   return s;    
   }

  func real[int] matId(real[int] &u) { return u;};
  
  b=0.; u=0.; // set  right hand side and initial gest
  LinearCG(DJ,u,eps=1.e-6,nbiter=20,precon=matId);
  cout << "LinearCG (Affin)  : J(u) = " << J(u) << " err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);

  b=1; u=0; // set  right hand side and initial gest
  LinearCG(DJ0,u,b,eps=1.e-6,nbiter=20,precon=matId);
  cout << "LinearCG (Linear) : J(u) = " << J(u) << " err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);

  b=1; u=0; // set  right hand side and initial gest
 //   LinearGMRES(DJ,u,eps=1.e-6,nbiter=20,precon=matId); // bug don't work compile error v1.44
  LinearGMRES(DJ0,u,b,eps=1.e-6,nbiter=20,precon=matId);
  cout << "LinearGMRES: J(u) = " << J(u) << " err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);

  b=1; u=0; // set  right hand side and initial gest
  NLCG(DJ,u,eps=1.e-6,nbiter=20,precon=matId);
  cout << "NLCG: J(u) = " << J(u) << " err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);

  // warning BFGS use a full matrix of size nxn (where n=u.n) 
  b=1; u=2; // set  right hand side and initial gest
  BFGS(J,DJ,u,eps=1.e-6,nbiter=20,nbiterline=20);
   cout << "BFGS: J(u) = " << J(u) << " err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);

  
};
{ // ---  a real non linear test ---
mesh Th=square(10,10);  // mesh definition of $\Omega$
fespace Vh(Th,P1);      // finite element space
fespace Ph(Th,P0);      // make optimization

Vh b=1;  // to defined b 
// $ J(u) = 1/2\int_\Omega f(|\nabla u|^2) - \int\Omega  u b $
// $ f(u) = a*u + u-ln(1+u), \quad f'(u) = a+\frac{u}{1+u}, \quad f''(u) =  \frac{1}{(1+u)^2}$
real a=0.001;
func real f(real u) { return u*a+u-log(1+u); }
func real df(real u) { return a+u/(1+u);}
func real ddf(real u) { return 1/((1+u)*(1+u));}

// the functionnal J 

func real J(real[int] & u)
  {
    Vh w;w[]=u; 
    real r=int2d(Th)(0.5*f( dx(w)*dx(w) + dy(w)*dy(w) ) - b*w) ;
    cout << "J(u) =" << r << " " << u.min <<  " " << u.max << endl;
    return r;
  }
// -----------------------

Vh u=0; //  the current value of the solution
Ph alpha; // of store  $df(|\nabla u|^2)$
int iter=0;
alpha=df( dx(u)*dx(u) + dy(u)*dy(u) ); // optimization 

func real[int] dJ(real[int] & u)
  {
    int verb=verbosity; verbosity=0; 
    Vh w;w[]=u; 
    alpha=df( dx(w)*dx(w) + dy(w)*dy(w) ); // optimization 
    varf au(uh,vh) = int2d(Th)( alpha*( dx(w)*dx(vh) + dy(w)*dy(vh) ) - b*vh)
    + on(1,2,3,4,uh=0);
    u= au(0,Vh);  
    verbosity=verb;
    return u; // warning no return of local array  
  }

varf alap(uh,vh)=  
   int2d(Th)( alpha *( dx(uh)*dx(vh) + dy(uh)*dy(vh) ))   + on(1,2,3,4,uh=0);

varf amass(uh,vh)=  int2d(Th)( uh*vh)  + on(1,2,3,4,uh=0);

matrix Amass = amass(Vh,Vh,solver=CG);

matrix Alap=  alap(Vh,Vh,solver=Cholesky,factorize=1);   

func real[int] C(real[int] & u)
{
   real[int] w=Amass*u;
   u = Alap^-1*w;
   return u; // no return of local array  variable 
}
   verbosity=5;
   int conv=0;
   real eps=1e-6; 
   for(int i=0;i<20;i++)
   {
     conv=NLCG(dJ,u[],nbiter=10,precon=C,veps=eps); 

     if (conv) break; 
     alpha=df( dx(u)*dx(u) + dy(u)*dy(u) ); // optimization 
     Alap = alap(Vh,Vh,solver=Cholesky,factorize=1);   
     cout << " restart with new preconditionner " << conv << " eps =" << eps << endl;
   }
   plot (u,wait=1,cmm="solution with NLCG");
   umax = u[].max; 

   Vh sss= df( dx(u)*dx(u) + dy(u)*dy(u) ) ;
   plot (sss,wait=0,fill=1,value=1);

// the  method of  Newton Ralphson to solve dJ(u)=0;
//  see Newton.edp example

}