/************************************************************************/
/* UDF for specifying fully-developed profile boundary conditions       */
/************************************************************************/

#include "udf.h"

#define YMIN 0.0
#define YMAX 0.2
#define UMEAN 3.7
#define B 1./7.
#define DELOVRH 0.5
#define VISC 1.7894e-05
#define CMU 0.09
#define VKC 0.41


/*  profile for x-velocity    */


DEFINE_PROFILE(x_velocity, thread, position)
{
  real y, del, h, x[ND_ND], ufree;
  face_t f;

  h = YMAX - YMIN;
  del = DELOVRH*h;
  ufree = UMEAN*(B+1.);

  begin_f_loop(f, thread)
    {
      F_CENTROID(x,f,thread);
      y = x[1];

      if (y <= del)
         F_PROFILE(f,thread,position) = ufree*pow(y/del,B);
      else
         F_PROFILE(f,thread,position) = ufree*pow((h-y)/del,B);
    }
  end_f_loop(f, thread)
}


/*  profile for kinetic energy  */


DEFINE_PROFILE(k_profile, thread, position)
{
  real y, del, h, ufree, x[ND_ND];
  real ff, utau, knw, kinf;
  face_t f;

  h = YMAX - YMIN;
  del = DELOVRH*h;
  ufree = UMEAN*(B+1.);
  ff = 0.045/pow(ufree*del/VISC,0.25);
  utau=sqrt(ff*pow(ufree,2.)/2.0);
  knw=pow(utau,2.)/sqrt(CMU);
  kinf=0.002*pow(ufree,2.);


  begin_f_loop(f, thread)
    {
      F_CENTROID(x,f,thread);
      y=x[1];

      if (y <= del)
        F_PROFILE(f,thread,position)=knw+y/del*(kinf-knw);
      else
        F_PROFILE(f,thread,position)=knw+(h-y)/del*(kinf-knw);
    }
  end_f_loop(f, thread)
}


/* profile for dissipation rate  */


DEFINE_PROFILE(dissip_profile, thread, position)
{
  real y, x[ND_ND], del, h, ufree;
  real ff, utau, knw, kinf;
  real mix, kay;
  face_t f;

  h = YMAX - YMIN;
  del = DELOVRH*h;
  ufree = UMEAN*(B+1.);
  ff = 0.045/pow(ufree*del/VISC,0.25);
  utau=sqrt(ff*pow(ufree,2.)/2.0);
  knw=pow(utau,2.)/sqrt(CMU);
  kinf=0.002*pow(ufree,2.);

  begin_f_loop(f, thread)
    {
      F_CENTROID(x,f,thread);
      y=x[1];

      if (y <= del)
        kay=knw+y/del*(kinf-knw);
      else 
        kay=knw+(h-y)/del*(kinf-knw);

      if (VKC*y < 0.085*del)
        mix = VKC*y;
      else
        mix = 0.085*del;

      F_PROFILE(f,thread,position)=pow(CMU,0.75)*pow(kay,1.5)/mix;
    }
  end_f_loop(f,thread)
}