Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Navier.hpp

//
// DataStructure for Navier-Stokes Simulations
//

#include <stdio.h>
#include <stdlib.h>

#include "AdaptiveGrid.hpp"
#include "AdaptiveData.hpp"
#include "UniformData.hpp"
#include "Function.hpp"

template<int DIM>
struct Navier {
 
       Navier() ;
      ~Navier() ;

  //
  // Parameters
  char ID[200]     ; // Name of problem required for directory data to be stored in
  int  prstep      ; // output velocity components of timestep % prstep == 0
  int  maxiter     ; // max number of iterations for solver 
  double maxres    ; // max residual ;
  int  J[DIM]      ; // level of refinement for full grid data, i.e. output/input of velocity fields
  
  double dt,T      ; // simulate t\in [0,T] with stepsize dt
  double nu        ; // viscosity
  double gamma     ; // additional parameter

  double umax      ; // max velocity, e.g. timestep size control

  int  BC[DIM][DIM][2] ; // Boundary Conditions for Velocity components u[i]: BC[i][][]

  double XL[DIM] ;
  double OutTimes[1000]   ;
  int         MovieStep   ; // Output of ... each MovieStep-th time step
  int         MovieJ[DIM] ; // Levels for Size of output
  UniformData<DIM> MovieTmp    ; // aux. matrix for output.WriteMatLab(...)

  Function *SetF[DIM],*SetUD[DIM] ;
  
  //
  // always required temp. variables
  int ts,a[DIM],e[DIM],op[DIM],res[DIM],Level0[DIM] ;
  int BG[DIM][2] ;  // 'general' boundary conditions
  double t       ;
  char st[200]   ;

  UniformData<DIM> U0[DIM],P0 ; // for IO

  void WriteParameters() ;
  void ReadParameters (char *NID) ;

  void ReadInitialValues(char *value , char *AddNameExtension , AdaptiveData<DIM> *U) ;

  void SetProblem1() ;
  void SetProblem2() ;
  void SetProblem3() ;
  void SetProblem4() ;
  void SetProblem5() ;
  void SetProblem6() ;
  void SetProblem7() ;
  void SetProblem8() ;
  void SetProblemStokes() ;
  void SetProblemStokes2() ;
  void SetProblemStokes3() ;

  void SetBG() ;

  void Print() ;

  Wavelets WC ;
  void     ReadWavelets() ;
  
} ;

//
// Implementations
//
template<int DIM>
Navier<DIM>::Navier()  { 
 for (int i=0; i<DIM; i++) { 
   XL[i]  =1.0 ;
   SetF[i]=SetUD[i]=NULL ;
   J[i]=0      ;
   MovieJ[i]=0 ;
  }
 MovieStep  =MAXINT ;
 prstep     =MAXINT ;
 OutTimes[0]=1e+100 ;
 T=dt=nu=0.0 ; 
}

template<int DIM>
Navier<DIM>::~Navier() { 
  for (int i=0; i<DIM; i++) { 
    delete SetF[i]; 
    delete SetUD[i] ; 
   }
}

template<int DIM>
void Navier<DIM>::WriteParameters() {
  int i,j ;
  FILE *file ;
  char name[300] ;

  sprintf(name,"%s/Data/%s/parameter.dat",_WAVELET_ROOT_,ID) ;
  if ((file=fopen(name,"w"))==NULL) {
    printf("Could not open file for write: %s\n",name) ; 
    exit(-1) ;
  }

  fprintf(file,"%s\n",ID) ;
  fprintf(file,"%d    /prstep\n",prstep) ;
  fprintf(file,"%d    /maxiter\n",maxiter) ;
  fprintf(file,"%e    /maxres\n",maxres) ;
  for (i=0; i<DIM; i++) 
     fprintf(file,"%d  ",J[i]) ; 
  fprintf(file,"/J[]\n") ;

  fprintf(file,"%e %e  /dt,T\n",dt,T) ;
  fprintf(file,"%e /nu\n",nu) ;
  
  for (i=0;i<DIM;i++) {
    for (j=0; j<DIM; j++) fprintf(file,"%d %d    ",BC[i][j][0],BC[i][j][1]) ; 
    fprintf(file,"\n") ;
  } 
  fclose(file) ;
}

template<int DIM>
void Navier<DIM>::ReadParameters(char *NID) {
  int i,j ;
  FILE *file ;
  char name[300] ;

  sprintf(name,"%s/Data/%s/parameter.dat",_WAVELET_ROOT_,NID) ;
  if ((file=fopen(name,"w"))==NULL) {
    printf("Could not open file for read: %s\n",name) ; 
    exit(-1) ;
  }

  fscanf(file,"%s\n",ID) ;
  fscanf(file,"%d\n",&prstep) ;
  fscanf(file,"%d\n",&maxiter) ;
  fscanf(file,"%le\n",&maxres) ;

  for (i=0; i<DIM; i++) 
     fscanf(file,"%d",&J[i]) ; 
  fscanf(file,"%s\n",name) ;

  fscanf(file,"%le %le\n",&dt,&T) ;
  fscanf(file,"%le\n",&nu) ;
  
  for (i=0;i<DIM;i++) {
    for (j=0; j<DIM; j++) fscanf(file,"%d %d    ",&BC[i][j][0],&BC[i][j][1]) ; 
    fscanf(file,"%s\n",name) ;
  } 
  fclose(file) ;

//
// Init. of some temp. Variables
//
  for (i=0;i<DIM;i++) {
    a[i]=0 ;
    e[i]=(1<<J[i]) ;
   }

  for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;
  
  for (i=0;i<DIM;i++) {
    sprintf(name,"%s/Data/%s/%s",_WAVELET_ROOT_,ID,(i==0) ? "U" : "V") ;
    U0[i].ba->ReadMatLab(name) ;
  }

printf("%d %d %e\n",prstep,maxiter,maxres) ;
for (i=0;i<DIM;i++) printf("J[%d]=%d\n",i,J[i]) ;

printf("%e %e %e\n",dt,T,nu) ;
for (i=0;i<DIM;i++) {
   for (j=0; j<DIM; j++) printf("%d %d    ",BC[i][j][0],BC[i][j][1]) ; 
   printf("\n") ;
 }

}

