[tstorm]

I'm working with several UDF turbulence models, all variations on the v2-f model. Would anyone mind helping me parallelize the model? The Fluent UDF manual is confusing to me.

Thanks in advance.

Code starts here:


#include "udf.h"
#include "math.h"

/* Turbulence model constants */
const real C_MU=0.22;
const real CE_D=0.050;
const real CE_2=1.92;
const real SIG_K=1.0;
const real SIG_E=1.3;
const real C_1=1.4;
const real C_2=0.3;
const real C_L=0.23;
const real C_ETA=70.0;
const real y_star_limit=1.e32;
const real C_MU_ke=0.09;
const real zeta=8.0/3.0;
const real I_inlet=0.1; /* turbulent intensity Inlet*/
const real I_backflow=0.1; /* turbulent intensity outlet*/
const real mu_t_mu_l_backflow=10.; /* turbulent viscosity ratio mu_t / mu_l Backflow*/
const real mu_t_mu_l_inlet=10.; /* turbulent viscosity ratio mu_t / mu_l Inlet*/



/* User-defined scalars */
enum
{
K,
E,
V2,
F,
R_11,
R_22,
R_33,
R_12,
N_REQUIRED_UDS
};


DEFINE_SOURCE(u_source, c, t, dS, eqn)
{
real source;

dS[eqn]= 0.;
source = - C_R(c,t) * ( C_UDSI_G(c,t,R_11)[0] + C_UDSI_G(c,t,R_12)[1] );
return source;
}



DEFINE_SOURCE(v_source, c, t, dS, eqn)
{
real source;

dS[eqn]= 0.;
source = - C_R(c,t) * ( C_UDSI_G(c,t,R_12)[0] + C_UDSI_G(c,t,R_22)[1] );
return source;
}




DEFINE_SOURCE(k_source, c, t, dS, eqn)
{
real source;


dS[eqn]=- C_R(c,t)*C_R(c,t)*C_MU*C_UDSI(c,t,V2)/C_UDMI(c,t,0);
source=C_UDMI(c,t,3) - C_R(c,t)*C_R(c,t)*C_MU*C_UDSI(c,t,V2)*C_UDSI(c,t,K )/C_UDMI(c,t,0);
return source;
}



DEFINE_SOURCE(e_source, c, t, dS, eqn)
{
real CE_1, d, L;
real source;

C_UDSI(c,t,K)=MAX(C_UDSI(c,t,K),1.0e-12);

d=C_WALL_DIST(c,t);
L=C_UDMI(c,t,2);

CE_1=1.4 * (1.0+CE_D*sqrt( C_UDSI(c,t,K) / C_UDSI(c,t,V2) ));


dS[eqn]= -CE_2*C_R(c,t)/C_UDMI(c,t,1);
source= CE_1*C_UDMI(c,t,3)/C_UDMI(c,t,1) - CE_2*C_R(c,t)*C_UDSI(c,t,E)/C_UDMI(c,t,1);
return source;
}






DEFINE_SOURCE(v2_source, c, t, dS, eqn)
{
real source;
real kf;

C_UDSI(c,t,E)=MAX(C_UDSI(c,t,E),1.0e-12);


kf=C_UDSI(c,t,K)*C_UDSI(c,t,F);

dS[eqn]=-6.*C_R(c,t)*C_UDSI(c,t,E)/C_UDSI(c,t,K);
source=C_R(c,t)*kf-6.*C_R(c,t)*C_UDSI(c,t,V2)*C_UDSI(c,t,E)/C_UDSI(c,t,K);
return source;
}



DEFINE_SOURCE(f_source, c, t, dS, eqn)
{
real source;
C_UDSI(c,t,V2)=MAX(C_UDSI(c,t,V2),1.0e-12);

dS[eqn]=-1.0/SQR(C_UDMI(c,t,2));
source=-C_UDSI(c,t,F)/SQR(C_UDMI(c,t,2))-1.0/(SQR(C_UDMI(c,t,2))*C_UDMI(c,t,1))*((C_1-6.)*C_UDSI(c,t,V2)/C_UDSI(c,t,K)-2./3.*(C_1-1.))+C_2*C_UDMI(c,t,3)/(C_R(c,t)*C_UDSI(c,t,K)*SQR(C_UDMI(c,t,2)));
return source;
}




