[Cooper24]

I want to loop over all the faces of all the cells along the z direction. I am using c_face_loop inside begin_c_loop. It's a two phase flow and I want to find the interfacial cells after every iteration. I am using INTERIOR_FACE_GEOMETRY to get the unit normal vector between the cell centroids. Please advise on the approach I am using. Here is my UDF:
mp_thread_loop_c(t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *ppt = pt[phase_domain_index];

begin_c_loop(c,t) // This will loop over all the cells in the domain
{
real xc[ND_ND];
C_CENTROID(xc,c,t);

int n;
int top;
C_UDMI(c,t,7) = xc[2];
C_UDMI(c,t,8) = C_VOF(c,ppt);

c_face_loop(c,t,n) // This will loop over all the faces of a cell and n is the local face index number
{
f = C_FACE(c,t,n);
tn = C_FACE_THREAD(c,t,n);
ca = F_C1(f,tn);
cb = F_C0(f,tn);
INTERIOR_FACE_GEOMETRY(f,tn,A,ds,es,A_by_es,dr0,dr 1);
ZDOT = NV_DOT(es,zunit);

if(ZDOT == 1 && C_VOF(cb,ppt)>0.05 && C_VOF(cb,ppt)< 1 && C_VOF(ca,ppt) == 0)
{
top = n;
ftop = C_FACE(c,t,top); // ftop is the global face index number. C_FACE converts local face index number, top to global face index number, ftop
tf = C_FACE_THREAD(c,t,top);
c1 = F_C1(ftop,tf); // c1 is the index of face's neighbouring C1 cell
c0 = F_C0(ftop,tf); // c0 is the index of face's neighbouring C0 cell
}
}

The UDF compiles but the simulation results are not right.

[AlexanderZ]

I think there is a logic problem in your algorithm

ftop = C_FACE(c,t,top); // ftop is the global face index number. C_FACE converts local face index number, top to global face index number, ftop

this statement is not true, as top = n;

where n is defined by the loop, where you are at this moment,
so

ftop = C_FACE(c,t,top); 
and 
f = C_FACE(c,t,n);

are equivalent

instead you could do

ftop = f
tf = tn
c0 = ca
c1 = cb

zunit is not defined in your code

ZDOT = NV_DOT(es,zunit);

and this code doesn't generate any data, so there is nothing to comment actually

[Cooper24]

Originally Posted by AlexanderZ

I think there is a logic problem in your algorithm

ftop = C_FACE(c,t,top); // ftop is the global face index number. C_FACE converts local face index number, top to global face index number, ftop

this statement is not true, as top = n;

where n is defined by the loop, where you are at this moment,
so

ftop = C_FACE(c,t,top); 
and 
f = C_FACE(c,t,n);

are equivalent

instead you could do

ftop = f
tf = tn
c0 = ca
c1 = cb

zunit is not defined in your code

ZDOT = NV_DOT(es,zunit);

and this code doesn't generate any data, so there is nothing to comment actually

Thanks for your reply. I have defined zunit in the previous lines which I forgot to mention. I will try to change the logical assignments as you have suggested. I wanted to ask about my approach. Is it a correct way to track the interfacial cells along the z direction?

[AlexanderZ]

as I said, the code you've showed doesn't make any output
so I can't comment is it the good way to organize this loop or not.

If I were you, I would make the loop over cell, check if it is the phase boundary but VOF criteria.
And if so, I would make a loop over faces here.

you start with faces and go vice versa, but I don't see problems here, It still could work.

[Cooper24]

Sorry that I didn't post the whole code earlier. Actually after identifying the interfacial cells with cell index c0 and c1 and storing them in UDMs, I am using them to apply source terms in those interfacial cells only as shown:

if(c == C_UDMI(c,t,9) && x[2]>0 && x[2]<0.08e-3)
{
if(r<=R)
{
source = (((2*eta*P)/(Pi*R*R))*exp((-2*(r*r))/(R*R)))*C_UDMI(c,t,3)*factor;
dS[eqn] = 0.0;
}
else
{
source = 0.0;
dS[eqn] = 0.0;
}
}
else

As for using the VOF criteria(between 0 and 1), from what I have learnt, it is partially correct way and I have seen that the convergence is not good with that. Now, my source is varying with time and it is a moving laser. Solid will melt to liquid and so the interface will continuously change and that's why I want to track those interfacial cells to apply the source terms.
The UDF compiles without any error and the simulation is also running but the results are completely different, infact wrong. From what I realised there might be some issue with order of lines in the UDF.
Here is the complete UDF:

DEFINE_ADJUST(adjust_interface, domain)
{
Thread *t;
Thread **pt;
cell_t c;
face_t f;
face_t fc;
face_t ftop;
cell_t c0;
cell_t c1;
cell_t ca;
cell_t cb;
Thread *tn;
Thread *tf;

int phase_domain_index = 1;

real zunit[ND_ND];
real ZDOT;

real es[ND_ND];
real ds;
real A[ND_ND];
real A_by_es;
real dr0[ND_ND];
real dr1[ND_ND];

NV_D(zunit, = ,0,0,1);


mp_thread_loop_c(t,domain,pt)
if (FLUID_THREAD_P(t))
{
Thread *ppt = pt[phase_domain_index];

begin_c_loop(c,t)
{
real xc[ND_ND];
C_CENTROID(xc,c,t);

int n;
int top;
C_UDMI(c,t,7) = xc[2];
C_UDMI(c,t,8) = C_VOF(c,ppt);

c_face_loop(c,t,n)
{
f = C_FACE(c,t,n);
tn = C_FACE_THREAD(c,t,n);
ca = F_C1(f,tn);
cb = F_C0(f,tn);
INTERIOR_FACE_GEOMETRY(f,tn,A,ds,es,A_by_es,dr0,dr 1);
ZDOT = NV_DOT(es,zunit);

if(ZDOT == 1 && C_VOF(cb,ppt)>0.05 && C_VOF(cb,ppt)< 1 && C_VOF(ca,ppt) == 0)
{
top = n;
ftop = C_FACE(c,t,top);
tf = C_FACE_THREAD(c,t,top);
c1 = F_C1(ftop,tf);
c0 = F_C0(ftop,tf);
}
}

/* This is the volume fraction condition */
if (C_VOF(c1,ppt) = 0.0 && C_VOF(c0,ppt)>0.05 && C_VOF(c0,ppt)< 1)
{
C_UDMI(c,t,9) = c0;
C_UDMI(c,t,10) = c1;
}
}
end_c_loop(c,t)
}
}

DEFINE_SOURCE(Laser, c, t, dS, eqn) // The name of the UDF is heat_source
{
Thread *pri_th;
Thread *sec_th;

real source;
real x[ND_ND], time; // Define face centroid vector, time
time = RP_Get_Real("flow-time"); // Acquire time from Fluent solver
C_CENTROID(x, c, t); // Acquire the cell centroid location
real T = C_T(c,t);

pri_th = THREAD_SUB_THREAD(t, 0);
sec_th = THREAD_SUB_THREAD(t, 1);

real alpha = C_VOF(c,sec_th); // cell volume fraction
real gamma;

if (T>293.0 && T<1658.0)
{
gamma = 0;
}
else if (T>=1658.0 && T<=1723.0)
{
gamma = (T-Ts)/(Tl-Ts);
}
else if (T>1723.0)
{
gamma = 1;
}


real Pv = 0.54*Pa*exp((Lv*M*(T-Tv))/(Rg*T*Tv));
real mv = (0.82*M*(Pv/0.54))/(sqrt(2*Pi*M*Rg*T));

real rhog = 1.662; // density of argon
real Cpg = 520.64; // specific heat of argon
real rhos = 7900; // density of solid SS316
real rhol = 7433 + 0.0393*T - 0.00018*pow(T,2); // density of liquid SS316
real rhom = rhol*gamma + rhos*(1-gamma); // density of SS316
real rho = alpha*rhom + rhog*(1-alpha); // density of cell containing metal and gas
real Cps = 462 + 0.134*T; // specific heat of solid SS316
real Cpl = 775; // specific heat of liquid SS316
real Cpm = Cpl*gamma + Cps*(1-gamma); // specific heat of SS316
real Cpn = alpha*Cpm + Cpg*(1-alpha); // specific heat of cell containing metal and gas

real factor = (2*rho*Cpn)/(rhom*Cpm + rhog*Cpg);

real r = sqrt(pow(x[0]-x0-v*time,2.0) + pow(x[1]-y0,2.0));

if(c == C_UDMI(c,t,9) && x[2]>0 && x[2]<0.08e-3)
{
if(r<=R)
{
source = (((2*eta*P)/(Pi*R*R))*exp((-2*(r*r))/(R*R)))*C_UDMI(c,t,3)*factor;
dS[eqn] = 0.0;
}
else
{
source = 0.0;
dS[eqn] = 0.0;
}
}
else
{
source = 0.0;
dS[eqn] = 0.0;
}
return source;
}

[AlexanderZ]

I think you may make you code more here:
was

if(ZDOT == 1 && C_VOF(cb,ppt)>0.05 && C_VOF(cb,ppt)< 1 && C_VOF(ca,ppt) == 0)
{
top = n;
ftop = C_FACE(c,t,top);
tf = C_FACE_THREAD(c,t,top);
c1 = F_C1(ftop,tf);
c0 = F_C0(ftop,tf);
}
}

/* This is the volume fraction condition */
if (C_VOF(c1,ppt) = 0.0 && C_VOF(c0,ppt)>0.05 && C_VOF(c0,ppt)< 1)
{
C_UDMI(c,t,9) = c0;
C_UDMI(c,t,10) = c1;
}

to be

if(ZDOT == 1 && C_VOF(cb,ppt)>0.05 && C_VOF(cb,ppt)< 1 && C_VOF(ca,ppt) == 0)
{
C_UDMI(c,t,9) = cb;
C_UDMI(c,t,10) = ca;
}

plot your C_UDMI(c,t,9) and check if it meets your expectations or not.

you've still provided not the full code, so I assume, that all parameters of laser is defined somewhere else and it's not a problem.

so start with plotting UDMI. If UDMs are correct, so the problem is not here.

I agree with the whole logic of your code. That should work.

[Cooper24]

Thanks a lot for your suggestions. I will pot C_UDMI(c,t,9) and see if I get the desired results. As I said that I am not an expert on writing UDFs and therefore was not confident enough with the logic.