//
// Output
template<int DIM>
void Navier<DIM>::Print() {
 int i,j ;
 printf("Problem: %s in D=%d\n",ID,DIM) ;
 printf("nu= %e\n",nu) ;
 printf("dt= %e\n",dt) ;
 printf("T = %e\n",T ) ;
 printf("gamma  = %e\n",gamma ) ;
 printf("maxres = %e\n",maxres) ;

 printf ("BC:\n") ;
 for (i=0; i<DIM; i++) {
   for (j=0; j<DIM; j++) printf("[%d,%d] ",BC[i][j][0] , BC[i][j][1]) ;
   printf ("\n") ;
 }
}

//
// Wavelets holen
template<int DIM>
void Navier<DIM>::ReadWavelets() {
  WC.GetCoefficients("Interpolet","left" ,6) ;
  WC.GetCoefficients("Interpolet","right",6) ;
  for (int i=0; i<DIM; i++) Level0[i]=WC.Level0 ;
  WC.LiftingFlag=false ;
}


//
// Set 'general' boundary conditions
template<int DIM>
void Navier<DIM>::SetBG() {
  int i,j ;
  for (j=0; j<DIM; j++) {

    BG[j][0]=BG[j][1] = -1 ;

    if (BC[0][j][1]==WAVEPERIODICBC) {
      BG[j][1]=WAVEPERIODICBC ;

      assert(BC[0][j][0]==-1) ;
      for (i=1; i<DIM; i++) assert((BC[i][j][0]==-1) && (BC[i][j][1]==WAVEPERIODICBC)) ;
    }
  }
}


// Lesen von Anfangswerten + interpolieren auf aktuelles Gitter
template<int DIM>
void Navier<DIM>::ReadInitialValues(char *value , char *ext , AdaptiveData<DIM> *U) {
  Matrix<double,DIM> M ;
  char name[300] ;
  int MJ[DIM],i,a[DIM],e[DIM],ee[DIM] ;
  
  if (ext) sprintf(name,"../TestDat/%s/%s%s",ID,value,ext) ;
     else  sprintf(name,"../TestDat/%s/%s",ID,value) ;

  M.ReadMatLab(name) ;

  M.CalcLevel(MJ) ;
  
  for (i=0;i<DIM;i++) { 
    op[i]=MJ[i]-Level0[i] ;
    printf("ReadInitialValues: %d %d\n",MJ[i] , Level0[i]) ;
    a[i]=0 ;
    e[i]=1<<J[i] ;
    ee[i]=1<<min(J[i],MJ[i]) ;
   }
  
  Matrix<double,DIM> K(a,e) ;
  K.Set(0.0) ;
  M.TensorMatrixOperation(BG , &M , BG , &WC , bwt , op , res) ;

  // copy to right sized matrix
  Matrix<double,DIM>::iterator it(a,ee) ;
  for (it.Set(a); it<=ee; ++it) K[it.i]=M[it.i] ;

  U->FromFull(&K,J) ;
}

  
//
// Some predefined Test-Problems
//
// einfaches Poblem aus Howell/Bell
//
template<int DIM>
void Navier<DIM>::SetProblem1() {
 int i,j ;
 this->ReadWavelets() ;

 // Parameters
 strcpy(ID,"Vortex") ;
 nu=1/15000. ;
 dt=1e-4 ;
 T =1.0  ;
 maxres =1e-5 ;
 maxiter=1000 ;

 prstep=10 ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=6 ;
   for (j=0;j<DIM;j++) BC[i][j][0]=BC[i][j][1]=0 ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }

 this->SetBG() ;

 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 // Berechne 1D Geschwingkeitsprofil u(r) aus gegebener Rotationsverteilung um Ursprung

 double si=0.04 ;

 Matrix<double,DIM>::iterator it(U0[0].ba) ;
 umax=0.0 ;
 for (it=(*U0[0].ba).begin(); it<=(*U0[0].ba).end(); ++it) {

   double x[DIM] , ra ,ur;
   for (i=0; i<DIM; i++) x[i]=((double)it.i[i])/e[i] ;

   ra=sqrt((x[0]-0.4)*(x[0]-0.4) + (x[1]-0.5)*(x[1]-0.5)) ;

   ur=(1-exp(-(ra*ra)/(2*si*si)))/ra ;

   U0[0][it]= -ur*(x[1]-0.5)/ra ;
   U0[1][it]=  ur*(x[0]-0.4)/ra ;
  
   if ((it.i[0]==0)||(it.i[1]==0)||(it.i[0]==e[0])||(it.i[1]==e[1])) U0[0][it]=U0[1][it]=0.0 ;

   umax = (umax < fabs(ur)) ? fabs(ur) : umax ; 
 }
 
 U0[0].SetBoundaryConditions(BC[0]) ;
 U0[1].SetBoundaryConditions(BC[1]) ;

 for (i=0; i<2; i++)
   for (j=0; j<2; j++) U0[i].ApplyOp(BC[i][j], &U0[i], j, TRANSFORM) ;

 double w=0.0 ;
 for (i=0; i<DIM; i++) {
   U0   [i].SetBoundaryConditions(BG) ;
   SetUD[i]=new ConstFunction() ; SetUD[i]->SetParameters(&w) ;
   SetF [i]=new ConstFunction() ; SetF [i]->SetParameters(&w) ;
  }
 
 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

}


