sandbox/Antoonvh/adaptLES.h

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    
    //Adapt initial grids
    
    void adaptLESinit(int MAXLEVEL)
    {
      unrefine(y>L0/4 && level>=(MAXLEVEL-1),all);
      unrefine(y>200 && level>=(MAXLEVEL-2),all); 
    }
    
    //Adapt grid according to solution while running: 
    //SGS-flux is compared to resolved-flux to check if the solution is really a LES. If not -> REFINE
    //If the SGS-flux is neglegible to resolved flux we are not taking full benefit from LES -> Coarsen
    //When fluxes vanish (laminarization) check if Ed.vis < Mol.vis 
    //Alse check if there is any resolved flux otherwise i am not sure what to do  
    
    
    void adaptLESrun(int MAXLEVEL, vector u,scalar T, scalar Evis, int lev[], double yc[]) 
    {
      fprintf(ferr,"hoi\n");
      double cur = 0.1, cr=0.5;
      scalar fmr[] , fms[],uh[];	
      foreach()
        {	
          uh[]=sqrt(sq(u.x[])+sq(u.z[]));
          fmr[] = sqrt(sq(u.y[]*(uh[0,1,0]-uh[0,-1,0])/Delta));
          fms[] = sqrt(sq(Evis[]*(uh[0,1,0]-(2*uh[])+uh[0,-1,0])/(sq(Delta))));  	
        }
    	
      fprintf(ferr,"hoi\n");
      int k = 0;
      double xp,zp,Url[1<<MAXLEVEL]; 
      int uri = 0;
      while (yc[k] != 0)	
        {
    #ifdef _OPENMP
    #pragma omp parallel for reduction(+:F) default(shared)
    #endif
          double afmr=0 , afms=0 ;
          double yp = yc[k];
          for (int i = 0; i < pow(2,lev[k]); i++)
    	{	
    	  xp=X0+L0/((pow(2,lev[k]+1))+i*(L0/(pow(2,lev[k]))));	
    	  for (int j = 0; j < pow(2,lev[k]); j++)
    	    {
    	      zp=Z0+L0/((pow(2,lev[k]+1))+i*(L0/(pow(2,lev[k]))));						Point point = locate (xp, yp,zp);
    	      afmr += fmr[];
    	      afms += fms[];
    	    }
    	}
          if (cur * afmr > afms )
    	{
    	  Url[uri] = yc[k];
    	  if (uri>0)
    	    {	
    	      if ( (Url[uri]-Url[uri-1]) == (L0/(pow(2,lev[k]))))
    		unrefine(y == ( (Url[uri]+Url[uri]) /2 ),all);
    	      
    	    }
    	  uri++;
    	}
          else if (cr * afmr < afms && (afmr/pow(2,2*lev[k])) > 0.01)
    	{
    	  refine(y>(lev[k]-0.1) && y<(lev[k]+0.1),all);
    	}
          k++;
          boundary(all); 
        }
    
    }