DEFINE_DIFFUSIVITY(ke_v2f_diffusivity, c, t, eqn)
{
real diff;

switch (eqn)
{
case K:
diff=C_UDMI(c,t,0)/SIG_K+C_MU_L(c,t);
break;

case E:
diff=C_UDMI(c,t,0)/SIG_E+C_MU_L(c,t);
break;

case V2:
diff=C_UDMI(c,t,0)+C_MU_L(c,t);
break;

case F:
diff=1.0;
break;

default:
diff=0.0;
}
return diff;
}

DEFINE_UDS_FLUX(user_flux, f, t, eqn)
{

switch (eqn)
{
case K:
return F_FLUX(f,t);
break;

case E:
return F_FLUX(f,t);
break;

case V2:
return F_FLUX(f,t);
break;

case F:
return 0.0;
break;

default:
return 0.0;
}

}



DEFINE_ADJUST(ke_adjust,domain)
{

Thread *t;
cell_t c;
real y_star, L, T, V, T1;
real mu_t;
real S11, S12, S13, S21, S22, S23, S31, S32, S33;
real W11, W12, W13, W21, W22, W23, W31, W32, W33;
real S, W, eta1, eta2, eta3;
real P11, P22, P33, Pkk;
real beta1, gamma1, cmu1, cmu2, cmu3;

/* Set the turbulent viscosity */
thread_loop_c(t,domain)
if (FLUID_THREAD_P(t))
{

if (NULL != THREAD_STORAGE(t,SV_UDSI_G(0)))
{


begin_c_loop(c,t)
{


y_star=C_R(c,t)*pow(C_MU_ke,0.25)*sqrt(C_UDSI(c,t, K))*C_WALL_DIST(c,t)/C_MU_L(c,t);

if (y_star <= y_star_limit)
{
T=MAX( C_UDSI(c,t,K)/C_UDSI(c,t,E), 6.0*sqrt( C_MU_L(c,t)/(C_R(c,t)*C_UDSI(c,t,E))) );
L=C_L*MAX(pow( (C_UDSI(c,t,K)),(3./2.))/C_UDSI(c,t,E),C_ETA*pow( ((pow(C_MU_L(c,t)/C_R(c,t),3.)/C_UDSI(c,t,E))),0.25));
}
else
{
T=C_UDSI(c,t,K)/C_UDSI(c,t,E);
L=C_L*pow( C_UDSI(c,t,K),(3./2.))/C_UDSI(c,t,E);
}



V=MAX( 2./3.-C_UDSI(c,t,V2)/C_UDSI(c,t,K) , 0.);

S11 = 0.5*( C_DUDX(c,t)+C_DUDX(c,t) );
S12 = 0.5*( C_DUDY(c,t)+C_DVDX(c,t) );
S13 = 0.5*( 0.+0. );
S21 = 0.5*( C_DVDX(c,t)+C_DUDY(c,t) );
S22 = 0.5*( C_DVDY(c,t)+C_DVDY(c,t) );
S23 = 0.5*( 0.+0. );
S31 = 0.5*( 0.+0. );
S32 = 0.5*( 0.+0. );
S33 = 0.5*( 0.+0. );
W11 = 0.5*( C_DUDX(c,t)-C_DUDX(c,t) );
W12 = 0.5*( C_DUDY(c,t)-C_DVDX(c,t) );
W13 = 0.5*( 0.-0. );
W21 = 0.5*( C_DVDX(c,t)-C_DUDY(c,t) );
W22 = 0.5*( C_DVDY(c,t)-C_DVDY(c,t) );
W23 = 0.5*( 0.-0. );
W31 = 0.5*( 0.-0. );
W32 = 0.5*( 0.-0. );
W33 = 0.5*( 0.-0. );

S=S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33;
W=W11*W11+W12*W12+W13*W13+W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33;

eta1=SQR(T)*S;
eta2=SQR(T)*W;
eta3=eta1-eta2;

beta1=1./(0.1+sqrt(eta1*eta2));
gamma1=1./(0.1+eta1);

cmu1=C_MU;
cmu2=6.0/5.0*sqrt(MAX(1.0-SQR(cmu1*C_UDSI(c,t,V2)/C_UDSI(c,t,K))*2.*eta1,0.))/(beta1+sqrt(eta1*eta2));
cmu3=6.0/5.0/(gamma1+eta1);

T1=MIN( T , 2./3.*C_UDSI(c,t,K)/(cmu1*C_UDSI(c,t,V2)*sqrt(S*zeta)) );
mu_t= C_R(c,t)*C_MU*C_UDSI(c,t,V2)*T;


C_UDSI(c,t,R_11)= -V*C_UDSI(c,t,K)*SQR(T)*
( cmu2*(S11*W11+S12*W21+S13*W31+S11*W11+S12*W21+S13* W31) - cmu3*(S11*S11+S12*S12+S13*S13-1./3.*S) );

C_UDSI(c,t,R_22)= -V*C_UDSI(c,t,K)*SQR(T)*
( cmu2*(S21*W12+S22*W22+S23*W32+S21*W12+S22*W22+S23* W32) - cmu3*(S12*S12+S22*S22+S23*S23-1./3.*S) );

C_UDSI(c,t,R_33)= -V*C_UDSI(c,t,K)*SQR(T)*
( cmu2*(S31*W13+S32*W23+S33*W33+S31*W13+S32*W23+S33* W33) - cmu3*(S13*S13+S23*S23+S33*S33-1./3.*S) );

C_UDSI(c,t,R_12)= -V*C_UDSI(c,t,K)*SQR(T)*
( cmu2*(S11*W12+S12*W22+S13*W32+S21*W11+S22*W21+S23* W31) - cmu3*(S11*S12+S12*S22+S13*S32) );



C_K(c,t)=C_UDSI(c,t,K);
C_D(c,t)=C_UDSI(c,t,E);

P11=2.0*C_R(c,t)*(-C_UDSI(c,t,R_11)*C_DUDX(c,t)-C_UDSI(c,t,R_12)*C_DUDY(c,t));
P22=2.0*C_R(c,t)*(-C_UDSI(c,t,R_12)*C_DVDX(c,t)-C_UDSI(c,t,R_22)*C_DVDY(c,t));
P33=0.;
Pkk=0.5*(P11+P22+P33);

C_UDMI(c,t,0)=mu_t;
C_UDMI(c,t,1)=T;
C_UDMI(c,t,2)=L;
C_UDMI(c,t,3)=mu_t*SQR(Strainrate_Mag(c,t))+Pkk;

}
end_c_loop(c,t)
}
}
}