//
// Channel-flow zwischen zwei Platten
//
template<int DIM>
void Navier<DIM>::SetProblem2() {
 int i ;
 this->ReadWavelets() ;

 // Parameters
 strcpy(ID,"Channel") ;
 nu=1/100. ;
 dt=1.0e-2 ;
 T =100.0  ; 
 maxres =1e-5 ;
 maxiter=1000 ;

 prstep=10 ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=6 ;
   BC[i][0][0]=-1 ;
   BC[i][0][1]=WAVEPERIODICBC ;

   BC[i][1][0]=0 ;
   BC[i][1][1]=0 ;

   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }

 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 Matrix<double,DIM>::iterator it(U0[0].ba) ;
 umax=0.0 ;
 for (it=(*U0[0].ba).begin(); it<=(*U0[0].ba).end(); ++it) {

   (*U0[0].ba)[it]= 2./3 ;
   (*U0[1].ba)[it]= 0.0 ;
  
   if ((it.i[1]==0)||(it.i[1]==e[1])) U0[0][it]=U0[1][it]=0.0 ;

   umax = 1.0 ;
 }

 double w=0.0 ;
 SetUD[0]=new ConstFunction() ;SetUD[0]->SetParameters(&w) ;
 SetUD[1]=new ConstFunction() ;SetF [1]->SetParameters(&w) ;
 
 SetF [1]=new ConstFunction() ;SetF [1]->SetParameters(&w) ;
 w=8*nu ;
 SetF [0]=new ConstFunction() ;SetF [0]->SetParameters(&w) ;

 this->SetBG() ;
}



//
// TestProblem aus A. Prohl
//
template<int DIM>
void Navier<DIM>::SetProblem3() {
 int i ;
 this->ReadWavelets() ;

 // Parameters
 strcpy(ID,"Prohl") ;
 nu   =1/10.  ;
 dt   =1.0e-3 ;
 T    =1.0    ;
 gamma=1.0    ;

 maxres =1e-5 ;
 maxiter=1000 ;

 prstep=10 ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=7 ;
   BC[i][0][0]=0 ;
   BC[i][0][1]=0 ;

   BC[i][1][0]=0 ;
   BC[i][1][1]=0 ;

   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }

 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 Matrix<double,DIM>::iterator it(U0[0].ba) ;
 umax=0.0 ;
 for (it=(*U0[0].ba).begin(); it<=(*U0[0].ba).end(); ++it) {

   (*U0[0].ba)[it]= 0.0 ;
   (*U0[1].ba)[it]= 0.0 ;
  
   if ((it.i[1]==0)||(it.i[1]==e[1])) U0[0][it]=U0[1][it]=0.0 ;

   umax = 1.0 ;
 }

 double w=0.0 ;
 SetUD[0]=new ConstFunction() ;SetUD[0]->SetParameters(&w) ;
 SetUD[1]=new ConstFunction() ;SetUD[1]->SetParameters(&w) ;
 
 SetF [0]=new ProhlFunctionF1() ;
 SetF [1]=new ProhlFunctionF2() ;
  
 SetBG() ;
}


