Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

WENO.cc

#include "WENO.hpp" 
#include "Adaptive1D.hpp"

void IndexSet::VecApplyWENO(double *to, double *v, double *roespeed, int *BC, int J, int k , double DX) {
 assert ((k==2) || (k==3)) ;
 assert (BC[1]==PERIODIC) ;

 int i,nJ=1<<J,r,m,s ;
 double vr[10],br[10],d[10],aa[10],sa,fm,fi,vi,vi1,vi2,vm1,vm2,rs ;
 
 switch (k) {
    case 2: {d[0]=2./3 , d[1]=1./3 ;}
    break ;
    case 3: {d[0]=3./10, d[1]=3./5, d[2]=1./10 ;}
    break ;
  }
  
 // reconstruction of f
 
 m=(1<<J)-1 ;
 for (r=0; r<nr; r++)
   for (i=a[r]; i<=e[r]; i++) {
     // rs =  roespeed[i] ;
     rs = (roespeed[i] + roespeed[(i+1)&m]) / 2  ;

     //
     // calculate v(r)_(i+1/2) and  v(r)_(i-1/2)
     //  
     if (k==2) {
       if (rs >= 0) {
         vm1=v[(i+m)&m], vi=v[i],vi1=v[(i+1)&m] ;
       } else {
         vm1=v[(i+2)&m], vi=v[(i+1)&m],vi1=v[i] ;
       }
       
       vr[0]=(vi+vi1)/2      ;
       vr[1]=-0.5*vm1+1.5*vi ;

       br[0]=sqr(vi1-vi) ;
       br[1]=sqr(vi-vm1) ;
     }

     if (k==3) {
       if (rs >= 0) {
         vm2=v[(i+m-1)&m],  vm1=v[(i+m)&m], vi=v[i],vi1=v[(i+1)&m], vi2=v[(i+2)&m] ;
       } else {
         vm2=v[(i+3)&m],  vm1=v[(i+2)&m], vi=v[(i+1)&m],vi1=v[i], vi2=v[(i+m)&m] ;
       }

       vr[0]=( 2*vi  + 5*vi1 -   vi2)/6 ;
       vr[1]=(  -vm1 + 5*vi  + 2*vi1)/6 ;
       vr[2]=( 2*vm2 - 7*vm1 +11*vi )/6 ;

       br[0]=13./12.*sqr(vi -2*vi1+vi2) + sqr(3*vi -4*vi1+  vi2)/4;
       br[1]=13./12.*sqr(vm1-2*vi +vi1) + sqr(  vm1-vi1        )/4;
       br[2]=13./12.*sqr(vm2-2*vm1+vi ) + sqr(  vm2-4*vm1+3*vi )/4;
     }

     sa=0 ;
     for (s=0; s<k; s++) { 
       aa[s]=d[s]/sqr(1e-6 + br[s]); 
       sa +=aa[s] ;
      }

     to[i]=0 ;  
     for (s=0; s<k; s++) to[i] += (aa[s]*vr[s])/sa ;
   }

 // conservative discretization 
 
 double q=nJ/DX ;
 for (r=0; r<nr; r++) {
   fm = (a[r]==0) ? to[nJ-1] : 0 ;
   for (i=a[r]; i<=e[r]; i++) {
     fi=to[i] ;
     to[i]=(fi-fm)*q ;
     fm=fi ; 
   }
 }

}



/*
void WENO_Apply(double *dv,int * , int *, double *v,int * , Wavelets *,int J,int ) {
 int i,nJ=1<<J,k=3,r ;
 double vs[10],vr[10],br[10],om[10],omt[10],d[10],a[10],at[10],sa,sat,vp,vm ;
 
 double *f=(double *)malloc(sizeof(double)*(nJ+1)) ;

 switch (k) {
    case 2: {d[0]=2./3 , d[1]=1./3 ;}
    break ;
    case 3: {d[0]=3./10, d[1]=3./5, d[2]=1./10 ;}
    break ;
    default: assert(0) ;
  }


 int m=(1<<J)-1 ;
 for (i=0; i<nJ; i++) {

   //
   // calculate v(r)_(i+1/2) and  v(r)_(i-1/2)
   //  
   if (k==2) {
     double vm1=v[(i+m)&m], vi=v[i],vi1=v[(i+1)&m] ;
   
     br[0]=sqr(vi1-vi) ;
     br[1]=sqr(vi-vm1) ;
  
     vr[0]=(vi+vi1)/2      ;
     vr[1]=-0.5*vm1+1.5*vi ;

     vs[0]=1.5*vi - 0.5*vi1 ;
     vs[1]=(vm1+vi)/2       ;
   }

   if (k==3) {
     double  vm2=v[(i+m-1)&m],  vm1=v[(i+m)&m], vi=v[i],vi1=v[(i+1)&m], vi2=v[(i+2)&m] ;
   
     br[0]=13./12.*sqr(vi -2*vi1+vi2) + sqr(3*vi -4*vi1+  vi2)/4;
     br[1]=13./12.*sqr(vm1-2*vi +vi1) + sqr(  vm1-vi1        )/4;
     br[2]=13./12.*sqr(vm2-2*vm1+vi ) + sqr(  vm2-4*vm1+3*vi )/4;

     vr[0]=( 2*vi  + 5*vi1 -   vi2)/6 ;
     vr[1]=(  -vm1 + 5*vi  + 2*vi1)/6 ;
     vr[2]=( 2*vm2 - 7*vm1 +11*vi )/6 ;

     vs[0]=(11*vi  - 7*vi1  + 2*vi2)/6 ;
     vs[1]=( 2*vm1 + 5*vi   -   vi1)/6 ;
     vs[2]=(  -vm2 + 5*vm1  + 2*vi )/6 ;
   }


   sa=0,sat=0 ;
   for (r=0; r<k; r++) {
     a [r]=d[    r]/sqr(1e-6 + br[r]) ;
     at[r]=d[k-1-r]/sqr(1e-6 + br[r]) ;
     sa += a[r] ; 
     sat+=at[r] ;
    }
   
   for (r=0; r<k; r++) om[r]=a[r]/sa, omt[r]=at[r]/sat ;

   vm=vp=0 ;
   for (r=0; r<k; r++) {
     vp += om [r]*vr[r] ;
     vm += omt[r]*vs[r] ;
   }
   
   f[i]=vp ;
 }

 for (i=0; i<nJ; i++) dv[i]=(f[i] - f[(i+m)&m])*nJ ;

 free (f) ;
}
*/

---