#include "udf.h"
#include "mem.h"
#include "metric.h"       
DEFINE_TURBULENT_VISCOSITY(sgsvre,cell,thread)
{
/*     By A.W. Vreman, Vreman Research 2005. 
       Implementation of the isotropic form of the LES-model published
       in Physics of Fluids 16, 3670-3681 (2004). 
       In the form of a user-defined function for Fluent.
*/

       real csma,coef,del,d1,d2,d3;
       real d1v1,d2v1,d3v1,d1v2,d2v2,d3v2,d1v3,d2v3,d3v3;
       real b11,b12,b13,b22,b23,b33;
       real abeta,bbeta,eddyvis;

/*      csma: Smagorinsky reference constant, e.g. 0.20, 0.17 or 0.10. 
*/
       csma=0.17e0;
       coef=2.5e0*csma*csma;
 

        del=C_DIAMETER(cell,thread);
        d1=del*del;
        d2=d1;
        d3=d1; 
/*
For anisotropic form prescribe d1, d2 and d3 by h1**2, h2**2 and h3**2, 
where h1, h2 and h3 represent the meshsizes in the three coordinate directions
in physical space.
*/
        d1v1=C_U_G(cell,thread)[0];
        d2v1=C_U_G(cell,thread)[1];
        d3v1=C_U_G(cell,thread)[2];

        d1v2=C_V_G(cell,thread)[0];
        d2v2=C_V_G(cell,thread)[1];
        d3v2=C_V_G(cell,thread)[2];

        d1v3=C_W_G(cell,thread)[0];
        d2v3=C_W_G(cell,thread)[1];
        d3v3=C_W_G(cell,thread)[2];

        b11=d1*d1v1*d1v1+d2*d2v1*d2v1+d3*d3v1*d3v1;
        b12=d1*d1v1*d1v2+d2*d2v1*d2v2+d3*d3v1*d3v2;        
        b13=d1*d1v1*d1v3+d2*d2v1*d2v3+d3*d3v1*d3v3;
        b22=d1*d1v2*d1v2+d2*d2v2*d2v2+d3*d3v2*d3v2;
        b23=d1*d1v2*d1v3+d2*d2v2*d2v3+d3*d3v2*d3v3;
        b33=d1*d1v3*d1v3+d2*d2v3*d2v3+d3*d3v3*d3v3;

        abeta=d1v1*d1v1+d1v2*d1v2+d1v3*d1v3+
              d2v1*d2v1+d2v2*d2v2+d2v3*d2v3+
              d3v1*d3v1+d3v2*d3v2+d3v3*d3v3;        
        bbeta=b11*b22-b12*b12+b11*b33-b13*b13+b22*b33-b23*b23;

        if (bbeta<1e-12) 
         {eddyvis=0.;}
        else
         eddyvis=coef*sqrt(bbeta/abeta);
/*  in case of single precision (real*4): "if (bbeta.lt.1e-6) " 
*/

return eddyvis;
}