Homec4science
Diffusion Gear (master) Browse

Gear/src/vanishing.c
f97dd21ff912master

Last Change
Raw File

vanishing.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
#ifdef VANISHING_PARTICLES
/*
* This routine compute the accretion of gas
*/
void vanishing_particles(void)
{
if (ThisTask==0)
printf("Start vanishing particles...\n");
/* flag the particles to be removed */
vanishing_particles_flag();
/* remove the flaged particles */
vanishing_particles_remove();
if (ThisTask==0)
printf("vanishing particles done.\n");
}
/*
* This routine flag the particles that needs
* to be removed.
* This routine can in principle be called anywere in the code.
*/
void vanishing_particles_flag(void)
{
int i;
for(i = 0; i < NumPart; i++)
{
/* here we flag the particles according to a condition */
if (P[i].Pos[1]>49.5)
{
P[i].VanishingFlag=1;
//printf("--> vanishing %d\n",i);
}
else
P[i].VanishingFlag=0;
}
}
int compact_StP(int oldN_stars)
{
int i,j;
int N_StP_removed=0;
for(i = 0, j = oldN_stars-1; i <oldN_stars; i++)
{
if (StP[i].PIdx==-1)
{
/* skip last particles that must be removed too */
while (StP[j].PIdx==-1)
{
N_StP_removed++;
j--;
}
if (i>=j)
{
break;
}
/* replace the particle */
P[StP[j].PIdx].StPIdx = i;
StP[i] = StP[j];
N_StP_removed++;
j--;
}
}
printf("N_StP_removed=%d\n",N_StP_removed);
return N_StP_removed;
}
/*
* This routine remove the particles that have been flagged
* to be removed.
*/
void vanishing_particles_remove(void)
{
int DomainDecompFlag,DomainDecompFlagMax;
int i,j;
int N_gas_removed=0;
int N_stars_removed=0;
int N_other_removed=0;
int NumPart_removed=0;
int oldN_stars;
oldN_stars=N_stars;
/* be sure that all particles are in the right block */
rearrange_particle_sequence();
//printf("NumPart=%d\n",NumPart);
//printf("N_gas=%d\n",N_gas);
/******************************************/
/* remove block 0
/******************************************/
for(i = 0, j = N_gas-1; i < N_gas; i++)
{
if (P[i].VanishingFlag)
{
/* skip last particles that must be removed too */
while (P[j].VanishingFlag)
{
N_gas_removed++;
j--;
}
if (i>=j)
{
break;
}
/* replace the particle */
if(P[i].Ti_endstep == All.Ti_Current) /* if the particle was active, we need to decrease NumForceUpdate */
NumForceUpdate--;
P[i] = P[j];
SphP[i] = SphP[j];
N_gas_removed++;
j--;
}
}
/* shift 1 */
for(i = N_gas-N_gas_removed,j = (N_gas+N_stars)-1; i <N_gas ; i++,j--)
{
P[i] = P[j];
StP[P[i].StPIdx].PIdx=i;
}
/* shift others */
for(i = N_gas-N_gas_removed+N_stars,j = NumPart-1; i <(N_gas-N_gas_removed+N_stars)+N_gas_removed ; i++,j--)
{
P[i] = P[j];
}
N_gas=N_gas-N_gas_removed;
NumPart=NumPart-N_gas_removed;
/******************************************/
/* remove 1
/******************************************/
for(i = N_gas,j = N_gas+N_stars-1; i < N_gas+N_stars; i++)
{
if (P[i].VanishingFlag)
{
/* skip last particles that must be removed too */
while (P[j].VanishingFlag)
{
N_stars_removed++;
j--;
}
if (i>=j)
{
break;
}
/* replace the particle */
if(P[i].Ti_endstep == All.Ti_Current) /* if the particle was active, we need to decrease NumForceUpdate */
NumForceUpdate--;
StP[P[i].StPIdx].PIdx=-1; /* set it as free */
P[i] = P[j];
StP[P[i].StPIdx].PIdx=i;
N_stars_removed++;
j--;
}
}
/* shift others */
for(i = N_gas+N_stars-N_stars_removed,j = NumPart-1; i < N_gas+N_stars ; i++,j--)
{
P[i] = P[j];
}
N_stars=N_stars-N_stars_removed;
NumPart=NumPart-N_stars_removed;
/******************************************/
/* remove others
/******************************************/
for(i = N_gas+N_stars,j = NumPart-1; i < NumPart; i++)
{
if (P[i].VanishingFlag)
{
/* skip last particles that must be removed too */
while (P[j].VanishingFlag)
{
//printf("skip j=%2d %2d\n",j,P[j]);
N_other_removed++;
j--;
}
if (i>=j)
{
break;
}
/* replace the particle */
if(P[i].Ti_endstep == All.Ti_Current) /* if the particle was active, we need to decrease NumForceUpdate */
NumForceUpdate--;
P[i] = P[j];
N_other_removed++;
j--;
}
}
NumPart=NumPart-N_other_removed;
compact_StP(oldN_stars);
NumPart_removed = N_gas_removed+N_stars_removed+N_other_removed;
//printf("N_gas_removed =%d\n",N_gas_removed);
//printf("N_stars_removed=%d\n",N_stars_removed);
//printf("N_other_removed=%d\n",N_other_removed);
//printf("NumPart_removed=%d\n",NumPart_removed);
if (NumPart_removed>0)
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 */
for(i = 0; i < N_gas; i++)
{
//printf("%07d %d %g %g %g %g\n",i,P[i].Type,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],SphP[i].Hsml);
if (P[i].Type!=0)
{
endrun(999888777);
}
}
for(i = N_gas; i < NumPart; i++)
{
//printf("%07d %d %g %g %g\n",i,P[i].Type,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2]);
if (P[i].Type==0)
{
endrun(999888778);
}
}
//printf("NumPart=%d\n",NumPart);
//printf("N_gas=%d\n",N_gas);
}
#endif
c4science · Help