//
// Driven Cavity flow
//
template<int DIM>
void Navier<DIM>::SetProblem4() {
 int i,j ;
 ReadWavelets() ;

 // Parameters
 strcpy(ID,"DrivC") ;
 nu   =1./1000. ;
 dt   =0.25e-3  ;
 T    =1000.0   ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=1000 ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=5 ;
   for (j=0; j<DIM; j++) BC[i][j][0]=BC[i][j][1]=0 ;

   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }
 SetBG() ;

 // Domain
 XL[0]=1.0  ;
 XL[1]=XL[0]  ;
 nu *=XL[0]*XL[0] ;

 // OutTimes
 for (i=0; i<=20 ;i++) OutTimes[i]=2*i ;
 OutTimes[i]=1e+100 ;
 prstep=20 ;

 // Movie
 MovieStep=MAXINT ;
 for (i=0; i<DIM; i++) {
   MovieJ[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Initial Conditions.

 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;
 P0.Init(a,e,&WC) ;

 for (j=0; j<DIM; j++) { 
   U0[j].Set(0.0) ; 
   U0[j].SetBoundaryConditions(BC[i]) ;
 }

 P0.Set(0.0) ;
 P0.SetBoundaryConditions(BG) ;

 // Values for Dirichlet-Boundary Conditions
 double w=0.0 ;
 SetUD[0]=new DrivCFunctionUD1() ;
 for (i=1; i<DIM; i++) { 
   SetUD[i]=new ConstFunction() ;
   SetUD[i]->SetParameters(&w)  ;
 }
 
 for (i=0; i<DIM; i++) { 
   SetF[i]=new ConstFunction() ;
   SetF[i]->SetParameters(&w) ;
 }
 
}



//
// Modons
//
template<int DIM>
void Navier<DIM>::SetProblem5() {
 int i,j,a[DIM],e[DIM] ;

 int    NV  =3 ;
 double px[]={3./8 , 5./8 , 5./8} ;
 double py[]={1./2 , 1./2 , (1+1/(2*sqrt(2.0)))/2 } ;
 //double am[]={-1*PI   , 2*PI   ,-1*PI} ;
 double am[]={PI   ,  PI   ,-0.5*PI} ;
 double si[]={1/(2*PI*PI) , 1/(2*PI*PI) ,1/(2*PI*PI)} ;
 double ommean ;

 // Parameters
 strcpy(ID,"Modons")  ;
 nu   =5e-5/(4*PI*PI) ;
 dt   =5.0e-3 ;
 T    =40.1   ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100 ;

 // OutTimes
 for (i=0; i<=20 ;i++) OutTimes[i]=2*i ;
 OutTimes[i]=1e+100 ;

 prstep=100 ;

 // Movie
 MovieStep=40 ;
 for (i=0; i<DIM; i++) {
   MovieJ[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0  ;
 XL[1]=XL[0]  ;
 for (i=0; i<NV; i++) {
   px[i] *=XL[0] ;
   py[i] *=XL[1] ;
   si[i] *=XL[0] ;
 }
 nu *=XL[0]*XL[0] ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=8 ;
   for (j=0; j<DIM; j++) BC[i][j][0]=-1 , BC[i][j][1]=WAVEPERIODICBC ;

   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 //printf("Allocated\n") ; fflush(stdout) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 Matrix<double,DIM>::iterator it(U0[0].ba) ;

 UniformData<DIM> Omega(a,e,&WC) , TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+2] , Psi(a,e,&WC) ;

 for (i=0;i<Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+2;i++) TMP[i].Init(a,e,&WC) ;

 ommean=0 ;
 for (it.Set(a); it<=e; ++it) {
   double x[DIM],r2,s2,dx,dy,om=0.0 ;
   for (i=0; i<DIM; i++) x[i]=(it.i[i]*XL[i])/e[i] ;
        
   for (i=0; i<NV; i++) {
     dx =x[0]-px[i] ;
     dy =x[1]-py[i] ;

     r2 =dx*dx + dy*dy ;
     s2 =si[i]*si[i];
     om+=am[i]*exp(-r2/s2) ;
   }

   for (i=0; i<DIM; i++) if (it.i[i]==e[i]) om=0 ;
   
   Omega[it]=om ;
   ommean  +=om ;
 }
 for (j=0; j<DIM; j++) Omega.ismultiscale[j]=false ;
 
 ommean /=(double)(1<<(2*J[0])) ;
 //printf("OmegaMean: %e\n",ommean) ;

 for (it.Set(a); it<=e; ++it) {
   bool f=true ;
   for (i=0; i<DIM; i++) if (it.i[i]==e[i]) f=false ;
   if (f) Omega[it] -=ommean ;
 }

 //ommean=0 ;
 //for (it.Set(a); it<=e; ++it) ommean +=Omega[it] ;
 //printf("Omega: %e\n",ommean) ;

 U0[0].SetBoundaryConditions(BG) ;
 U0[1].SetBoundaryConditions(BG) ;
 Omega.SetBoundaryConditions(BG) ;
   Psi.SetBoundaryConditions(BG) ;

 double f=0.0 ;
 for (i=0; i<DIM; i++) { 
   SetF [i]=new ConstFunction() ;
   SetUD[i]=new ConstFunction() ;
   SetF [i]->SetParameters(&f) ;
   SetUD[i]->SetParameters(&f) ;
 }
   
 // Berechnung von U,V uber Laplace Problem
 // Operatoren ...
 HelmholtzOperator< UniformData<DIM> , DIM> MO ;
 MO.par[0]=0.0 ;
 MO.Tmp =&TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2] ;
 MO.Tmp2=&TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+1] ;
 

 for (i=0; i<DIM; i++) MO.par[i+1]=1/(XL[i]*XL[i]) ;

 Solver<UniformData<DIM>,DIM> MS ;
 MS.eps          =1e-10 ;
 MS.maxiter      =150 ;
 MS.Op           =&MO ;
 MS.StopCriterion=StoppingCriterionAbsoluteResidual ;
 MS.print        ="MatBiCG:" ;
 MS.prstep       =1 ;

 for (i=0; i<Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2; i++) MS.Tmp[i]=&TMP[i] ;

 // Lap(Psi)=-omega ;
 for (j=0; j<DIM; j++) Omega.ApplyOp(BG[j] , &Omega , j , TRANSFORM) ;
 Omega.Mul(-1) ;

 MS.BiCGSTAB2(&Psi,&Omega) ;
 MS.BiCGSTAB2(&Psi,&Omega) ;

 //for (j=0; j<DIM; j++) Psi.ApplyOp(BG[j] , &Psi , j , INVERSETRANSFORM) ;
 Psi.ba->WriteMatLab("PSI") ;

 for (i=0; i<2; i++) {
   j= (i==0) ? 1 : 0 ;
   U0[i].ApplyOp(BG[j] , &Psi , j , FIRSTDERIVATIVE) ;

   if (i==0) U0[i].Mul( 1/XL[1]) ;
   if (i==1) U0[i].Mul(-1/XL[0]) ;
   for (j=0; j<DIM; j++) TMP[0].ApplyOp(BG[j] , (j==0) ? &U0[i] : &TMP[0] , j , INVERSETRANSFORM) ;
   printf("Umax: %e\n",TMP[0].MaxAbs()) ;

   if (i==0) TMP[0].ba->WriteMatLab("U0") ;
   if (i==1) TMP[0].ba->WriteMatLab("U1") ;

 }
 //exit(-1) ;
}



