[Siba11]

Hello.

I'm simulating the flow of blood through a stenosed artery. I'm treating it as a single-phase continuum flow. But I want to track the trajectories of all the hypothetical blood corpuscles that will travel in the artery (along the flow streamlines). For this, I have assumed that all the (hypothetical) particles will enter the domain from the inlet. Also, their positions at each iteration will update according to the Euler explicit formula. (I have attached a small 2-page PDF to better explain what I want to achieve).

Below is the UDF that I wrote for it:
Here, N_ITER is a macro that gives the iteration number.
p = is the number of particles.
T = total number of iterations.

#include "udf.h"
#include "cxndsearch.h"
#include "surf.h"
static ND_Search *domain_table = NULL;


 DEFINE_ADJUST (particle_tracking, d)
{
	#if !RP_HOST           
	
	CX_Cell_Id *cx_cell;   
	real P[2]={0,0};
	face_t f;
	cell_t c;
	Thread *t;
	FILE *str;
	real cen[ND_ND];
	int zoneid = 7; // inlet ID
	int i, p=0; 
	double h=0.001; 
	int T = 1500; 
	Thread *f_thread = Lookup_Thread(d,zoneid);

//-------Calculating number of inlet faces-------------
	
	begin_f_loop(f,f_thread)
	{
		p = p + 1; 
	}
	end_f_loop(f,f_thread)

//-----------------------------------------------------	
	
//----Declaring arrays to store particle trajectories---------	

	real x[T][p], y[T][p]; 
	
//---------------------------------------------------------------	
	if (N_ITER == 0)  // this 'if' block gives us the inlet face centroid coordinates.
	{
		begin_f_loop(f,f_thread)
		{
		F_CENTROID(cen,f,f_thread);
		x[N_ITER][f] = cen[0];
		y[N_ITER][f] = cen[1];
		}
	        end_f_loop(f,f_thread)
	}
//------------------------------------------------------------------	
	if (N_ITER == 1)  // this 'if' block is for the initial particle projection (shown by blue arrow in attached figure)
	{
		begin_f_loop(f,f_thread)
		{
		x[N_ITER][f] = x[N_ITER-1][f] + F_U(f,f_thread)*h;
		y[N_ITER][f] = y[N_ITER-1][f] + F_V(f,f_thread)*h;
		}
		end_f_loop(f,f_thread)
	}	
	
//-------------------------------------------------------------------
	
if (N_ITER>1)
{
	for (i=0;i<p;i++)
	{
		P[0] = x[N_ITER-1][i];
		P[1] = y[N_ITER-1][i];
		domain_table = CX_Start_ND_Point_Search(domain_table,TRUE,-1);
		cx_cell = CX_Find_Cell_With_Point(domain_table,P,0.0);
		c = RP_CELL(cx_cell);
		t = RP_THREAD(cx_cell);
		domain_table = CX_End_ND_Point_Search(domain_table);
		x[N_ITER][i] = x[N_ITER-1][i] + C_U(c,t)*h;
		y[N_ITER][i] = y[N_ITER-1][i] + C_V(c,t)*h;
	}
}
//------------Writing to file-----------------------------------	
			
		for (i=0;i<p;i++)
		{
			str = fopen("Output.txt","a");
			if(str==NULL)
			Message("error\n");
			fprintf(str,"%g,%g,",x[N_ITER][i],y[N_ITER][i]);
			fclose(str);
		}
		str = fopen("Output.txt","a");
		fprintf(str,"\n");
		fclose(str);
//----------------------------------------------------------	
#endif	
}

This code works fine for about T/2 iterations. Then it suddenly crashes, giving SIGSEGV (segmentation fault).

An observation that I made:
For T = 400, it crashes at 144th iteration.
For T = 500, it crashes at 244th iteration.
For T = 1500, it crashes at 819th iteration.

As far as I know, SIGSEGV occurs when there is some faulty memory allocation. But, I don't know what's causing it to crash midway through the solution.

I'll appreciate any help!