[NonStopEagle]

Hey guys
I want to know if there is a way to calculate velocity gradients without using macros like C_U_G(),

I have been trying to implement maxwell First order slip condition to a wall of a supersonic micro-nozzle. The code I have works fine initially but diverges after about a 1000 iterations. I am guessing that this happens because of the gradients. So if there is any way of calculating the gradients manually please let me know. The udf is as follows:
#include "udf.h" /* must be at the beginning of every UDF */
/*#include "mem.h" /* must be at the beginning of this UDF */
//#include <math.h> /* must be at the beginning of this UDF */
/*#include "metric.h" /* must be at the beginning of this UDF */
#define k_B 1.3805e-23 /* Boltzman constant (j/K)*/
#define sigma 419e-12 /* molecular diameter*/
#define sigma_v 1 /* Tangential momentum accommodation coefficient*/
#define sigma_T 1 /*Thermal accommodation coefficient*/
#define pi 3.14159 /* pi number*/
#define sqrt_2 1.41421 /* sqrt(2) number */
#define len = 20e-5;
double relax = 0.1;
int counter = 0;
//#define relax 50e-3
DEFINE_PROFILE(wall_slip_velocity, t, index)
{
double T, P, density, MU, K, landa, a_m, du_da, dT_dx, u_f, DUY, DVY, DUX, DVX, u, v;
double A[ND_ND];
face_t f;
Thread* t0;
cell_t c0;
double DUR, DVR;


begin_f_loop(f, t)
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
T = C_T(c0, t0); /* temperature of the cell */
P = C_P(c0, t0); /* peressure of the cell */
density = C_R(c0, t0); /* density of the cell*/
MU = C_MU_L(c0, t0); /* laminar viscosity of the cell */
u = C_U(c0, t0);
v = C_V(c0, t0);
/*DUR = C_U_RG(c0, t0)[1];
DVR = C_V_RG(c0, t0)[1];*/

K = C_K_L(c0, t0); /* thermal conductivity of the cell */
landa = (k_B * T) / (sqrt_2 * pi * sigma * sigma * P); /* mean free path line */

DUX = C_DUDX(c0, t0);
DVX = C_DVDX(c0, t0);
DUY = C_DUDY(c0, t0);
DVY = C_DVDY(c0, t0);
F_AREA(A, f, t);
a_m = NV_MAG(A); /* magnitude of the cell area vector */
du_da = NV_DOT(A, A) / a_m; /* n- component of the cell x-velocity reconstruction gradient(RG) vector */
dT_dx = 1; /* x-component of the cell temperature reconstruction gradient(RG) vector */
//F_PROFILE(f, t, index) = relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
//F_UDMI(f,t,1)= relax * (((2 - sigma_v) / sigma_v) * landa * (u / sqrt(u * u + v * v) * (((DUX + DVX) * (v / sqrt(u * u + v * v))) + ((DUY + DVY) * (u / sqrt(u * u + v * v)))) + 0.75 * (MU / (density * T)) * dT_dx);
F_UDMI(f, t, 0) = relax * (((2 - sigma_v) / sigma_v) * landa * (((DUX + DVX) * (v / sqrt(u * u + v * v))) - ((DUY + DVY) * (u / sqrt(u * u + v * v)))));
// F_PROFILE(f, t, index) = F_UDMI(f, t, 0);
//F_PROFILE(f, t, index) = F_UDSI(f, t, 0);
F_UDMI(f, t, 1) = landa;

}
end_f_loop(f, t)

}
DEFINE_PROFILE(temperature_jump, t, index)
{
face_t f;
cell_t c0;
Thread* t0;
double p, tem, g, lamda, num_den, temp, vis, temp1, x, i, mag, DT2;
double DT, prd, tc, tj, cp, a,u,v;
double QUANT;
g = 9.81;
a = 4.29866e-19;
begin_f_loop(f, t) // for each face: get face-id and thread
{
t0 = THREAD_T0(t); /* adjacent cell thread to f */
c0 = F_C0(f, t); /* adjacent cell to f */
p = C_P(c0, t0); // getting the pressure
tem = C_T(c0, t0); // getting the temperature
tc = C_K_L(c0, t0);
vis = C_MU_L(c0, t0);
cp = C_CP(c0, t0);
u = C_U(c0, t0);
v = C_V(c0, t0);
DT = C_T_G(c0, t0)[1];
DT2 = C_T_G(c0, t0)[0];
prd = cp * vis / (g * tc);

num_den = p / ((1.38064852e-23) * tem); // Calculation of the number density
temp1 = num_den * a;
lamda = 1 / (temp1);
QUANT = (((DT2) * (v / sqrt(u * u + v * v))) - ((DT) * (u / sqrt(u * u + v * v))));
tj = 380+relax *((((2 - sigma_v) / sigma_v) *1.14566*lamda * QUANT ));//0.15596 (380 + ((1.14566 * 1 * (lamda / prd) * ((0.9877 * DT)))) * 0.9877);
F_PROFILE(f, t, index) = tj;


}
end_f_loop(f, t)
}
DEFINE_EXECUTE_AT_END(define_relax)
{
if ((counter % 200) == 0)
{
relax = relax + 0.01;
printf("Relaxation : %f", relax);
}
counter = counter + 1;
}