Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

2DMergingInitial.cc

//
// generate initial conditions for 2D merging vortices simulation
// the initial velocities are computed from an initial vorticity 
// distribution
//

#include "UniformData.hpp"
#include "Operators.hpp"
#include "Solver.hpp"

#define D 2 
int debugRefine ;

int main() {
  printf("2DMergingInitial\n") ;
  // which wavelets
  Wavelets WC("Interpolet",6) ;
  // box
  double XA[2]={0,0},XE[2]={1,1} ;
  // refinement level for initial values
  int L[D]={8,8} ;
  // boundary conditions
  int i,j,BC[D][2] ;
  for (i=0; i<D; i++) BC[i][0]=-1, BC[i][1]=PERIODIC ;

  // initial distribution of vortices: 3 Gaussian blobs
  int    NV  =3 ; // number of vortices
  double px[]={3./8 , 5./8 , 5./8                  },    // x coordinates
         py[]={1./2 , 1./2 , (1+1/(2*sqrt(2.0)))/2 },    // y coordinates
         am[]={PI   ,  PI   ,-0.5*PI},                   // amplitude
         si[]={1/(2*PI*PI) , 1/(2*PI*PI) ,1/(2*PI*PI)},  // variance
         ommean=0;

  // allocate uniform data structures
  int a[D]={0},e[D] ;
  for (i=0; i<D; i++) e[i]=1<<L[i] ;
  UniformData<D> U[D],Omega(a,e,&WC,XA,XE),Psi(a,e,&WC,XA,XE), TMP[Solver<UniformData<D>,D>::NUMTMP_BICGSTAB2+2],P(a,e,&WC,XA,XE) ;
  for (i=0; i<D; i++) U[i].Init(a,e,&WC,XA,XE) ;
  for (i=0; i<Solver<UniformData<D>,D>::NUMTMP_BICGSTAB2+2; i++) TMP[i].Init(a,e,&WC,XA,XE) ;

  // set initial vorticity distribution
  Matrix<double,D>::iterator s(a,e) ;
  for (s.Set(a); s<=e; ++s) {
    double x[D],r2,s2,dx,dy,om=0.0 ;
    for (i=0; i<D; i++) x[i]=(s.i[i]*XE[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<D; i++) if (s.i[i]==e[i]) om=0 ;
    
    Omega[s]=om ;
    ommean +=om ;
  }
  // set flags/bcs
  for (i=0; i<D; i++) Omega.Ext.ismultiscale[i]=false ,  Omega.Ext.islifting[i]=false ;
  Omega.SetBoundaryConditions(BC) ;

  // enforce vanishing mean vorticity
  ommean /=(1<<(2*L[0])) ;
  for (s.Set(a); s<=e; ++s) {
    bool f=true ;
    for (i=0; i<D; i++) if (s.i[i]==e[i]) f=false ;
    if (f) Omega[s] -=ommean ;
  }

  // compute stream function psi by  \Delta(psi) - Omega
  HelmholtzOperator< UniformData<D> , D> MO ;
  MO.Tmp =&TMP[Solver<UniformData<D>,D>::NUMTMP_BICGSTAB2] ;
  MO.Tmp2=&TMP[Solver<UniformData<D>,D>::NUMTMP_BICGSTAB2+1] ;
  
  MO.par[0]=0.0 ;
  for (i=0; i<D; i++) MO.par[i+1]=1/(XE[i]*XE[i]) ;
  
  Solver<UniformData<D>,D> MS ;
  MS.eps          =1e-10 ;
  MS.maxiter      =150 ;
  MS.Op           =&MO ;
  MS.StopCriterion=StoppingCriterionAbsoluteResidual ;
  MS.print        ="BiCG:" ;
  MS.prstep       =10 ;
  for (i=0; i<Solver<UniformData<D>,D>::NUMTMP_BICGSTAB2; i++) MS.Tmp[i]=&TMP[i] ;

  for (i=0; i<D; i++) Omega.ApplyOp(BC[i] , &Omega , i , TRANSFORM) ;
  Omega.Mul(-1) ;

  // initial guess for Psi
  ConstFunction F ;
  ommean=0.0 ;
  F.SetParameters(&ommean) ;
  Psi.SetBoundaryConditions(BC) ;
  Psi.SetFunction(&F) ;

  MS.BiCGSTAB2(&Psi,&Omega) ;
  MS.BiCGSTAB2(&Psi,&Omega) ;
  
  // compute velocity components
  for (i=0; i<2; i++) {
    j= (i==0) ? 1 : 0 ;
    U[i].ApplyOp(BC[j] , &Psi , j , FIRSTDERIVATIVE) ;
    
    if (i==0) U[i].Mul( 1/XE[1]) ;
    else      U[i].Mul(-1/XE[0]) ;
  }
  
  // set pressure
  P.SetBoundaryConditions(BC) ;
  P.SetFunction(&F) ;

   // write 
  UniformData<D> *MM[D+1]={&U[0],&U[1],&P} ;
  U[0].WriteSparse("../../Data/Merging2D/UVWP0", MM, D+1) ;
  
}

---