//
// mixing layer
//
template<int DIM>
void Navier<DIM>::SetProblem6() {
 int i,j,a[DIM],e[DIM] ;
 
 // Parameters
 strcpy(ID,"Mixing")  ;
 nu   =5e-5/(4*PI*PI) ;
 // dt   =2.5e-3 ;
 dt   =3*2.5e-3 ;
 T    =80.0   ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // IO
 prstep=100 ;

 for (i=0; i<=(int)(T/4+1e-10); i++) OutTimes[i]=4*i ;
 OutTimes[i++]=12.5   ;
 OutTimes[i++]=25     ;
 OutTimes[i++]=37.5   ;
 OutTimes[i  ]=1e+100 ;
 qsort(OutTimes , i , sizeof(double) , doublecmp) ;

 printf("OutTimes\n") ;
 for (j=0; j<=i ;j++) printf("%d %e\n",j,OutTimes[j]) ;

 // Movie
 MovieStep=50 ;
 MovieJ[0]=8 ;
 MovieJ[1]=9 ;
 for (i=0; i<DIM; i++) {
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0     ;
 XL[1]=2*XL[0] ;

 nu *=XL[0]*XL[0] ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   if (i==0) J[i]=8 ;
      else   J[i]=J[0]+1 ;

   for (j=0; j<DIM; j++) BC[i][j][0]=-1 , BC[i][j][1]=WAVEPERIODICBC ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }
 

 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 Matrix<double,DIM>::iterator it(U0[0].ba) ;

 UniformData<DIM> Omega(a,e,&WC) , TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+2] , Psi(a,e,&WC) ,
           *TMP1=&TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+1] ;

 for (i=0; i<Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2+2; i++) TMP[i].Init(a,e,&WC) ;  
 
 char name[1000] ;
 sprintf(name,"%s/../../tex/Berlin99/ICshear.256",_WAVELET_ROOT_) ;
 FILE *fid=fopen(name,"r") ;
 if (fid==NULL) {
   printf("Could not open: %s\n",name) ;
   exit(-1) ;
  }
 
 Omega.Set(0.0) ;
 double om ;
 for (it.Set(a); it<=e; ++it) {
   if ((it.i[0]<e[0])&&(it.i[1]<e[1]))
     if (it.i[1]<256) {
       int in[2]={it.i[1] , it.i[0] } ;
       fscanf(fid,"%le",&om) ;
       Omega[in]=om ;
     } else {
       int in[2]={it.i[0] , 511-it.i[1]} ;
       Omega[it]=-Omega[in] ;
     }
 }
 fclose(fid) ; 

 om=0 ;
 for (it.Set(a); it<=e; ++it) om +=Omega[it] ;
 printf("Omega: %e\n",om) ;
 double XA[2]={0},XE[2]={1,2} ;
 Omega.ba->WriteUDF("Omega",XA,XE,"omega",NULL,1,false,false) ;
 //exit(-1) ;

 U0[0].SetBoundaryConditions(BG) ;
 U0[1].SetBoundaryConditions(BG) ;
 Omega.SetBoundaryConditions(BG) ;
   Psi.SetBoundaryConditions(BG) ;

 double f=0.0 ;
 for (i=0; i<DIM; i++) { 
   SetF [i]=new ConstFunction() ;
   SetUD[i]=new ConstFunction() ;
   SetF [i]->SetParameters(&f) ;
   SetUD[i]->SetParameters(&f) ;
 }
   
 // Berechnung von U,V uber Laplace Problem
 // Operatoren ...
 HelmholtzOperator<UniformData<DIM>,DIM> MO ;
 MO.par[0]=0.0 ;
 MO.Tmp=&TMP[Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2] ;

 MO.par[1]=1/(XL[0]*XL[0]) ;
 MO.par[2]=1/(XL[1]*XL[1]) ;

 Solver<UniformData<DIM>,DIM> MS ;
 MS.eps          =1e-20 ;
 MS.maxiter      =150 ;
 MS.Op           =&MO ;
 MS.StopCriterion=StoppingCriterionAbsoluteResidual ;
 MS.print        ="MatBiCG:" ;
 MS.prstep       =1 ;

 for (i=0; i<Solver<UniformData<DIM>,DIM>::NUMTMP_BICGSTAB2; i++) MS.Tmp[i]=&TMP[i] ;

 // Lap(Psi)=-omega ;
 for (j=0; j<DIM; j++) TMP1->ApplyOp(BG[j] , (j==0) ? &Omega : TMP1 , j , TRANSFORM) ;
 TMP1->PPlus(-2.0 , TMP1)  ;

 MS.BiCGSTAB2(&Psi,TMP1) ;
 MS.BiCGSTAB2(&Psi,TMP1) ;

 for (j=0; j<DIM; j++) Psi.ApplyOp(BG[j] , &Psi , j , INVERSETRANSFORM) ;
 Psi.ba->WriteUDF("Psi",XA,XE,"psi",NULL,1,false,false) ;

 for (i=0; i<2; i++) {
   j= (i==0) ? 1 : 0 ;
   U0[i].ApplyOp(BG[j] , &Psi , j , FIRSTDERIVATIVE) ;

   if (i==0) U0[i].Mul( 1/XL[1]) ;
   if (i==1) U0[i].Mul(-1/XL[0]) ;
   printf("Umax: %e\n",U0[i].MaxAbs()) ;

   sprintf(name,"%s/Data/%s/U%s",_WAVELET_ROOT_, ID , (i==0) ? "0" : "1" ) ;
   U0[i].ba->WriteMatLab(name) ;
   U0[i].ba->WriteUDF(name,XA,XE,"u",NULL,1,false,false) ;
 }

 // Test der Rotation:
 Omega.ApplyOp(BG[0] , &U0[1] , 0 , FIRSTDERIVATIVE) ;
 Psi  .ApplyOp(BG[1] , &U0[0] , 1 , FIRSTDERIVATIVE) ;
 Omega.Add(1/XL[0],&Omega , -1/XL[1],&Psi) ;
 printf("OmegaMax: %e\n",Omega.MaxAbs()) ;
 Omega.ba->WriteUDF("Omega.dat",XA,XE,"omega",NULL,1,false,false) ;
}


