gas_accretion.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Oct 20, 09:38

gas_accretion.c
View Options

This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.

	#include <stdio.h>
	#include <stdlib.h>
	#include <string.h>
	#include <math.h>
	#include <mpi.h>

	#include "allvars.h"
	#include "proto.h"


	#ifdef GAS_ACCRETION

	static struct mass_data
	{
	FLOAT key;
	int index;
	}
	*mp;

	static int *Id;

	/*! This is a comparison kernel for a sort routine, which is used to group
	* particles that are going to be exported to the same CPU.
	*/
	//int compare_mass(const void a, const void b)
	//{
	// if(((struct particle_data ) a)->Mass < (((struct particle_data ) b)->Mass))
	// return +1;
	// if(((struct particle_data ) a)->Mass > (((struct particle_data ) b)->Mass))
	// return -1;
	// return 0;
	//}

	int compare_mass(const void a, const void b)
	{

	/* perfom the comparison, only if particles have negative masses */

	if((((struct mass_data ) a)->key <= 0) && ((((struct mass_data ) b)->key <= 0 )))
	{

	if(((struct mass_data ) a)->key < (((struct mass_data ) b)->key))
	return +1;

	if(((struct mass_data ) a)->key > (((struct mass_data ) b)->key))
	return -1;
	}


	if((((struct mass_data ) a)->key <= 0) && ((((struct mass_data ) b)->key > 0 )))
	{
	return -1;
	}

	if((((struct mass_data ) a)->key > 0) && ((((struct mass_data ) b)->key <= 0 )))
	{
	return +1;
	}



	return 0;
	}




	/*! This function brings the gas particles into the same order as the sorted
	* keys. (The sort is first done only on the keys themselves and done
	* directly on the gas particles in order to reduce the amount of data that
	* needs to be moved in memory. Only once the order is established, the gas
	* particles are rearranged, such that each particle has to be moved at most
	* once.)
	*/
	void reorder_gas_using_mass(void)
	{
	int i;
	struct particle_data Psave, Psource;
	struct sph_particle_data SphPsave, SphPsource;

	int idsource, idsave, dest;

	for(i = 0; i < N_gas; i++)
	{
	if(Id[i] != i)
	{
	Psource = P[i];
	SphPsource = SphP[i];

	idsource = Id[i];
	dest = Id[i];

	do
	{
	Psave = P[dest];
	SphPsave = SphP[dest];
	idsave = Id[dest];

	P[dest] = Psource;
	SphP[dest] = SphPsource;
	Id[dest] = idsource;

	if(dest == i)
	break;

	Psource = Psave;
	SphPsource = SphPsave;
	idsource = idsave;

	dest = idsource;
	}
	while(1);
	}
	}
	}

	#ifdef STELLAR_PROP
	/*! This function brings the stars into the same order as
	* the sorted keys. (The sort is first done only on the keys themselves and
	* done directly on the particles in order to reduce the amount of data that
	* needs to be moved in memory. Only once the order is established, the
	* particles are rearranged, such that each particle has to be moved at most
	* once.)
	*/
	void reorder_stars_using_mass(void)
	{
	int i;
	struct particle_data Psave, Psource;
	int idsource, idsave, dest;

	for(i = N_gas; i < N_gas+N_stars; i++)
	{
	if(Id[i] != i)
	{
	Psource = P[i];
	idsource = Id[i];

	dest = Id[i];

	do
	{
	Psave = P[dest];
	idsave = Id[dest];

	P[dest] = Psource;
	Id[dest] = idsource;

	/* restore the link with Stp */
	StP[ P[dest].StPIdx ].PIdx = dest;

	if(dest == i)
	break;

	Psource = Psave;
	idsource = idsave;

	dest = idsource;
	}
	while(1);
	}
	}
	}


	void reorder_st_using_mass(void)
	{
	int i;
	struct st_particle_data StPsave, StPsource;
	int idsource, idsave, dest;

	for(i = 0; i < N_stars; i++)
	{
	if(Id[i] != i)
	{
	StPsource = StP[i];
	idsource = Id[i];
	dest = Id[i];

	do
	{
	StPsave = StP[dest];
	idsave = Id[dest];
	StP[dest] = StPsource;
	Id[dest] = idsource;

	/* restore the link with P */
	P[ StP[dest].PIdx ].StPIdx = dest;
	if(dest == i)
	break;

	StPsource = StPsave;
	idsource = idsave;
	dest = idsource;
	}
	while(1);
	}
	}

	}
	#endif


	/*! This function brings the collisionless particles into the same order as
	* the sorted keys. (The sort is first done only on the keys themselves and
	* done directly on the particles in order to reduce the amount of data that
	* needs to be moved in memory. Only once the order is established, the
	* particles are rearranged, such that each particle has to be moved at most
	* once.)
	*/
	void reorder_particles_using_mass(void)
	{
	int i;
	struct particle_data Psave, Psource;
	int idsource, idsave, dest;

	#ifdef STELLAR_PROP
	for(i = N_gas+N_stars; i < NumPart; i++)
	#else
	for(i = N_gas; i < NumPart; i++)
	#endif
	{
	if(Id[i] != i)
	{
	Psource = P[i];
	idsource = Id[i];

	dest = Id[i];

	do
	{
	Psave = P[dest];
	idsave = Id[dest];

	P[dest] = Psource;
	Id[dest] = idsource;

	if(dest == i)
	break;

	Psource = Psave;
	idsource = idsave;

	dest = idsource;
	}
	while(1);
	}
	}
	}



















	/*! This function sorts particles
	* according to their mass. Since gas particles need to
	* stay at the beginning of
	* the particle list, they are sorted as a separate block.
	*/
	void order_particles_using_mass(void)
	{
	int i;


	if(N_gas)
	{
	mp = malloc(sizeof(struct mass_data) * N_gas);
	Id = malloc(sizeof(int) * N_gas);

	for(i = 0; i < N_gas; i++)
	{
	mp[i].index = i;
	mp[i].key = P[i].Mass;
	}

	qsort(mp, N_gas, sizeof(struct mass_data), compare_mass);


	for(i = 0; i < N_gas; i++)
	Id[mp[i].index] = i;

	reorder_gas_using_mass();

	free(Id);
	free(mp);
	}



	#ifdef STELLAR_PROP
	if(N_stars>0)
	{
	mp = malloc(sizeof(struct mass_data) * (N_stars));
	mp -= (N_gas);

	Id = malloc(sizeof(int) * (N_stars));
	Id -= (N_gas);

	for(i = N_gas; i < N_gas+N_stars; i++)
	{
	mp[i].index = i;
	mp[i].key = P[i].Mass;
	}

	qsort(mp + N_gas, N_stars, sizeof(struct mass_data), compare_mass);

	for(i = N_gas; i < N_gas+N_stars; i++)
	Id[mp[i].index] = i;

	reorder_stars_using_mass();

	Id += N_gas;
	free(Id);
	mp += N_gas;
	free(mp);
	}


	if(NumPart - N_gas - N_stars > 0)
	{
	mp = malloc(sizeof(struct mass_data) * (NumPart - N_gas - N_stars));
	mp -= (N_gas+N_stars);

	Id = malloc(sizeof(int) * (NumPart - N_gas - N_stars));
	Id -= (N_gas+N_stars);

	for(i = N_gas+N_stars; i < NumPart; i++)
	{
	mp[i].index = i;
	mp[i].key = P[i].Mass;
	}

	qsort(mp + N_gas+N_stars, NumPart - N_gas - N_stars, sizeof(struct mass_data), compare_mass);

	for(i = N_gas+N_stars; i < NumPart; i++)
	Id[mp[i].index] = i;

	reorder_particles_using_mass();

	Id += N_gas+N_stars;
	free(Id);

	mp += N_gas+N_stars;
	free(mp);
	}

	/* here, we need to initializate the link
	* between P and StP, as it is usually done only in init() further.
	*/
	#ifndef CHIMIE_INPUT_ALL
	endrun(1212008);
	#else
	for(i = N_gas; i < N_gas+N_stars; i++)
	{
	P[i].StPIdx = i-N_gas;
	StP[P[i].StPIdx].PIdx = i;
	#ifdef CHECK_ID_CORRESPONDENCE
	StP[P[i].StPIdx].ID = P[i].ID;
	#endif
	}
	#endif


	/*
	if(N_stars>0)
	{
	Id = malloc(sizeof(int) * (N_stars));

	for(i = N_gas; i < N_gas+N_stars; i++)
	{
	Id[P[i].StPIdx] = i - N_gas;
	}

	reorder_st_using_mass();
	free(Id);
	}
	*/


	#else



	if(NumPart - N_gas > 0)
	{
	mp = malloc(sizeof(struct mass_data) * (NumPart - N_gas));
	mp -= (N_gas);

	Id = malloc(sizeof(int) * (NumPart - N_gas));
	Id -= (N_gas);

	for(i = N_gas; i < NumPart; i++)
	{
	mp[i].index = i;
	mp[i].key = Key[i];
	}

	qsort(mp + N_gas, NumPart - N_gas, sizeof(struct mass_data), compare_mass);

	for(i = N_gas; i < NumPart; i++)
	Id[mp[i].index] = i;

	reorder_particles_using_mass();

	Id += N_gas;
	free(Id);
	mp += N_gas;
	free(mp);
	}


	#endif




	#ifdef CHECK_ID_CORRESPONDENCE

	if (ThisTask==0)
	printf("Check id correspondence after mass order...\n");

	for(i = N_gas; i < N_gas+N_stars; i++)
	{

	//printf("-> i=%d ID=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n",i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);

	if(P[i].ID != StP[i-N_gas].ID)
	{
	printf("\nP/StP ID correspondance error\n");
	printf("(%d) (in mass) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d StP[P[i].StPIdx].ID=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx,StP[P[i].StPIdx].ID);
	endrun(1212001);
	}


	if( StP[P[i].StPIdx].PIdx != i )
	{
	printf("\nP/StP correspondance error\n");
	printf("(%d) (in gas accretion) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
	endrun(1212001);
	}


	if(StP[P[i].StPIdx].ID != P[i].ID)
	{
	printf("\nP/StP correspondance error\n");
	printf("(%d) (in gas accretion) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
	endrun(1212002);
	}
	}

	if (ThisTask==0)
	printf("Check id correspondence after mass order...\n");

	#endif



	}





	void allocate_gas_accretion(void)
	{
	size_t bytes;
	double bytes_tot = 0;


	if(NumPart_acc > 0)
	{
	bytes_tot = 0;

	if(!(Acc = malloc(bytes = NumPart_acc * sizeof(struct acc_particle_data))))
	{
	printf("failed to allocate memory for `Acc' (%g MB) %d.\n", bytes / (1024.0 * 1024.0),sizeof(struct acc_particle_data));
	endrun(1);
	}
	bytes_tot += bytes;

	if(ThisTask == 0)
	printf("\nAllocated %g MByte for accretion storage. %d\n\n", bytes_tot / (1024.0 * 1024.0), sizeof(struct acc_particle_data));
	}



	if(N_gas_acc > 0)
	{
	bytes_tot = 0;

	if(!(SphAcc = malloc(bytes = N_gas_acc * sizeof(struct gas_acc_particle_data))))
	{
	printf("failed to allocate memory for `SphAcc' (%g MB) %d.\n", bytes / (1024.0 * 1024.0),sizeof(struct gas_acc_particle_data));
	endrun(1);
	}
	bytes_tot += bytes;

	if(ThisTask == 0)
	printf("\nAllocated %g MByte for gas accretion storage. %d\n\n", bytes_tot / (1024.0 * 1024.0), sizeof(struct gas_acc_particle_data));
	}





	}


	/*
	* This routine compute the accretion of gas
	*/
	void init_gas_accretion(void)
	{
	int i,j,k;
	int type;
	int *temp;
	double MaxMass[6];

	if(ThisTask == 0)
	{
	printf("Init gas accretion...\n");
	fflush(stdout);
	}



	/*

	As we do not have the value for the mass, we take
	the maximum of the particle class.

	*/

	for(i = 0, N_gas_acc = 0, NumPart_acc; i < NumPart; i++)
	{

	type = P[i].Type;
	MaxMass[type] = dmax(MaxMass[type],P[i].Mass);

	if (P[i].Mass <= 0)
	{
	NumPart_acc++;
	if (type==0)
	{
	N_gas_acc++;
	SphP[i].ActiveFlag=0;
	}

	#ifdef STELLAR_PROP
	if (type==1)
	{
	N_stars_acc++;
	}
	#endif



	}


	}

	//printf("(%d) NumPart_acc=%d N_gas_acc=%d\n",ThisTask,NumPart_acc,N_gas_acc);


	MPI_Allreduce(&MaxMass, &All.AccretionParticleMass, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

	if(ThisTask == 0)
	{
	printf("(%d) AccretionParticleMass[0]=%g\n",ThisTask,All.AccretionParticleMass[0]);
	printf("(%d) AccretionParticleMass[1]=%g\n",ThisTask,All.AccretionParticleMass[1]);
	printf("(%d) AccretionParticleMass[2]=%g\n",ThisTask,All.AccretionParticleMass[2]);
	printf("(%d) AccretionParticleMass[3]=%g\n",ThisTask,All.AccretionParticleMass[3]);
	printf("(%d) AccretionParticleMass[4]=%g\n",ThisTask,All.AccretionParticleMass[4]);
	printf("(%d) AccretionParticleMass[5]=%g\n",ThisTask,All.AccretionParticleMass[5]);
	fflush(stdout);
	}



	/* because ntot_type_all[] is of type `long long', we cannot do a simple
	* MPI_Allreduce() to sum the total particle numbers
	*/

	temp = malloc(NTask * sizeof(int));
	MPI_Allgather(&NumPart_acc, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);

	for(i = 0, All.TotNumPart_acc=0; i < NTask; i++)
	All.TotNumPart_acc += temp[i];


	temp = malloc(NTask * sizeof(int));
	MPI_Allgather(&N_gas_acc, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);

	for(i = 0, All.TotN_gas_acc=0; i < NTask; i++)
	All.TotN_gas_acc += temp[i];



	if(ThisTask == 0)
	{
	printf("Number of accretion particles found : %6d%09d \n",(int) (All.TotNumPart_acc / 1000000000), (int) (All.TotNumPart_acc % 1000000000));
	printf("Number of gas accretion particles found : %6d%09d \n",(int) (All.TotN_gas_acc / 1000000000), (int) (All.TotN_gas_acc % 1000000000));
	fflush(stdout);
	}


	allocate_gas_accretion();


	/*

	---> ok, jusque ici c'est bon


	pour la suite:

	1) on trie toutes les particules selon leur masse, par type...
	---> un peu embetant si on utilise Stp... mais bon, comme dans peano (peano_hilbert_order)

	-> est la suite est simple


	2) - on extrait les particules sans les trier
	- on fait une quantité de memshift (mais c'est pas trops grave)
	- on trie uniquement Acc et SphAcc ....

	*/


	/* sort gas particles according to their Mass (blocks are respected) */
	/* negative mass are first */
	order_particles_using_mass();


	if (ThisTask==0)
	printf("mass order done.\n");




	//for(i=0;i<NumPart;i++)
	// printf("(%d) i=%d id=%d type=%d m=%g \n",ThisTask,i,P[i].ID,P[i].Type,P[i].Mass);




	/* copy the particles to Acc and SphAcc*/
	for(i = 0, j = 0; i < NumPart; i++)
	{

	if (P[i].Mass <= 0)
	{

	Acc[j].Pos[0] = P[i].Pos[0];
	Acc[j].Pos[1] = P[i].Pos[1];
	Acc[j].Pos[2] = P[i].Pos[2];

	Acc[j].Vel[0] = P[i].Vel[0];
	Acc[j].Vel[1] = P[i].Vel[1];
	Acc[j].Vel[2] = P[i].Vel[2];

	Acc[j].Mass = All.AccretionParticleMass[P[i].Type];

	Acc[j].ID = P[i].ID;
	Acc[j].Type = P[i].Type;

	Acc[j].Time = -P[i].Mass;


	/* gas particle */
	if (P[i].Type==0)
	{
	SphAcc[j].Entropy = SphP[i].Entropy;

	#ifdef CHIMIE
	for(k=0;k<NELEMENTS;k++)
	SphAcc[j].Metal[k] = SphP[i].Metal[k];
	#endif
	}

	j++;

	}
	}


	//printf("test 0\n");
	//for(i=0;i<NumPart;i++)
	// printf("(%d) i=%d id=%d m=%g type=%d\n",ThisTask,i,P[i].ID,P[i].Mass,P[i].Type);


	/* shift gas */

	//memmove(&P[0],&P[N_gas_acc],(NumPart-(N_gas-N_gas_acc))*sizeof(struct particle_data));
	//memmove(&SphP[0],&SphP[N_gas_acc],(N_gas-(N_gas-N_gas_acc))*sizeof(struct sph_particle_data));
	memmove(&P[0],&P[N_gas_acc],(NumPart-N_gas_acc)*sizeof(struct particle_data));
	memmove(&SphP[0],&SphP[N_gas_acc],(N_gas-N_gas_acc)*sizeof(struct sph_particle_data));

	N_gas = N_gas-N_gas_acc;
	NumPart = NumPart-N_gas_acc;
	All.TotN_gas=All.TotN_gas - All.TotN_gas_acc;
	All.TotNumPart = All.TotNumPart - All.TotN_gas_acc;

	#ifdef STELLAR_PROP
	if (N_stars_acc>0)
	{
	printf("gas_accretion: this is not implemented\n");
	endrun(6612345);
	}
	/* shift stars */
	/* not implemented */


	/* shift particles */
	//memmove(&P[N_gas+N_stars],&P[N_gas+N_stars+(NumPart_acc-N_gas_acc)],(NumPart-(N_gas+N_stars+(NumPart_acc-N_gas_acc)))*sizeof(struct particle_data));
	memmove(&P[N_gas],&P[N_gas+(NumPart_acc-N_gas_acc-N_stars_acc)],(NumPart-(N_gas+(NumPart_acc-N_gas_acc-N_stars_acc)))*sizeof(struct particle_data));

	NumPart = NumPart-(NumPart_acc-N_gas_acc-N_stars_acc);
	All.TotNumPart = All.TotNumPart - (All.TotNumPart_acc-All.TotN_gas_acc);

	#else

	/* shift particles */
	memmove(&P[N_gas],&P[N_gas+(NumPart_acc-N_gas_acc)],(NumPart-(N_gas+(NumPart_acc-N_gas_acc)))*sizeof(struct particle_data));

	NumPart = NumPart-(NumPart_acc-N_gas_acc);
	All.TotNumPart = All.TotNumPart - (All.TotNumPart_acc-All.TotN_gas_acc);

	#endif



	//printf("test 1\n");
	//for(i=0;i<NumPart;i++)
	// printf("(%d) i=%d id=%d m=%g type=%d\n",ThisTask,i,P[i].ID,P[i].Mass,P[i].Type);



	for (i=0;i<NumPart;i++)
	if (P[i].Mass<=0)
	endrun(12453243);


	//for (i=0;i<NumPart_acc;i++)
	// {
	// printf("(%d) i=%06d ID=%d Accretion Time=%10.5g Type=%d xyz=(%g %g %g)\n ",ThisTask,i,Acc[i].ID,Acc[i].Time,Acc[i].Type,Acc[i].Pos[0],Acc[i].Pos[1],Acc[i].Pos[2]);
	// }

	if(ThisTask == 0)
	{
	printf("gas accretion init done.\n\n");
	fflush(stdout);
	}

	}









	/*
	* This routine compute the accretion of gas
	*/
	void gas_accretion(void)
	{
	int i,j,p;
	int np,ng;
	int accNumPart,accN_gas,accN_stars;
	int accTotNumPart,accTotN_gas,accTotN_stars;
	int *temp;
	float accM,accTotM;
	int DomainDecompFlag,DomainDecompFlagMax;



	if(ThisTask == 0)
	{
	printf("Start gas accretion...\n");
	fflush(stdout);
	}



	/* init */

	accNumPart=accN_gas=accN_stars=0;
	accTotNumPart=accTotN_gas=accTotN_stars=0;
	accM=accTotM=0;


	for(i = 0; i < NumPart_acc; i++)
	{
	if (Acc[i].Time<All.Time)
	{
	accNumPart++;
	accM+=Acc[i].Mass;
	printf(" (%d) time=%g acreting i=%d id=%d time=%g type=%d\n",ThisTask,All.Time,i,Acc[i].ID,Acc[i].Time,Acc[i].Type);


	if (Acc[i].Type==0)
	accN_gas++;


	#ifdef STELLAR_PROP
	if (Acc[i].Type==1)
	accN_stars++;
	#endif

	}

	//else
	// {
	// //printf("(%d) -- %g %g\n",ThisTask,All.Time,Acc[i].Time);
	// break;
	// }

	}


	temp = malloc(NTask * sizeof(int));
	MPI_Allgather(&accNumPart, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);

	for(i = 0; i < NTask; i++)
	accTotNumPart += temp[i];

	free(temp);


	temp = malloc(NTask * sizeof(int));
	MPI_Allgather(&accN_gas, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);

	for(i = 0; i < NTask; i++)
	accTotN_gas += temp[i];

	free(temp);

	#ifdef STELLAR_PROP
	temp = malloc(NTask * sizeof(int));
	MPI_Allgather(&accN_stars, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);

	for(i = 0; i < NTask; i++)
	accTotN_stars += temp[i];

	free(temp);
	#endif



	if(ThisTask == 0)
	{
	printf("Number of accretion particles this time : %9d \n",accTotNumPart);
	printf("Number of gas accretion particles this time : %9d \n",accTotN_gas);
	printf("Number of other accretion particles this time : %9d \n",accTotNumPart-accTotN_gas);

	fflush(stdout);
	}

	All.TotN_gas += accTotN_gas;
	All.TotNumPart += accTotNumPart;



	DomainDecompFlag=0;

	if (accNumPart>0)
	{

	if (NumPart + accNumPart > All.MaxPart)
	{
	printf("No memory left for new particle : NumPart + accNumPart=%d All.MaxPart=%d\n\n",NumPart + accNumPart,All.MaxPart);
	endrun(877001);
	}

	/* need to insert particles */

	if (accNumPart-accN_gas>0)
	{

	/* shift */
	#ifndef STELLAR_PROP
	memmove(&P[N_gas+(accNumPart-accN_gas)], &P[N_gas], (NumPart-N_gas) * sizeof(struct particle_data));
	NumPart = NumPart + (accNumPart-accN_gas);
	#else
	memmove(&P[N_gas+N_stars+(accNumPart-accN_gas)], &P[N_gas+N_stars], (NumPart-N_gas-N_stars) * sizeof(struct particle_data));
	NumPart = NumPart + (accNumPart-accN_gas);
	#endif
	}


	/* need to insert gas */
	if (accN_gas>0)
	{


	if (N_gas + accN_gas > All.MaxPartSph)
	{
	printf("No memory left for new gas particle : N_gas + accN_gas=%d All.MaxPartSph=%d\n\n",N_gas + accN_gas,All.MaxPartSph);
	endrun(877002);
	}

	memmove(&P[accN_gas], &P[0], (NumPart) * sizeof(struct particle_data));
	NumPart = NumPart + accN_gas;

	memmove(&SphP[accN_gas], &SphP[0], (N_gas) * sizeof(struct sph_particle_data));
	N_gas = N_gas + accN_gas;


	}



	#ifdef STELLAR_PROP
	/* recover the link between P and StP */
	for(i = N_gas; i < N_gas+N_stars; i++)
	StP[P[i].StPIdx].PIdx=i;
	#endif


	/*
	now, insert particles

	i = index of Acc and SphAcc
	p = index in P and SphP

	*/

	ng=0; /* index of first gaseous particle */

	#ifndef STELLAR_PROP
	np=N_gas; /* index of first other particle */
	#else
	np=N_gas+N_stars; /* index of first other particle */
	#endif



	for(i = 0; i < NumPart_acc; i++) /* particles are not in the correct order */
	//for (i=0;i<accNumPart;i++)
	{

	if (Acc[i].Time<All.Time)
	{


	if (Acc[i].Type==0)
	{
	p=ng;
	ng++;
	}
	else
	{
	p=np;
	np++;
	}


	P[p].ID = Acc[i].ID;

	for (j=0;j<3;j++)
	P[p].Pos[j] = Acc[i].Pos[j];

	for (j=0;j<3;j++)
	P[p].Vel[j] = Acc[i].Vel[j];

	P[p].Mass = Acc[i].Mass;
	P[p].Type = Acc[i].Type;

	/* force to be active */
	P[p].Ti_endstep=All.Ti_Current;


	#ifdef VANISHING_PARTICLES
	P[p].VanishingFlag=0;
	#endif



	if (Acc[i].Type==0)
	{
	SphP[p].ActiveFlag=1;

	SphP[p].Entropy = SphAcc[i].Entropy; /* initial internal energy*/
	SphP[p].Hsml=All.SofteningTable[P[p].Type]; /* this can be better */


	#ifdef CHIMIE
	for(j=0;j<NELEMENTS;j++)
	{
	SphP[p].Metal[j] = SphAcc[i].Metal[j];
	}
	#endif
	}


	//printf("(%d)ACCRETION : id=%d x=%g y=%g z=%g\n",ThisTask,P[p].ID,P[p].Pos[0],P[p].Pos[1],P[p].Pos[2]);


	}


	}


	/*
	finally, shift Acc and SphAcc
	*/


	if (accNumPart-accN_gas>0)
	{
	/* shift accreated particles */
	memmove(&Acc[N_gas_acc],&Acc[N_gas_acc+(accNumPart-accN_gas)],(NumPart_acc-N_gas_acc-(accNumPart-accN_gas))*sizeof(struct acc_particle_data));
	NumPart_acc = NumPart_acc - (accNumPart-accN_gas);
	}

	if (accN_gas>0)
	{
	/* shift accreated gas particles */
	memmove(&Acc[0],&Acc[accN_gas],(NumPart_acc-accN_gas)*sizeof(struct acc_particle_data));
	NumPart_acc = NumPart_acc - accN_gas;

	memmove(&SphAcc[0],&SphAcc[accN_gas],(N_gas_acc-accN_gas)*sizeof(struct gas_acc_particle_data));
	N_gas_acc = N_gas_acc - accN_gas;
	}




	/* force to update the domain and build the tree */
	DomainDecompFlag=1;


	}


	/* gather flags */
	MPI_Allreduce(&DomainDecompFlag, &DomainDecompFlagMax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);

	if (DomainDecompFlagMax>0)
	All.NumForcesSinceLastDomainDecomp = All.TotNumPart * All.TreeDomainUpdateFrequency + 1;

	/* some checks */


	#ifdef CHECK_ID_CORRESPONDENCE
	if (ThisTask==0)
	printf("Check id correspondence after memmove in gas_accretion...\n");

	rearrange_particle_sequence();

	for(i = N_gas; i < N_gas+N_stars; i++)
	{
	if( StP[P[i].StPIdx].PIdx != i )
	{
	printf("\nP/StP correspondance error\n");
	printf("(%d) (in domain before) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
	endrun(877003);
	}


	if(StP[P[i].StPIdx].ID != P[i].ID)
	{
	printf("\nP/StP correspondance error\n");
	printf("(%d) (in domain before) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
	endrun(877004);
	}
	}
	if (ThisTask==0)
	printf("Check id correspondence after memmove in gas_accretion done.\n");
	#endif



	#ifdef CHECK_BLOCK_ORDER
	int typenow;
	typenow = 0;
	for(i = 0; i < NumPart; i++)
	{

	if ((P[i].Type<typenow)&&(typenow<=ST))
	{
	printf("\nBlock order error\n");
	printf("(%d) i=%d id=%d Type=%d typenow=%d\n\n",ThisTask,i,P[i].ID,P[i].Type,typenow);
	endrun(877005);
	}

	typenow = P[i].Type;
	}
	#endif




	/* some info */


	if (ThisTask==0)
	{
	fprintf(FdGasAccretion, "%15g %15g %09d %g\n", All.Time,All.TimeStep, accTotNumPart,accTotM);
	fflush(FdGasAccretion);
	}

	if(ThisTask == 0)
	{
	printf("gas accretion done.\n");
	fflush(stdout);
	}


	printf("(%d) gas accretion done.\n",ThisTask);


	}


	/*
	* This routine compute the accretion of gas
	*/
	void update_entropy_for_accreated_particles(void)
	{
	int i;

	for (i=0;i<N_gas;i++)
	if (SphP[i].ActiveFlag)
	{
	//printf("------------------------ end of density id=%d rho=%g Hsml=%g entr=%g %g\n",P[i].ID,SphP[i].Density,SphP[i].Hsml,SphP[i].Entropy,P[i].Mass);
	SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / 1, GAMMA_MINUS1);
	//printf("------------------------ end of density id=%d rho=%g Hsml=%g entr=%g %g\n",P[i].ID,SphP[i].Density,SphP[i].Hsml,SphP[i].Entropy,P[i].Mass);
	SphP[i].ActiveFlag=0;
	}
	}


	#endif /* GAS_ACCRETION */

gas_accretion.cNo OneTemporaryActions

File Metadata

gas_accretion.cView Options

Event Timeline

gas_accretion.c
No OneTemporary
Actions

gas_accretion.c
View Options