NavierRK4.cc

00001 //
00002 // Chorin Projection Method
00003 // with CONSISTENT pressure Laplacian
00004 // for 2D or 3D Navier-Stokes equations
00005 //
00006 // periodic BCs only! 
00007 //
00008 
00009 # include "Wavelet.hpp"
00010 
00011 # include "AdaptiveData.hpp"
00012 # include "Solver.hpp"
00013 
00014 int debugRefine ;
00015 
00016 # define D 2
00017 
00018 int main(int argc,char **argv) {
00019  int i,j,k ;
00020  // basis functions 
00021  Wavelets WC("Interpolet",6) ;
00022  bool     useLifting=false ;
00023 
00024  // problem parameters
00025  char   id[200] ;
00026  double nu,dt,maxdt,T=400.0,eps0,epsfak=1.3,cfl=0.3,printdt=0.01,hmin,printt,starttime=0,q1=0,qinf=0,starteps=-1.0 ;
00027 
00028  int    prstep=10,ts=0,largestlevel[3],dof,printn,refinestep=25,EZstep=0 ; // output of U[],P each <prstep> time step
00029 
00030  // domain
00031  int BC[D][2],MaxLevel,MaxDOF=3000000 ;
00032  for (i=0; i<D; i++) { BC[i][0]=-1, BC[i][1]=PERIODIC ;}
00033 
00034  // read arguments
00035  if (argc<7) {
00036    printf("NavierRK4 (%dD) -ID   <problem>   // AWFD/Data/<problem> is the directory where the initial condition file UVWP0\n",D);
00037    printf("                                  // is expected and in which the current velocity files are written\n") ; 
00038    printf("          -nu   <nu>              // The viscosity\n") ;
00039    printf("          -starttime  <start> (0) //\n") ;
00040    printf("          -T    <T_end>   (400.0) // The time goes from <starttime> till <T_end> \n") ;
00041    printf("          -dt   <max-dt>          // Maximal time step. However smaller time steps may be used to meet the CFL condition.\n")              ;
00042    printf("          -CFL  <cfl>       (0.3) // Maximal allowed CFL number\n")   ;
00043    printf("          -rs   <refine-step> (25)// Refine the grid each <refine-step>-th time step\n")    ;
00044    printf("          -ML   <MaxLevel>        // Maximal refinement level. Limiting is useful to avoid very small time steps.\n") ;
00045    printf("          -DOF  <max-DOF>   (3e+6)// Maximal number of degrees of freedom. \n") ;
00046    printf("                                  // Useful to avoid 'out of memory' error in long running simulations.\n") ;
00047    printf("          -eps  <threshold>       // (Minimal) threshold parameter\n") ;
00048    printf("          -starteps <threshold>   // If the actual number of DOF is close to <max-DOF> the threshold-parameter is adjusted,\n") ;
00049    printf("                                  // i.e. multiplied by a factor (>1) to relax the refinement criterion. If the actual number\n");
00050    printf("                                  // of DOF then gets significantly  smaller than <max-DOF>, the current threshold is divided\n") ;
00051    printf("                                  // by this factor, until is reaches <threshold>. <starteps> is just the initial value of the\n") ;
00052    printf("                                  // current threshold parameter. Useful for restarts of a simulation.\n");
00053    printf("          -epsfactor <adjustment for epsilon> (1.3) // The above factor.\n") ;
00054    printf("          -q1   <weight 2^{q1  |l|_1  } for coefficients>  (0) // Parameters for level-dependent weights for the threshold,\n") ;
00055    printf("          -qinf <weight 2^{qinf|l|_inf} for coefficients>  (0) // see Refine()\n") ;
00056    printf("          -printdt <print-timestep> (0.01) // Write current velocity in regular time intervals. The files are numbered. \n") ;
00057    printf("                                           // This is the size of these intervals.\n");
00058    printf("          -prstep     (10)        // Write current velocity each <refine-step>*<prstep>-th timestep in the file UVWP.\n");
00059    printf("          -EZstep     (0)         // If >0 compute and print enstrophy and energy each <EZstep>-th time step.\n");
00060    exit(-1) ;
00061  }
00062 
00063  // problem / nu / dt 
00064  for (i=1; i<argc; i++) {
00065    if (strcmp(argv[i],"-ID")==0) sscanf(argv[++i],"%s" ,id ) ;    
00066    else
00067      if (strcmp(argv[i],"-nu")==0) sscanf(argv[++i],"%le",&nu) ;
00068    else 
00069      if (strcmp(argv[i],"-T")==0) sscanf(argv[++i],"%le",&T) ;
00070    else 
00071      if (strcmp(argv[i],"-dt")==0) sscanf(argv[++i],"%le",&maxdt) ;
00072    else 
00073      if (strcmp(argv[i],"-CFL")==0) sscanf(argv[++i],"%le",&cfl) ;
00074    else 
00075      if (strcmp(argv[i],"-ML")==0) sscanf(argv[++i],"%d" ,&MaxLevel) ;
00076    else 
00077      if (strcmp(argv[i],"-eps")==0) sscanf(argv[++i],"%le",&eps0) ;
00078    else 
00079      if (strcmp(argv[i],"-starteps")==0) sscanf(argv[++i],"%le",&starteps) ;
00080    else 
00081      if (strcmp(argv[i],"-rs")==0) sscanf(argv[++i],"%d",&refinestep) ;
00082    else 
00083      if (strcmp(argv[i],"-DOF")==0) sscanf(argv[++i],"%d",&MaxDOF) ;
00084    else 
00085      if (strcmp(argv[i],"-epsfactor")==0) sscanf(argv[++i],"%le",&epsfak) ;
00086    else 
00087      if (strcmp(argv[i],"-q1")==0) sscanf(argv[++i],"%le",&q1) ;
00088    else 
00089      if (strcmp(argv[i],"-qinf")==0) sscanf(argv[++i],"%le",&qinf) ;
00090    else 
00091      if (strcmp(argv[i],"-printdt")==0) sscanf(argv[++i],"%le",&printdt) ;
00092    else 
00093      if (strcmp(argv[i],"-prstep")==0) sscanf(argv[++i],"%d",&prstep) ;
00094    else 
00095      if (strcmp(argv[i],"-EZstep")==0) sscanf(argv[++i],"%d",&EZstep) ;
00096    else 
00097      if (strcmp(argv[i],"-starttime")==0) sscanf(argv[++i],"%le",&starttime) ;
00098    else { printf("syntax-error in command line: %s\n",argv[i]); exit(-1) ; }
00099  }
00100 
00101  dt=maxdt ;
00102  printf("Problem: %s  with dimension %d\n",id,D) ; 
00103  printf("viscosity nu: %e\n",nu) ;
00104  printf("starttime: %e endtime: %e     max. timestep: %e\n",starttime,T,maxdt) ;
00105 
00106  // adaptivity
00107  if (starteps<0) starteps=eps0 ; 
00108  AdaptivityCriterion R(AdaptivityCriterion::Hq_THRESHOLD,starteps) ;
00109  R.MaxLevel=MaxLevel ;
00110  R.q1  =q1   ;
00111  R.qinf=qinf ;
00112  
00113  if ((R.q1!= 0) && (R.qinf!= 0)) printf("threshold criterion:  |u_(l,t)|2^{|l|_1*%e + |l|_inf*%e} >= %e  ,  MaxLevel: %d  MaxDOF: %d   epsfaktor:%e\n",
00114     R.q1,R.qinf,R.eps,R.MaxLevel, MaxDOF, epsfak) ;
00115 
00116  if ((R.q1!= 0) && (R.qinf== 0)) printf("threshold criterion:  |u_(l,t)|2^{|l|_1*%e} >= %e  ,  MaxLevel: %d  MaxDOF: %d   epsfaktor:%e\n",
00117     R.q1,R.eps,R.MaxLevel, MaxDOF, epsfak) ;
00118 
00119  if ((R.q1== 0) && (R.qinf!= 0)) printf("threshold criterion:  |u_(l,t)|2^{|l|_inf*%e} >= %e  ,  MaxLevel: %d  MaxDOF: %d   epsfaktor:%e\n",
00120     R.qinf,R.eps,R.MaxLevel, MaxDOF, epsfak) ;
00121 
00122  if ((R.q1== 0) && (R.qinf== 0)) printf("threshold criterion:  |u_(l,t)| >= %e  ,  MaxLevel: %d  MaxDOF: %d   epsfaktor:%e\n",
00123     R.eps,R.MaxLevel, MaxDOF, epsfak) ;
00124 
00125  // io
00126  printf("output each %e seconds  and  each %dth timestep, EZstep: %d\n",printdt,prstep,EZstep) ;
00127 
00128  // other
00129  int    which[4]={0,1,2,3} ;
00130  double t,umax ;
00131 
00132  // grid alloc,  Level0 must be known
00133  AdaptiveGrid<D> G(&WC) ;
00134  AdaptiveData<D> P(&G) , DP(&G) , U[D] , UN[D] , RHS[D] ,
00135         *Tmp    =G.tmp[AdaptiveGrid<D>::NUM_TMP-1] ,
00136         *Tmp1   =G.tmp[AdaptiveGrid<D>::NUM_TMP-2] ,
00137         *Tmp2   =G.tmp[AdaptiveGrid<D>::NUM_TMP-3] ,
00138         *K[4]   ={G.tmp[AdaptiveGrid<D>::NUM_TMP-4], G.tmp[AdaptiveGrid<D>::NUM_TMP-5], G.tmp[AdaptiveGrid<D>::NUM_TMP-6], G.tmp[AdaptiveGrid<D>::NUM_TMP-7] },
00139         *V      =G.tmp[AdaptiveGrid<D>::NUM_TMP-8] ,
00140         *UVWP[D+1] ;
00141 
00142  UVWP[D]=&P ;
00143  for (i=0; i<D; i++) {
00144    UVWP[i]=&U[i] ;
00145    U  [i].Attach(&G) ;
00146    UN [i].Attach(&G) ;
00147    RHS[i].Attach(&G) ;
00148 
00149    U[i].SetRefine() ;
00150  
00151     U[i].SetBoundaryConditions(BC) ;
00152    UN[i].SetBoundaryConditions(BC) ;
00153  }
00154  
00155   P.SetRefine() ;
00156  DP.SetRefine() ;
00157  
00158   P.SetBoundaryConditions(BC) ;
00159  DP.SetBoundaryConditions(BC) ;
00160 
00161  // read initial conditions
00162  char name[1000] ;
00163  sprintf(name,"%s/Data/%s/UVWP0",_WAVELET_ROOT_,id) ;
00164  G.ReadSparse(name) ;
00165  for (i=0; i<D; i++)
00166    printf("XA/E[%d]=%e %e\n",i,G.XA[i],G.XE[i]) ;
00167 
00168  U[0].ReadSparse(name,UVWP,D+1,which) ;
00169  for (i=0; i<D; i++) {
00170    printf("|U[%d]| %e\n",i,U[i].MaxAbs()) ;
00171    U[i].SetBoundaryConditions(BC) ;
00172  }
00173 
00174  // refine Grid for the first time
00175  j=G.Size() ;
00176 
00177  for (i=0; i<D; i++) {
00178    if (i==0) Tmp->Abs(&U[i]) ;
00179    else {
00180      RHS[0].Abs(&U[i]) ;
00181      Tmp->PPlus(&RHS[0]) ;
00182     }
00183   }
00184  G.Refine(Tmp,&R)  ;
00185  G.LargestActiveLevel(largestlevel) ;
00186  dof=G.Size() ;
00187 
00188  printf("DOF before Refine %d  after %d, active Levels: %d %d %d\n",j,dof,largestlevel[0], largestlevel[1], largestlevel[2]) ;
00189 
00190  for (i=0; i <D; i++) printf("|U[%d]| %e\n",i,U[i].MaxAbs()) ;
00191 
00192  DP.Set(0.0) ;
00193 
00194  // initialize operators 
00195  PressureOperatorIII<AdaptiveData<D>,D> POISSON ;
00196  HelmholtzOperator  <AdaptiveData<D>,D> DIFF    ;
00197  ConvectionOperator <AdaptiveData<D>,D> CONV    ;
00198 
00199  // for Laplace(P)
00200  POISSON.Tmp[0]=Tmp  ;
00201  POISSON.Tmp[1]=Tmp1 ;
00202  POISSON.do_restrict=false ;
00203  POISSON.do_patch   =false ;
00204  POISSON.use_lifting_preconditioner=useLifting ;
00205 
00206  // (1/dt -nu*Lap)V
00207  DIFF.Tmp =Tmp  ;
00208  DIFF.Tmp2=Tmp1 ;
00209  DIFF.par[0]=0.0 ;
00210  for (i=0; i<D; i++) DIFF.par[i+1]=-nu ;
00211 
00212  // K(U,.)X
00213  CONV.Tmp =Tmp    ;
00214  CONV.Tmp1=Tmp1   ;
00215  for (i=0; i<D; i++) CONV.U [i]=&UN[i] ;
00216 
00217  Solver<AdaptiveData<D>,D> S ;
00218  S.eps          =1e-6  ;
00219  S.maxiter      =100   ;
00220  S.Op           =NULL  ;
00221  S.StopCriterion=StoppingCriterionAbsoluteResidual ;
00222  S.print        =NULL  ;
00223  S.prstep       =1     ;
00224  for (i=0; i<Solver<AdaptiveData<D>,D>::NUMTMP_BICGSTAB2; i++) S.Tmp[i]=G.tmp[i] ;
00225 
00226  //
00227  // main loop
00228  //
00229 
00230  t     =starttime ;
00231  printt=t ;
00232  printn=(int)(starttime/printdt+0.9) ;
00233 
00234  while (t<T) {
00235    int advcgit=0 , poscgit=0 ;
00236 
00237    // output
00238    if (t>=printt) {
00239      sprintf(name,"%s/Data/%s/UVWP.%d",_WAVELET_ROOT_,id,printn) ;
00240      U[0].WriteSparse(name,UVWP,D+1,true) ;
00241      printn++ ;
00242      printt += printdt ;
00243    }
00244 
00245    // energy  and enstrophy
00246    if (EZstep>0)
00247      if ( (ts%EZstep) == 0) {
00248        double eng=0, ens=0 ;
00249        for (i=0; i<D; i++) eng += sqr(U[i].L2Norm(Tmp)) ; 
00250 
00251        int ia= (D==2) ? 2 : 0 ;
00252 
00253        for (i=ia ; i<=2; i++) {
00254          j=(i+1)%3 , k=(i+2)%3 ;
00255          Tmp1->ApplyOp(BC[j],&U[k] , j , FIRSTDERIVATIVE) ;
00256          Tmp ->ApplyOp(BC[k],&U[j] , k , FIRSTDERIVATIVE) ;
00257          Tmp1->Sub(Tmp1 ,  Tmp) ;
00258          ens += sqr(Tmp1->L2Norm(Tmp)) ; 
00259         } 
00260        printf("EZ %e : %e %e\n",t,eng,ens/2) ;
00261      }
00262 
00263    //
00264    // transform U[i] to nodal values and calculation of maximal velocity
00265    //
00266    umax=0.0 ;
00267    for (i=0; i<D; i++) {
00268      for (j=0; j<D; j++) 
00269        UN[i].ApplyOp( (j==0) ? &U[i] : &UN[i] , j , INVERSETRANSFORM) ;
00270 
00271      if (i==0) V->Sqr(&UN[0]) ;
00272      else { 
00273        K[0]->Sqr(&UN[0]) ; 
00274        V->PPlus(K[0])    ;
00275       }
00276     }
00277    umax=sqrt(V->Max()) ;
00278 
00279    // time-step size induced by CFL condition
00280    hmin= 1.0/(1<<max(largestlevel[0],max(largestlevel[1],largestlevel[2]))) ;
00281    dt  =cfl*min(maxdt , hmin/umax) ;
00282 
00283    // Tilde(U)
00284    for (i=0; i<D; i++) {
00285 
00286      // Runge-Kutta 4 scheme
00287      CONV.ApplyConservativeWENO(&U[i]  , K[0]) ;
00288      //CONV.ApplyConservative(&U[i],K[0]) ; worse results for Molenkamp and 2D mixing layer
00289 
00290      V->Add(1.0,&U[i], -0.5*dt, K[0]) ;
00291      CONV.ApplyConservativeWENO(V  , K[1]) ;
00292      //CONV.ApplyConservative(V,K[1]) ;worse results for Molenkamp and  2D mixing layer
00293 
00294      V->Add(1.0,&U[i], -0.5*dt, K[1]) ;
00295      CONV.ApplyConservativeWENO(V  , K[2]) ; 
00296      //CONV.ApplyConservative(V,K[2]) ;worse results for Molenkamp and 2D mixing layer 
00297 
00298      V->Add(1.0,&U[i],     -dt, K[2]) ;
00299      CONV.ApplyConservativeWENO(V  , K[3]) ;
00300      //CONV.ApplyConservative(V,K[3]) ;worse results for Molenkamp and 2D mixing layer
00301 
00302      Tmp2->Add(1.0/6,K[0],  2.0/6,K[1],  2.0/6,K[2],  1.0/6,K[3]) ;
00303 
00304      // Tmp   = projection of grad p
00305      Tmp->ApplyOp(BC[i] , &P , i , GRAD) ;
00306      // Rhs =   - K(U,U)        -grad P       + U/dt
00307      RHS[i].Add(-1.0,Tmp2 , -1/dt,Tmp,  1/dt,&U[i]) ;
00308 
00309      S.print=NULL ;
00310      if (/*(i==0)&&*/((ts<100)||((ts%10)==0))) S.print="ADVCG:" ;
00311      DIFF.par[0]=1./dt ;
00312      S.Op=&DIFF ;
00313      //S.maxiter      =60 ;
00314      //if (ts==5) debugRefine =1 ;
00315      advcgit +=S.BiCGSTAB(&U[i],&RHS[i]) ; 
00316    }
00317 
00318    // divTilde(U)
00319    for (i=0; i<D; i++) {
00320      RHS[i].ApplyOp(BC[i] , &U[i] , i , DIV) ;
00321      if (i>0) RHS[0].PPlus(&RHS[i]) ;
00322    }
00323   
00324    S.print=NULL ;
00325    if ((ts<100)||((ts%10)==0)) S.print="PSCG:" ;
00326    S.Op=&POISSON ;
00327    if (ts>100) S.maxiter=20 ;
00328    //S.maxiter=20 ;
00329    poscgit =S.BiCGSTAB2(&DP,&RHS[0]) ;
00330 
00331    P.PPlus(&DP) ;
00332 
00333    // U=Tilde(U)-grad(DP)
00334    for (i=0; i<D; i++) {
00335      Tmp->ApplyOp(BC[i] , &DP , i ,  GRAD) ;
00336      U[i].MMinus(Tmp) ;
00337    }
00338 
00339 // divergence free ?
00340 for (i=0; i<D; i++) {
00341   RHS[i].ApplyOp(BC[i] , &U[i] , i , DIV) ;
00342   if (i>0) RHS[0].PPlus( &RHS[i] ) ;
00343 }
00344 printf ("divergence: %e  %e\n",RHS[0].MaxAbs() , RHS[0].InnerProd(&RHS[0])) ;
00345 
00346 
00347    printf("loop %d  %e  %d  %d umax: %e\n",ts,t, advcgit , poscgit , umax) ;
00348    fflush(stdout) ;
00349   
00350    // refine grid
00351    if ((ts%refinestep)==0) {
00352 
00353      if ( ((ts/refinestep)%prstep) == 0 ) U[0].WriteSparse("UVWP",UVWP,D+1,true) ;
00354 
00355      for (i=0; i<D; i++) {
00356        if (i==0) Tmp->Abs(&U[i]) ;
00357        else {
00358          RHS[0].Abs(&U[i]) ;
00359          Tmp->PPlus(&RHS[0]) ;
00360         }
00361       }
00362 
00363      //  save error indicators
00364      Tmp1->Copy(Tmp) ;
00365      Tmp1->SetRefine() ;
00366      
00367      // refine with adjustment of threshold value
00368      bool repeat ;
00369      do {
00370        repeat=false ;
00371        G.Refine(Tmp,&R) ;
00372        
00373            G.LargestActiveLevel(largestlevel) ;
00374        dof=G.Size() ;
00375        printf("eps=%e | DOF %d | largest active level %d %d %d\n",R.eps,dof,largestlevel[0],largestlevel[1],largestlevel[2]) ;
00376  
00377        if (dof >     MaxDOF) { R.eps *= epsfak                   ; repeat=true ; Tmp->Copy(Tmp1) ; }
00378        if (dof < 0.9*MaxDOF) { R.eps  = max(eps0 , R.eps/epsfak) ; }
00379      } while (repeat) ;
00380 
00381    }
00382    
00383    ts++ ;
00384    t += dt ;
00385 
00386   }
00387 }
Generated at Mon Aug 19 10:02:32 2002 for AWFD by