# resolution Poiseuille Turbulent 2D
#
#  PYL 12/10
#
# http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC5/
#
#
Define Rstar 500
Define ulog(y)  log(0.82*y + pow(1 + 4*pow(0.41,2)*pow(y,2),0.5))*pow(0.41,-1) + (pow(0.41,-2)*pow(y,-1)*(1 - pow(1 + 4*pow(0.41,2)*pow(y,2),0.5)))/2.

#
1 1 GfsSimulation GfsBox GfsGEdge {
	 y=0.5 } { 
	Time {end = 60 }
		
# attention la boite est L0xL0  
    PhysicalParams { L = 1. }
	
 
      Refine {
      if (fabs(y) <= 0.05 )
	  return 9;
      if (fabs(y) > 0.05 &&  fabs(y) <= 0.1)
	  return 8;
	  if (fabs(y) > 0.10 &&  fabs(y) <= 0.2)
	  return 7;
      if (fabs(y) > 0.2 &&  fabs(y) < 0.4)
      return 6;
      return 5;
  }

      Init { istep = 1 } {
      DY = dy("U");	
      moinsupvp=	( 0.41*0.41*(y*y) *sqrt(2.)*D2)*sqrt(2.)*D2
      }   
      
# use backward Euler to avoid Crank-Nicholson oscillations in time
     SourceViscosity {} {
// # longueur de Prandtl   	 
    	double Mu = 1; 
    	 Mu = 1./Rstar + 0.41*0.41*sqrt(2.)*(y*y)*D2;
    	return Mu;
    } {
    	beta = 1 
    	tolerance = 1e-4
    }
    

    Init {} {U =  ulog(y*Rstar) }
     
    Source U 1
     
     
# we need this so that acceleration can be balanced by viscous stress
    AdvectionParams { gc = 1 }
    AdvectionParams { scheme = none }    
#   AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }
    EventStop { istep = 1 } U 1.5e-6 DU
    ProjectionParams { tolerance = 1e-6 }
    OutputSimulation { istep = 5 } stdout
    OutputTime { step= 0.25  } stderr
    OutputLocation { step = 0.25 } vals.data cuty.dat
    OutputSimulation { step = 0.25 } SIM/sim-%g.txt { format = text } 
    EventScript { step = 0.25 } {
	 cp SIM/sim-$GfsTime.txt sim.data }
}
GfsBox {
    bottom = Boundary {
        BcDirichlet U 0
        BcDirichlet V 0 
    }
    top = Boundary {
        BcNeumann U 0
        BcDirichlet V 0 
    }
}
1 1 right

#  p"< awk '{if(($1<0.01)&&($1>0)){print $0}}' sim.data" u 2:6 w l
# set logscale x;p[0.005:]"< awk '{if($1>30){print $0}}' vals.data" u 3:7 w p