ELBM: 熵格子Alpha求解C++程序

Posted by on 05/5/2012 in Generalized LBM, LBM理论 | 0 comments

这分享我使用的求解ELBM中非线性Alpha方程的核心代码,希望这对感兴趣熵格子方法的朋友有用!欢迎大家使用!
//-----------------------------------------------------------
#include <math.h>
//-----------------------------------------------------------
//External variables used in the computation.
//-----------------------------------------------------------
#define TOL            5e-16
#define   MAXIT        200
#define   TOLL          1e-7
//
double     log_[9]={1./8,0.5,0.5,0.5,0.5,2,2,2,2};
double     F_eval_r(double);
void         FD_eval_r(double,double *,double *);
double     NewtonMod(void);
double     inf_, sup_;
double     Finf_, Fsup_;
double     sum_f;
//-----------------------------------------------------------
double    alpha_(void)
{
   //-------------------------------------------------------
   Finf_=(*F_eval)(1.6);
   Fsup_=(*F_eval)(2.2);
   double    val_=Finf_*Fsup_;
   if(val_>=0)
       return    2;
   //-------------------------------------------------------
   return NewtonMod();
}
//-----------------------------------------------------------
//    Function that evaluate the following function:
//        F(alfa)=H(f_i)-H(f_i+alfa*(f_i-f_i^eq))
//-----------------------------------------------------------
double    F_eval_r(double    x)
{
   double    fn_loc=0;
   double    Delta_i=0, F_alfa=0;
   //-------------------------------------------------------
   //        H(f+alfa*Delta) - H(f) = 0
   //-------------------------------------------------------
   double    app_log=0;
   sum_f=0;
   for(int k=0;k<9;k++)
   {
       Delta_i=(equilib.link[k]-node_loc.link[k]);
       F_alfa=(node_loc.link[k]+x*Delta_i);
       app_log=(log(log_[k]*fabs(F_alfa)));
       //---------------------------------------------------
       sum_f=sum_f
           +node_loc.link[k]*(log(log_[k]*node_loc.link[k]));
       fn_loc=fn_loc+(F_alfa*app_log);
   }
   fn_loc=fn_loc-sum_f;
   return    fn_loc;
}
//-----------------------------------------------------------
//    Function that evaluate the following function:  
//            F(alfa)=H(f_i)-H(f_i+alfa*(f_i-f_i^eq))
//    and its derivative.
//-----------------------------------------------------------
void    FD_eval_r(double    x, double * fn, double * df)
{
   double    fn_loc=0, df_loc=0;
   double    Delta_i=0, F_alfa=0;
   //-------------------------------------------------------
   //    H(f+alfa*Delta)-H(f) = 0
   //-------------------------------------------------------
   double    app_log=0;
   for(int k=0;k<9;k++)
   {
       Delta_i=(equilib.link[k]-node_loc.link[k]);
       F_alfa=(node_loc.link[k]+x*Delta_i);
       app_log=(log(log_[k]*fabs(F_alfa)));
       //---------------------------------------------------
//        fn_loc=fn_loc+(F_alfa*app_log-calcolo.F_i);
       fn_loc=fn_loc+(F_alfa*app_log);
//        df_loc=df_loc+Delta_i*app_log;         //  (1)
       df_loc=df_loc+(Delta_i*(app_log+1)); //  (2)
   }
   fn_loc=fn_loc-cmpt.sum_f;//
   *fn=fn_loc;
   *df=df_loc;
}
//-----------------------------------------------------------
//
//-----------------------------------------------------------
double NewtonMod(void)
{
   int j;
   double        df,dx,dxold,f;
   double        temp,xh,xl,rts;
   double        x1=1.6, x2=2.2;
   //-------------------------------------------------------
   //ATTENTION: the values Finf_ and Fsup_ (that are global variables)  
   //         have been evaluated in function alpha_().
   if (Finf_==0.0)  
       return x1;
   if (Fsup_==0.0)  
       return x2;
   //-------------------------------------------------------
   if(Fsup_> 0.0)
   {
       //Orient the search so that f(xl) < 0.
       xl=x1;
       xh=x2;
   }
   else
   {
       xh=x1;
       xl=x2;
   }
   //-------------------------------------------------------
   rts=2; //this is the starting point for the root finder
        //the optimization strategies have to be applied here
   (*FD_eval)(rts,&f,&df);
   if(fabs(f)<=TOL) return rts;
   //
   dxold=fabs(x2-x1);    //(1)
//    dxold=f/df;            //(2)
   dx=dxold;                //and the last step.
   //-------------------------------------------------------
   for (j=1;j<=MAXIT;j++)
   {
       //---------------------------------------------------
       //Loop over allowed iterations.
       if((((rts-xh)*df-f)*((rts-xl)*df-f)>0.0)||(fabs(2.0*f)>fabs(dxold*df)))
       {
           //Bisect if Newton out of range,or not decreasing fast enough.
           dxold=dx;
           dx=0.5*(xh-xl);
           rts=xl+dx;
           if (xl == rts) return rts; //Change in root is negligible.
       }
       else
       {
           //Newton step acceptable. Take it.
           dxold=dx;
           dx=f/df;
           temp=rts;
           rts -= dx;
           if (temp == rts) return rts;
       }
       //---------------------------------------------------
       if (fabs(dx)<TOLL) return rts; //Convergence criterion.
       //---------------------------------------------------
       (*FD_eval)(rts,&f,&df);
       //---------------------------------------------------
       //The one new function evaluation per iteration.
       if (f < 0.0) //Maintain the bracket on the root.
           xl=rts;
       else
           xh=rts;
   }
   printf("Error. We exceed the maximum number of iterations!\n");
   return 0.0;
}
//-----------------------------------------------------------

 

    3,275 views

Submit a Comment