//
// mixing layer mit Einlesen der Anfangsgeschwindigkeiten
//
template<int DIM>
void Navier<DIM>::SetProblem7() {
 int i,j,a[DIM],e[DIM];
 double XA[2]={0,0}   ;
 
 // Parameters
 strcpy(ID,"Mixing")  ;
 nu   =5e-5/(4*PI*PI) ;
 // dt   =2.5e-3 ;
 dt   =4*2.5e-3 ;
 T    =80.0     ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // IO
 prstep=100 ;

 for (i=0; i<=(int)(T/4+1e-10); i++) OutTimes[i]=40*i ;
 OutTimes[i++]=12.5   ;
 OutTimes[i++]=25     ;
 OutTimes[i++]=37.5   ;
 OutTimes[i  ]=1e+100 ;
 qsort(OutTimes , i , sizeof(double) , doublecmp) ;

 printf("OutTimes\n") ;
 for (j=0; j<=i ;j++) printf("%d %e\n",j,OutTimes[j]) ;

 // Movie
 MovieStep=200 ;
 MovieJ[0]=9   ;
 MovieJ[1]=10  ;
 for (i=0; i<DIM; i++) {
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0     ;
 XL[1]=2*XL[0] ;

 nu *=XL[0]*XL[0] ;

 // Boundary Conditions
 J[0]=9 ;
 J[1]=10 ;
 for (i=0;i<DIM;i++) {
   for (j=0; j<DIM; j++) BC[i][j][0]=-1 , BC[i][j][1]=WAVEPERIODICBC ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 double f=0.0 ;
 for (i=0; i<DIM; i++) { 
   SetF [i]=new ConstFunction() ;
   SetUD[i]=new ConstFunction() ;
   SetF [i]->SetParameters(&f) ;
   SetUD[i]->SetParameters(&f) ;
 }

 // aux. Matrix for Read
 int e1[2]={511 , 1023} ; 
 Matrix<double , DIM>         M(a,e1) ;
 Matrix<double,DIM>::iterator s(&M)   ;

 for (i=0; i<2; i++) {
   // Read Matrix
   char name[1000] ;
   //sprintf(name,"%s/../../C/Spectral/jobs/Mixing/U%s-1024",_WAVELET_ROOT_, (i==0) ? "1" : "0") ;
sprintf(name,"../../Data/Mixing/U%s-1024", (i==0) ? "1" : "0") ;
   M.ReadMatLab(name) ;
   
   // copy into larger initial matrix
   U0[i].Set(0.0) ;
   U0[i].SetBoundaryConditions(BC[i]) ;
   for (j=0; j<DIM; j++) U0[i].Ext.ismultiscale[j]=false ;

   for (s=M.begin(); s<= M.end(); ++s) (*U0[i].ba)[s.i]=M[s] ;
 
   
   // periodize 
   int in[2] , im[2] ,j ;
   in[0]=0 ;
   im[0]=e[0] ;
   for (j=0; j<=e[1]; j++) {
     in[1]=im[1]=j ;
     (*U0[i].ba)[im]=(*U0[i].ba)[in] ;
    }

   in[1]=0 ;
   im[1]=e[1] ;
   for (j=0; j<=e[0]; j++) {
     in[0]=im[0]=j ;
     (*U0[i].ba)[im]=(*U0[i].ba)[in] ;
    }

   in[0]=in[1]=0 ;
   im[0]=e[0] , im[1]=e[1] ;
   (*U0[i].ba)[im]=(*U0[i].ba)[in] ; 
   
   // Test-ausgaben:
   printf("Umax: %e\n",U0[i].MaxAbs()) ;

   sprintf(name,"%s/Data/%s/U%s",_WAVELET_ROOT_, ID , (i==0) ? "0" : "1" ) ;
   U0[i].ba->WriteUDF(name,XA,XL,"u",NULL,1,false,false) ;
 }

 // Trafo der Anfangsbedingungen
 for (i=0; i<DIM; i++)
   for (j=0; j<DIM; j++) U0[i].ApplyOp(BC[i][j] , &U0[i] , j , TRANSFORM) ; 

 // Test der Rotation
  
 UniformData<DIM> Tmp1(a,e,&WC) , Tmp2(a,e,&WC) ;
 Tmp1.SetBoundaryConditions(BG) ;
 Tmp2.SetBoundaryConditions(BG) ;

 Tmp1.ApplyOp(BG[0] , &U0[1] , 0 , FIRSTDERIVATIVE) ;
 Tmp2.ApplyOp(BG[1] , &U0[0] , 1 , FIRSTDERIVATIVE) ;
 Tmp1.Sub(&Tmp1 , &Tmp2) ;
 for (j=0; j<DIM; j++) Tmp1.ApplyOp(BG[0] , &Tmp1 , j , INVERSETRANSFORM) ;
 
 printf("Omega: %e  %e\n",Tmp1.MaxAbs(),Tmp1.InnerProd(&Tmp1)/(e[0]*e[1])/2)   ;
 Tmp1.ba->WriteUDF("Omega.dat",XA,XL,"omega",NULL,1,false,false) ;

}



//
// Merging Modons layer mit Einlesen der Anfangsgeschwindigkeiten
//
template<int DIM>
void Navier<DIM>::SetProblem8() {
 int i,j,a[DIM],e[DIM] ;
 
 // Parameters
 strcpy(ID,"Merging")  ;
 nu   =5e-5/(4*PI*PI) ;
 // dt   =2.5e-3 ;
 dt   =1.25e-3/2 ;
 T    =40.1      ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // IO
 prstep=1000 ;

 for (i=0; i<=(int)(T/5+1e-10); i++) OutTimes[i]=5*i ;
 OutTimes[i  ]=1e+100 ;

 printf("OutTimes\n") ;
 for (j=0; j<=i ;j++) printf("%d %e\n",j,OutTimes[j]) ;

 // Movie
 MovieStep=MAXINT ;
 MovieJ[0]=9     ;
 MovieJ[1]=10    ;
 for (i=0; i<DIM; i++) {
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0 ;
 XL[1]=1.0 ;

 nu *=XL[0]*XL[0] ;

 // Boundary Conditions
 J[0]=9 ;
 J[1]=9 ;
 for (i=0;i<DIM;i++) {
   for (j=0; j<DIM; j++) BC[i][j][0]=-1 , BC[i][j][1]=WAVEPERIODICBC ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 double f=0.0 ;
 for (i=0; i<DIM; i++) { 
   SetF [i]=new ConstFunction() ;
   SetUD[i]=new ConstFunction() ;
   SetF [i]->SetParameters(&f) ;
   SetUD[i]->SetParameters(&f) ;
 }

 // aux. Matrix for Read
 int px=2048 , rx=px/(1<<J[0]) ;
 int e1[2]={px-1 , px-1} ;
 Matrix<double , DIM>         M(a,e1) ;
 Matrix<double,DIM>::iterator s(&M)   ;

 for (i=0; i<2; i++) {
   // Read Matrix
   char name[1000] ;
   sprintf(name,"%s/../../C/Spectral/jobs/Merging/U%s-%d-Intel",_WAVELET_ROOT_, (i==0) ? "0" : "1",px) ;
   M.ReadMatLab(name) ;
   
   // copy into larger initial matrix
   U0[i].Set(0.0) ;
   U0[i].SetBoundaryConditions(BC[i]) ;

   for (s=M.begin(); s<= M.end(); ++s) {
     bool cp=true ;
     int is[DIM]  ;
     for (j=0;j<DIM; j++) {
       is[1-j]=s.i[j]/rx ;
       if ( (s.i[j]%rx)!=0 )  cp=false ;
      }
     if (cp) (*U0[i].ba)[is]=M[s] ;
   }
   
   // periodize 
   int in[2] , im[2] ,j ;
   in[0]=0 ;
   im[0]=e[0] ;
   for (j=0; j<=e[1]; j++) {
     in[1]=im[1]=j ;
     (*U0[i].ba)[im]=(*U0[i].ba)[in] ;
    }

   in[1]=0 ;
   im[1]=e[1] ;
   for (j=0; j<=e[0]; j++) {
     in[0]=im[0]=j ;
     (*U0[i].ba)[im]=(*U0[i].ba)[in] ;
    }

   in[0]=in[1]=0 ;
   im[0]=e[0] , im[1]=e[1] ;
   (*U0[i].ba)[im]=(*U0[i].ba)[in] ; 
   
   // Test-ausgaben:
   printf("Umax: %e\n",U0[i].MaxAbs()) ;

   sprintf(name,"%s/Data/%s/U%s",_WAVELET_ROOT_, ID , (i==0) ? "0" : "1" ) ;
   U0[i].ba->WriteMatLab(name) ;
 }

 // Test der Rotation
  
 UniformData<DIM> Tmp1(a,e,&WC) , Tmp2(a,e,&WC) ;
 Tmp1.SetBoundaryConditions(BG) ;
 Tmp2.SetBoundaryConditions(BG) ;

 Tmp1.ApplyOp(BG[0] , &U0[1] , 0 , FIRSTDERIVATIVE) ;
 Tmp2.ApplyOp(BG[1] , &U0[0] , 1 , FIRSTDERIVATIVE) ;
 Tmp1.Add(1/XL[0],&Tmp1 , -1/XL[1],&Tmp2) ;
 printf("Omega: %e  %e\n",Tmp1.MaxAbs(),Tmp1.InnerProd(&Tmp1)/(e[0]*e[1])/2)   ;
 Tmp1.ba->WriteMatLab("Omega.dat") ;

}


//
// Stokes Test Problem
//
template<int DIM>
void Navier<DIM>::SetProblemStokes() {
 int i,j,a[DIM],e[DIM] ;

 // Parameters
 strcpy(ID,"Stokes")  ;
 nu   =1e-1   ;
 dt   =1.0e-1 ;
 T    =400.1  ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=8 ;
   for (j=0; j<DIM; j++) BC[i][j][0]=BC[i][j][1]=0 ;

   a[i]=0 ;
   e[i]=(1<<J[i]) ;
  }

 // OutTimes
 for (i=0; i<=20 ;i++) OutTimes[i]=2*i ;
 OutTimes[i]=1e+100 ;
 prstep=20 ;

 // Movie
 MovieStep=MAXINT ;
 for (i=0; i<DIM; i++) {
   MovieJ[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0  ;
 XL[1]=XL[0]  ;
 nu *=XL[0]*XL[0] ;
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 SetF[0]=new StokesFunctionF1 ;
 SetF[1]=new StokesFunctionF2 ;

 double f=0.0 ;
 for (i=0; i<DIM; i++) { 
   SetUD[i]=new ConstFunction() ;
   SetUD[i]->SetParameters (&f) ;
 }

 for (i=0; i<DIM; i++) {
   SetF[i]->SetParameters(&nu)        ;
   U0[i].Set(0.0)                     ;
   U0[i].SetBoundaryConditions(BC[i]) ;
 }
}




template<int DIM>
void Navier<DIM>::SetProblemStokes2() {
 int i,j,a[DIM],e[DIM] ;

 // Parameters
 strcpy(ID,"Stokes")  ;
 nu   =0.000100   ;
 dt   =1.0e-3 ;
 T    =1.0    ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
   BC[i][0][0]=-1 ;
   BC[i][0][1]=WAVEPERIODICBC ;
   BC[i][1][0]=-1 ;
   BC[i][1][1]=-1 ;
  }
 


 // OutTimes
 for (i=0; i<=20 ;i++) OutTimes[i]=2*i ;
 OutTimes[i]=1e+100 ;
 prstep=20 ;

 // Movie
 MovieStep=MAXINT ;
 for (i=0; i<DIM; i++) {
   MovieJ[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0  ;
 XL[1]=XL[0]  ;
 nu *=XL[0]*XL[0] ;
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 double pp[2]={nu,0.0} ;
 SetF[0]=new StokesFunction2F1 ;
 SetF[1]=new StokesFunction2F2 ;
 SetUD[0]=new StokesFunction2U1 ;
 SetUD[1]=new StokesFunction2U2 ;

 SetF [0]->SetParameters(pp) ;
 SetF [1]->SetParameters(pp) ;
 SetUD[0]->SetParameters(pp) ;
 SetUD[1]->SetParameters(pp) ;

 for (i=0; i<DIM; i++) {
   U0[i].SetFunction(SetUD[i]) ;
   U0[i].SetBoundaryConditions(BC[i]) ;
   for (j=0; j<DIM; j++) U0[i].ApplyOp(U0[i].BC[j], &U0[i], j, TRANSFORM) ;
 }
}


template<int DIM>
void Navier<DIM>::SetProblemStokes3() {
 int i,j,a[DIM],e[DIM] ;

 // Parameters
 strcpy(ID,"Stokes")  ;
 nu   =1.0   ;
 dt   =1.0e-3 ;
 T    =1.0    ;

 gamma=1e+100 ;

 maxres =1e-5 ;
 maxiter=100  ;

 // Boundary Conditions
 for (i=0;i<DIM;i++) {
   J[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<J[i]) ;
   BC[i][0][0]=-1 ;
   BC[i][0][1]=-1 ;
   BC[i][1][0]=-1 ;
   BC[i][1][1]=-1 ;
  }
 


 // OutTimes
 for (i=0; i<=20 ;i++) OutTimes[i]=2*i ;
 OutTimes[i]=1e+100 ;
 prstep=20 ;

 // Movie
 MovieStep=MAXINT ;
 for (i=0; i<DIM; i++) {
   MovieJ[i]=8 ;
   a[i]=0 ;
   e[i]=(1<<MovieJ[i]) ;
  }
 MovieTmp.Init(a,e,&WC) ;

 // Domain
 XL[0]=1.0  ;
 XL[1]=XL[0]  ;
 nu *=XL[0]*XL[0] ;
 
 this->ReadWavelets() ;
 this->SetBG() ;
 this->Print() ;

 // Initial Conditions.
 for (i=0; i<DIM; i++) U0[i].Init(a,e,&WC) ;

 P0.Init(a,e,&WC) ;
 P0.SetBoundaryConditions(BG) ;
 P0.Set(0.0) ;

 double pp[2]={nu,0.0} ;
 SetF[0]=new StokesFunction3F1 ;
 SetF[1]=new StokesFunction3F2 ;
 SetUD[0]=new StokesFunction3U1 ;
 SetUD[1]=new StokesFunction3U2 ;

 SetF [0]->SetParameters(pp) ;
 SetF [1]->SetParameters(pp) ;
 SetUD[0]->SetParameters(pp) ;
 SetUD[1]->SetParameters(pp) ;

 for (i=0; i<DIM; i++) {
   U0[i].SetFunction(SetUD[i]) ;
   U0[i].SetBoundaryConditions(BC[i]) ;
   for (j=0; j<DIM; j++) U0[i].ApplyOp(U0[i].BC[j], &U0[i], j, TRANSFORM) ;
 }
}

---