Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

CP.cc

//
// solves the linear convection problem    du/dt + nabla*(Au) = 0
// u is the advected scalar and A(x,y)=(-y,x) is the given velocity (a rotation in our case)
// 
// An explicit multistep scheme (3rd order Adams-Bashforth)
// and a WENO discretization of the convection term are used. 
//
// At regular intervals, the current solution is evaluated on a uniform grid and written to a 
// file with advancing time mark. These files can be read and visualized using MATLAB
// 

# include "AdaptiveData.hpp"
# include "Operators.hpp"
# include "Solver.hpp"

int debugRefine ; // needed for dummy reasons

struct HillSharpFunction : public Function {
  double x0[3] ;
  double sigma ;
  void   SetParameters(double *p) { int    i  ; for (    i=0; i<2; i++) x0[i]=p[i]; sigma=p[i]; }
  double Eval         (double *x) { double r=0; for (int i=0; i<2; i++) r+=pow(x[i]-x0[i],2)  ; return max(4*exp(sigma*r)-3.0 , 0.0); }
} ;

int main(int argc, char **argv) {

  Wavelets WC("Interpolet",4) ; // use 4th order Interpolets

  int    i,prstep,dof  ; 
  int    BC[2][2],BG[2][2] ; // boundary conditions markers 
  for    (i=0; i<2; i++) BC[i][0]=-1 , BC[i][1]=PERIODIC ;
  
  int    Prints=100 ;
  double dt=3.9e-4, T=2*PI,t=0 ; // time step, stop time 
  
  prstep=(int)(T/(Prints*dt)+0.9999) ;
  dt=T/(prstep*Prints) ;  // adjust dt such that T is exactly hit.
  
  // refinement criterion
  AdaptivityCriterion R(AdaptivityCriterion::L2_THRESHOLD,1e-3) ;
  R.MaxLevel=8 ; 
  assert(dt<0.1/(1<<R.MaxLevel)) ;
  
  char   name[300]  ;
  
  // set functions for velocity
  double pars[20] ;
  TensorNPolynomFunction<2> fv[2] ;
  
  i=0;            // u(x,y)= 1.0*(1  )*(0 -y)
  pars[i++]=1.0 ;
  pars[i++]= 0  ; // degree          of x-polynomial 
  pars[i++]= 1  ; // degree          of y-    "
  pars[i++]= 1  ; // 0th coefficient of x-    "
  pars[i++]= 0  ; // 0th coefficient of y-    "
  pars[i++]=-1  ; // 1th coefficient of y-    "
  fv[0].SetParameters(pars) ;
  
  i=0;            // v(x,y)= 1.0*(0+x)*(1)
  pars[i++]=1.0 ;
  pars[i++]= 1  ; // degree          of x-polynomial 
  pars[i++]= 0  ; // degree          of y-    "
  pars[i++]= 0  ; // 0th coefficient of x-    "
  pars[i++]= 1  ; // 1th coefficient of x-    "
  pars[i++]= 1  ; // 0th coefficient of y-    "
  fv[1].SetParameters(pars) ;
  
  // initial condition is a sharp edged hill
  HillSharpFunction fu0   ;
  pars[0]=1./4  ;
  pars[1]=1./4  ;
  pars[2]=-30/2 ;
  fu0.SetParameters(pars) ;
  
  // AdaptiveData
  AdaptiveGrid<2> G(&WC) ;
  G.XA[0]=G.XA[1]=-1 ;
  G.XE[0]=G.XE[1]= 1 ;
  
  AdaptiveData<2> U[4] , UN[2] , K(&G) , Tmp(&G) , Tmp1(&G), Tmp2(&G), *Uold,*Unew ;
  
  // BoundaryConditions
  for (i=0; i<2; i++) GetGeneralBoundaryConditions(BC[i] , BG[i]) ;
  for (i=0; i<2; i++) {
    UN[i].Attach(&G) ;
    UN[i].SetBoundaryConditions(BG) ;
  }
  for (i=0; i<4; i++) {
    U [i].Attach(&G) ;
    U [i].SetBoundaryConditions(BC) ;
    U [i].SetRefine() ;
  }
  
  K.SetBoundaryConditions(BC) ;
  K.SetRefine() ;
  
  // initial Grid
  G.SetFull(7)     ;
  
  // initial Conditions
  U[0].SetFunction(&fu0) ;
  for (i=0; i<2; i++) U[0].ApplyOp(BC[i] , &U[0] , i ,PROJECTION) ;
  Tmp.Copy(&U[0]) ;
  
  G.Refine(&Tmp,&R) ;
  dof=G.Size() ;
  
  for (i=0; i<2; i++) UN[i].SetFunction(&fv[i], false, true) ;
  
  // operators
  ConvectionOperator<AdaptiveData<2>,2> CONV ;
  CONV.Tmp =&Tmp      ;
  CONV.Tmp1=&Tmp1     ;
  for (i=0; i<2; i++) CONV.U [i]=&UN[i] ;
  
  //
  // Main Loop
  //
  int tsk=0, ts=0  ;
  Unew=&U[0] ;
  
  // T=dt ;
  while (t<T+dt) {
    
    // Output
    if (ts%prstep==0) {
      i=ts/prstep ;
      sprintf(name,"../../Data/Test/U.%d", i) ;
      Unew->WriteUDF(name,7,&Tmp2) ;
    }
    
    //AB3
    Uold=&U[0] ;
    Unew=&U[0] ;
    AdaptiveData<2> *KK[3]={&U[1],&U[2],&U[3]} ;
    CONV.ApplyConservativeWENO(Uold , KK[tsk%3]) ;
    
    switch (tsk) {
    case 0 :
      Tmp.Copy(KK[tsk%3]) ;
      break ;
    case 1 :
      Tmp.Add(2.0,KK[tsk%3] , -1.0,KK[(tsk+2)%3]) ;
      break ;
    default : 
      Tmp.Add(23./12,KK[tsk%3]  , -16./12,KK[(tsk+2)%3] , 5./12,KK[(tsk+1)%3]) ;
      break ;
    }
    
    Unew->Add(1.0,Uold , -dt,&Tmp) ;
    
    for (i=0; i<2; i++) Tmp.ApplyOp(Unew->Ext.BC[i] , (i==0) ? Unew : &Tmp, i , INVERSETRANSFORM) ;
    if ((ts%10)==0) printf("ts= %4d   t= %e   max U= %e   DOF= %d\n",ts ,t , Tmp.MaxAbs() , dof) ;
    fflush(stdout) ;
    
    // refine grid
    if ((ts%10==0) && (R.strategy!=AdaptivityCriterion::FIXED_LEVEL)) {
      Tmp.Copy(Unew) ;
     
      G.Refine(&Tmp,&R) ;
      dof=G.Size() ;
      
      // evaluate velocity functions on the new adaptive grid
      for (i=0; i<2; i++) UN[i].SetFunction(&fv[i], false, true) ;
    }
    
    tsk++ ;
    ts++  ;
    t+=dt ;
    
  } // next t
}

---