[brucelqs]

I encountered a problem while parallelizing the UDF code. For the following code, the results of using serial (1 CPU thread) calculation and parallel (15 CPU threads) calculation are different. The following is the original code (no parallelization operation):

#include "udf.h"

double CLP=0;
double T=0;
double TER=(110.0*M_PI/180.0);
double DCA=110.0;
int FT=1;
double CV=0;
FILE *file;
double S=0.0,R;


DEFINE_ADJUST(CAP, domain)
{
	Thread *thread = Lookup_Thread(domain, 1);
	Thread **pt = THREAD_SUB_THREADS(thread); 
	cell_t cell;
	face_t f;
	real x[ND_ND];
	double MX=0, FHI, XH, TE, Ca, VO;
	int n;

	S=0.0;
	
	begin_c_loop_all (cell,pt[1])
	{
		if(C_VOF(cell,pt[1])!=0)
		{
			C_CENTROID(x,cell,pt[1]);
			
			if(x[0]>MX)
			MX=x[0];
		}
		S+=C_VOF(cell,pt[1]);
	}
	end_c_loop_all (cell,pt[1])

	R=sqrt(S*1.0e-08/M_PI);
	printf("MX=%f R=%f CV=%f T=%f DCA=%f S=%f\n",MX,R,CV,T,DCA,S);
 
	if(FT==0)
	{
		CV= (R-CLP)/(RP_Get_Real("flow-time")-T);
		CLP = R;
		T=RP_Get_Real("flow-time"); 
		Ca = CV*1e-3/0.0728;
		TE= 0.5-0.5 * cos(TER);
		TE= 0.5 * log( (1.0+TE)/(1.0-TE) );
		FHI = -(9.78 * pow( TE,1.41))/(12.81 * pow(TE,1.41)-100.0);
		XH= Ca +  FHI ;
		if(CV>=0)
			DCA= acos( 1-2*tanh(5.16*pow((XH/(1+1.31*pow(XH,.99))),0.706)   ) ) *180.0/M_PI  ;
		else
			DCA= 2*TER*180.0/M_PI -acos( 1-2*tanh(5.16*pow((XH/(1+1.31*pow(XH,.99))),0.706)   ) ) *180.0/M_PI  ;
		file = fopen("file.txt", "a+");
		fprintf(file,"MX=%f R=%f CV=%f T=%f DCA=%f\n",MX,R,CV,T,DCA);
		fclose(file);
		
	 if(R>1e-18 && FT==1)
	 {
		CLP = R;
		T=RP_Get_Real("flow-time");
		FT=0;
		file = fopen("file.txt", "a+");
		fprintf(file,"MX=%f R=%f CV=%f T=%f\n",MX,R,CV,T);
	 	fclose(file);
	}
}

DEFINE_PROFILE(CAS,t,i)
{
	face_t f;
	begin_f_loop(f,t)
	{
		F_PROFILE(f,t,i) = DCA;
    }
	end_f_loop(f,t)
}

[AlexanderZ]

#include "udf.h"

double contact_line_position = 0;
double time = 0;
double theta_e_radian = (110.0 * M_PI / 180.0);
double dynamic_contact_angle = 110.0;
int first_time = 1;
double contact_velocity = 0;
FILE* file;
double sum = 0.0, R;

DEFINE_ADJUST(Contact_Angle_Update, domain)
{
	Thread* thread = Lookup_Thread(domain, 9);
	Thread** pt = THREAD_SUB_THREADS(thread);
	cell_t cell;
	face_t f;
	real x[ND_ND];
	double max_x = 0, f_Hoff_inverse, x_hoff, temp, Ca, volume;
	int n;

	sum = 0.0;

	begin_c_loop_all(cell, pt[1])
	{
		if (C_VOF(cell, pt[1]) != 0)
		{
			C_CENTROID(x, cell, pt[1]);
			
			if (x[0] > max_x)
				max_x = x[0];
		}
		sum += C_VOF(cell, pt[1]);
	}
	end_c_loop_all(cell, pt[1])

	#if RP_NODE
		sum = PRF_GRSUM1(sum);
	#endif

	R = sqrt(sum * 1.0e-08 / M_PI);
	Message0("max_x=%f R=%f contact_velocity=%f time=%f dynamic_contact_angle=%f sum=%f\n", max_x, R, contact_velocity, time, dynamic_contact_angle, sum);

	if (first_time == 0)
	{
		contact_velocity = (R - contact_line_position) / (RP_Get_Real("flow-time") - time);
		contact_line_position = R;
		time = CURRENT_TIME;
		Ca = contact_velocity * 1e-3 / 0.0728;
		temp = 0.5 - 0.5 * cos(theta_e_radian);
		temp = 0.5 * log((1.0 + temp) / (1.0 - temp));
		f_Hoff_inverse = -(9.78546 * pow(temp, 1.416430594900850)) / (12.819 * pow(temp, 1.41643) - 100.0);
		x_hoff = Ca + f_Hoff_inverse;
		if (contact_velocity >= 0)
			dynamic_contact_angle = acos(1. - 2. * tanh(5.16 * pow((x_hoff / (1. + 1.31 * pow(x_hoff, .99))), 0.706))) * 180.0 / M_PI;
		else
			dynamic_contact_angle = 2. * theta_e_radian * 180.0 / M_PI - acos(1. - 2. * tanh(5.16 * pow((x_hoff / (1. + 1.31 * pow(x_hoff, .99))), 0.706))) * 180.0 / M_PI;
	
	node_to_host_real_5(max_x, R, contact_velocity, time, dynamic_contact_angle);
	
	#if RP_HOST
			file = fopen("file.txt", "a+");
			fprintf(file, "max_x=%f R=%f contact_velocity=%f time=%f dynamic_contact_angle=%f\n", max_x, R, contact_velocity, time, dynamic_contact_angle);
			fclose(file);
	#endif
		if (R > 1e-18 && first_time == 1)
		{
			contact_line_position = R;
			time = CURRENT_TIME;
			first_time = 0;
	#if RP_HOST
				file = fopen("file.txt", "a+");
				fprintf(file, "max_x=%f R=%f contact_velocity=%f time=%f\n", max_x, R, contact_velocity, time);
				fclose(file);
	#endif
		}
	}
}

DEFINE_PROFILE(Contact_Angle_Set_Profile, t, i)
{
	face_t f;
	begin_f_loop(f, t)
	{
		F_PROFILE(f, t, i) = dynamic_contact_angle;
	}
	end_f_loop(f, t)
}

[brucelqs]

Thank you very much for your reply. I tried this code and found that in the case of 1 CPU thread (serial) and 16 CPU threads (parallel), the results obtained in parallel will diverge and terminate in the later stage of the calculation. This problem should not be a UDF code problem, the specific reason remains to be resolved.