DEFINE_TURBULENT_VISCOSITY(user_mu_t,c,t)
{

return C_UDMI(c,t,0);
}



DEFINE_PROFILE(e_bc, t, position)
{
real dy;

face_t f;
cell_t c0;
Thread *t0=t->t0;
real xw[ND_ND], xc[ND_ND], dx[ND_ND];

begin_f_loop(f,t)
{
c0=F_C0(f,t);
F_CENTROID(xw,f,t);
C_CENTROID(xc,c0,t0);
NV_VV(dx, =, xc, -, xw);
dy=ND_MAG(dx[0], dx[1], dx[2]);

F_PROFILE(f,t,position)=2.*C_MU_L(c0,t0)/C_R(c0,t0)*C_UDSI(c0,t0,K)/SQR(dy);

}
end_f_loop(f,t)
}



DEFINE_ON_DEMAND(on_demand_calc)
{
Domain *d; /* declare domain pointer since it is not passed as an argument to the DEFINE macro */

Thread *t;
cell_t c;

d = Get_Domain(1); /* Get the domain using Fluent utility */

thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
C_UDSI(c,t,K) = C_K(c,t);
C_UDSI(c,t,E) = C_D(c,t);
C_UDSI(c,t,V2) = 2.0/3.0*C_K(c,t);
C_UDSI(c,t,F) = 0.01;

}
end_c_loop(c,t)
}
}

[ghazanfarzahedi]

Hi Tstorm

I used a udf very similar to your udf (written by Dr. Heschl) in serial for a 2D flat plate.
It could run in parallel too without. (I checked it)
If you need it and its settings please mail to me:
ghazanfar.zahedi@gmail.com

Regards
Ghazanfar

[krisp]

Unless and until you have customized your UDF to write data you need not bother about modifying the code. On a quick look, I do not see such functions in your routine and you are good to go.

[friends.sk]

Hi Guys, I too stucked up parallezing my udf especially in the source term. Can someone help me in parallezing my udf.

Thanks.
-SK