Homec4science
Diffusion Gear (master) Browse

Gear/src/integral_conserving_dissipation.c
aefcf288a995master

Last Change
Raw File

integral_conserving_dissipation.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
#ifdef INTEGRAL_CONSERVING_DISSIPATION
#define DISSIPATIVE_TYPE 0
static double M,X[3],V[3],I,J,T;
static double Mtot,Xtot[3],Vtot[3],Itot,Jtot,Ttot;
static double alpha,DeltaT,DeltaTmin;
static int N_ngb;
static int *numlist;
static int LocalDissipation=1;
static int des_nngb=5;
static double xc[3];
static int Nmax_icd;
static struct icd_data
{
FLOAT key;
int index;
int proc;
}
*icdd,*icddGet;
static struct icd_particle_data
{
FLOAT Pos[3];
unsigned int ID;
int proc;
}
*Picd,*PicdGet;
/*! setup parameters for integral_conserving_dissipation
*
* LocalDissipation = 1, means that only des_nngb
* particles dissipates altogether
*
* LocalDissipation = 0, means that all particles dissipates
*
*/
int setup_integral_conserving_dissipation_local(int flag,int n,double x[3])
{
LocalDissipation = flag;
des_nngb = n;
xc[0] = x[0];
xc[1] = x[1];
xc[2] = x[2];
}
void icd_set_LocalDissipationOn(void)
{
LocalDissipation = 1;
}
void icd_set_LocalDissipationOff(void)
{
LocalDissipation = 0;
}
void icd_set_des_nngb(int n)
{
des_nngb = n;
}
void icd_set_xc(double x[3])
{
xc[0] = x[0];
xc[1] = x[1];
xc[2] = x[2];
}
/*! This function is the driver routine for the calculation of dissipation
*/
void integral_conserving_dissipation(void)
{
double Ekin1,Ekin2;
if(ThisTask == 0)
{
printf("Start integral conserving dissipation...\n");
fflush(stdout);
}
Ekin1 = get_Gas_energy_kin();
if (LocalDissipation)
local_ic_dissipation_m1();
else
global_ic_dissipation_m1();
Ekin2 = get_Gas_energy_kin();
if (ThisTask==0)
printf("\nEkin1Gas=%g Ekin2Gas=%g dE/E=%g\n\n", Ekin1,Ekin2,(Ekin2-Ekin1)/Ekin1 );
LocalSysState.EnergyICDissipation -= Ekin2-Ekin1;
if(ThisTask == 0)
{
printf("Integral conserving dissipation done.\n");
fflush(stdout);
}
}
/*!
*/
void global_ic_dissipation_m1(void)
{
compute_integrals_m1();
compute_alpha_m1();
apply_transformation_m1();
compute_energy_kin();
}
/*!
*/
void local_ic_dissipation_m1(void)
{
int i,j,n;
int N_icd,N_icdGet;
MPI_Status status;
Nmax_icd = (int)(0.1*N_gas);
Picd = malloc(Nmax_icd*sizeof(struct icd_particle_data));
PicdGet = malloc(Nmax_icd*sizeof(struct icd_particle_data));
N_icd = set_dissipation_list();
/* loop over all Tasks */
for (n=0;n<NTask;n++)
{
/* master get the list of particles from task n */
if (n==0)
{
N_icdGet = N_icd;
memcpy(PicdGet, Picd, N_icd*sizeof(struct icd_particle_data));
}
else
{
if (ThisTask==0)
{
MPI_Recv(&N_icdGet, 1, MPI_INT , n, TAG_N, MPI_COMM_WORLD, &status);
MPI_Recv(PicdGet, N_icdGet*sizeof(struct icd_particle_data), MPI_BYTE , n, TAG_N, MPI_COMM_WORLD, &status);
}
if (ThisTask==n)
{
MPI_Send(&N_icd, 1, MPI_INT, 0, TAG_N, MPI_COMM_WORLD);
MPI_Send(Picd, N_icd*sizeof(struct icd_particle_data), MPI_BYTE, 0, TAG_N, MPI_COMM_WORLD);
}
}
/* master send N_icdGet to all */
MPI_Bcast(&N_icdGet, 1, MPI_INT, 0, MPI_COMM_WORLD);
for(i=0;i<N_icdGet;i++)
{
//set_dissipation_center();
if (ThisTask==0)
{
for (j=0;j<3;j++)
xc[j] = PicdGet[i].Pos[j];
if (ThisTask==0)
printf("choosing particle %d (%g %g %g) from proc %d\n",PicdGet[i].ID,xc[0],xc[1],xc[2],PicdGet[i].proc);
}
/* send to all procs */
MPI_Bcast(&xc, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
icd_allocate_numlist();
icd_find_ngbs();
compute_integrals_m1();
compute_alpha_m1();
apply_transformation_m1();
//compute_energy_kin();
icd_free_numlist();
}
}
free(Picd);
free(PicdGet);
}
void icd_allocate_numlist()
{
numlist = malloc(des_nngb * sizeof(int)); /* index list of local particles belonging to neighbors */
N_ngb = 0;
}
void icd_free_numlist()
{
free(numlist);
}
/*
*/
int set_dissipation_list(void)
{
int i,j;
int N_icd;
int index;
N_icd = (int) (32./ NTask);
if (N_icd>Nmax_icd)
{
printf("N_icd (%d) > Nmax_icd (%d), you should increase Nmax_icd \n\n",N_icd,Nmax_icd);
endrun(10101088);
}
for (i=0;i<N_icd;i++)
{
index = (int) ( gsl_rng_uniform(random_generator) * (N_gas-1) );
for(j=0;j<3;j++)
Picd[i].Pos[j] = P[index].Pos[j];
Picd[i].proc = ThisTask;
Picd[i].ID = P[index].ID;
}
return N_icd;
}
/*
*/
void set_dissipation_center(void)
{
int i;
if (ThisTask==0)
i = (int) ( gsl_rng_uniform(random_generator) * (N_gas-1) );
xc[0] = P[i].Pos[0];
xc[1] = P[i].Pos[1];
xc[2] = P[i].Pos[2];
if (ThisTask==0)
printf("choosing particle %d (%g %g %g)\n",i,xc[0],xc[1],xc[2]);
/* send to all procs */
MPI_Bcast(&xc, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
/*!
*/
void compute_integrals_m1(void)
{
int i,j,k;
M=0;
T=0;
I=0;
J=0;
Mtot=0;
Ttot=0;
Itot=0;
Jtot=0;
for (j=0;j<3;j++)
{
X[j]=0;
V[j]=0;
Xtot[j]=0;
Vtot[j]=0;
}
if (LocalDissipation)
{
for (k=0;k<N_ngb;k++)
{
i = numlist[k];
M += P[i].Mass;
for (j=0;j<3;j++)
{
X[j] += P[i].Pos[j]*P[i].Mass;
V[j] += P[i].Vel[j]*P[i].Mass;
T += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
}
}
}
else
{
/* compute the (local) total mass M */
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == DISSIPATIVE_TYPE)
{
M += P[i].Mass;
for (j=0;j<3;j++)
{
X[j] += P[i].Pos[j]*P[i].Mass;
V[j] += P[i].Vel[j]*P[i].Mass;
T += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
}
}
}
}
T *= 0.5;
MPI_Allreduce(&M, &Mtot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&T, &Ttot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&X, &Xtot, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&V, &Vtot, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
/* compute X and V */
for (j=0;j<3;j++)
{
Xtot[j] /= Mtot;
Vtot[j] /= Mtot;
}
/* compute I and J */
if (LocalDissipation)
{
for (k=0;k<N_ngb;k++)
{
i = numlist[k];
I+= P[i].Mass * ( (P[i].Pos[0] - Xtot[0])*(P[i].Pos[0] - Xtot[0])
+ (P[i].Pos[1] - Xtot[1])*(P[i].Pos[1] - Xtot[1])
+ (P[i].Pos[2] - Xtot[2])*(P[i].Pos[2] - Xtot[2]) ) ;
J+= P[i].Mass * ( (P[i].Vel[0] - Vtot[0])*(P[i].Pos[0] - Xtot[0])
+ (P[i].Vel[1] - Vtot[1])*(P[i].Pos[1] - Xtot[1])
+ (P[i].Vel[2] - Vtot[2])*(P[i].Pos[2] - Xtot[2]) ) ;
}
}
else
{
/* compute the (local) total mass M */
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == DISSIPATIVE_TYPE)
{
I+= P[i].Mass * ( (P[i].Pos[0] - Xtot[0])*(P[i].Pos[0] - Xtot[0])
+ (P[i].Pos[1] - Xtot[1])*(P[i].Pos[1] - Xtot[1])
+ (P[i].Pos[2] - Xtot[2])*(P[i].Pos[2] - Xtot[2]) ) ;
J+= P[i].Mass * ( (P[i].Vel[0] - Vtot[0])*(P[i].Pos[0] - Xtot[0])
+ (P[i].Vel[1] - Vtot[1])*(P[i].Pos[1] - Xtot[1])
+ (P[i].Vel[2] - Vtot[2])*(P[i].Pos[2] - Xtot[2]) ) ;
}
}
}
J *=2;
MPI_Allreduce(&I, &Itot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&J, &Jtot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
DeltaTmin = -Jtot*Jtot/(8*Itot);
if (ThisTask==0)
printf("I=%g J=%g \n", Itot,Jtot );
}
/*!
*/
void compute_alpha_m1(void)
{
/* compute alpha */
double J2I;
double sgn;
DeltaT = -0.01*T;
J2I = 0.5*Jtot/Itot;
/*
if (DeltaT<=DeltaTmin)
{
alpha = (-J2I);
}
else
{
if (Jtot>0)
sgn=+1;
else
sgn=-1;
alpha = (-J2I) + sgn * pow( (J2I*J2I) + 2* DeltaT/Itot , 0.5) ;
}
*/
/* we force to dissipate the maximum */
alpha = (-J2I);
if (ThisTask==0)
printf("alpha = %g\n",alpha);
}
/*!
*/
void apply_transformation_m1(void)
{
int i,j,k;
/* apply transformation */
if (LocalDissipation)
{
for (k=0;k<N_ngb;k++)
{
i = numlist[k];
for (j=0;j<3;j++)
P[i].Vel[j] += alpha * (P[i].Pos[j] - Xtot[j]);
}
}
else
{
/* compute the (local) total mass M */
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == DISSIPATIVE_TYPE)
{
for (j=0;j<3;j++)
P[i].Vel[j] += alpha * (P[i].Pos[j] - Xtot[j]);
}
}
}
}
/*!
*/
void compute_energy_kin(void)
{
int i,j,k;
Ttot=0;
T=0;
/* apply transformation */
if (LocalDissipation)
{
for (k=0;k<N_ngb;k++)
{
i = numlist[k];
for (j=0;j<3;j++)
T += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
}
}
else
{
/* compute the (local) total mass M */
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == DISSIPATIVE_TYPE)
{
for (j=0;j<3;j++)
T += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
}
}
}
T *= 0.5;
MPI_Allreduce(&T, &Ttot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
}
double get_Total_Gas_energy_kin(void)
{
int i,j;
double totTtot,totT;
totTtot=0;
totT=0;
for(i = 0; i < N_gas; i++)
for (j=0;j<3;j++)
totT += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
totT *= 0.5;
MPI_Allreduce(&totT, &totTtot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
return totTtot;
}
double get_Gas_energy_kin(void)
{
int i,j;
double totT;
totT=0;
for(i = 0; i < N_gas; i++)
for (j=0;j<3;j++)
totT += (P[i].Vel[j]*P[i].Vel[j])*P[i].Mass;
totT *= 0.5;
return totT;
}
double icd_get_alpha(void)
{
return alpha;
}
double icd_get_M(void)
{
return Mtot;
}
double icd_get_X(double x[3])
{
x[0] = Xtot[0];
x[1] = Xtot[1];
x[2] = Xtot[2];
return 0;
}
double icd_get_V(double v[3])
{
v[0] = Vtot[0];
v[1] = Vtot[1];
v[2] = Vtot[2];
return 0;
}
double icd_get_T(void)
{
return Ttot;
}
double icd_get_I(void)
{
return Itot;
}
double icd_get_J(void)
{
return Jtot;
}
double icd_get_DeltaTmin(void)
{
return DeltaTmin;
}
/*!
*/
void icd_find_ngbs(void)
{
int i,j;
MPI_Status status;
if (ThisTask==0)
printf("searching %d ngbs around (x=%g y=%g z=%g)\n",des_nngb,xc[0],xc[1],xc[2]);
/* chose 1 particle */
/* !!! here this is a particles from proc 0 !!! */
//int central_particle_index=100000;
//if (ThisTask==0)
// for (j=0;j<3;j++)
// xc[j] = P[central_particle_index].Pos[j];
/* send to all procs */
MPI_Bcast(&xc, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
/* now, sort particles */
icdd = malloc(sizeof(struct icd_data) * N_gas);
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == DISSIPATIVE_TYPE)
{
icdd[i].index = i;
icdd[i].proc = ThisTask;
icdd[i].key = (P[i].Pos[0]-xc[0])*(P[i].Pos[0]-xc[0])
+ (P[i].Pos[1]-xc[1])*(P[i].Pos[1]-xc[1])
+ (P[i].Pos[2]-xc[2])*(P[i].Pos[2]-xc[2]);
}
}
qsort(icdd, N_gas, sizeof(struct icd_data), compare_icd_key);
//for(i = 0; i < des_nngb; i++)
// {
// printf("(%d) %g (%d) (id=%d) proc=%d\n",ThisTask,icdd[i].key,icdd[i].index,P[icdd[i].index].ID,icdd[i].proc);
// }
/* now, we need to get particles from other procs */
icddGet = malloc(sizeof(struct icd_data) * 2*des_nngb);
if (ThisTask==0)
{
for (j=0;j<des_nngb;j++)
icddGet[j] = icdd[j];
for(j=1;j<NTask;j++)
{
MPI_Recv(&icddGet[des_nngb], des_nngb*sizeof(struct icd_data), MPI_BYTE, j, TAG_N, MPI_COMM_WORLD, &status);
/* get particles to the second part of icddGet */
//printf("\n>>> particles from %d\n\n",j);
//for (i=0;i<2*des_nngb;i++)
// printf(" (%d) %g (%d) (id=%d) proc=%d\n",ThisTask,icddGet[i].key,icddGet[i].index,P[icddGet[i].index].ID,icddGet[i].proc);
/* now, sort */
qsort(icddGet, 2*des_nngb, sizeof(struct icd_data), compare_icd_key);
//printf("\n>>> particles from %d\n\n",j);
//for (i=0;i<2*des_nngb;i++)
// printf(" (%d) %g (%d) (id=%d) proc=%d\n",ThisTask,icddGet[i].key,icddGet[i].index,P[icddGet[i].index].ID,icddGet[i].proc);
}
}
else
{
MPI_Send(icdd, des_nngb*sizeof(struct icd_data), MPI_BYTE, 0, TAG_N, MPI_COMM_WORLD);
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast(icddGet, des_nngb*sizeof(struct icd_data), MPI_BYTE, 0, MPI_COMM_WORLD);
/* at this point, the first des_nngb in icddGet are the closest neighbors */
//for (i=0;i<des_nngb;i++)
// printf("(%d) %g (%d) (id=%d) proc=%d\n",ThisTask,icddGet[i].key,icddGet[i].index,P[icddGet[i].index].ID,icddGet[i].proc);
/* setup a list of local particles belonging to neighbors */
for (i=0;i<des_nngb;i++)
{
if (icddGet[i].proc==ThisTask)
{
numlist[N_ngb]=icddGet[i].index;
N_ngb++;
}
}
free(icdd);
free(icddGet);
/* now, we can compute physical quantities */
//if (N_ngb>0)
// for (i=0;i<N_ngb;i++)
// printf("(%d) (%d) \n",ThisTask,numlist[i]);
}
int compare_icd_key(const void *a, const void *b)
{
if(((struct icd_data *) a)->key < (((struct icd_data *) b)->key))
return -1;
if(((struct icd_data *) a)->key > (((struct icd_data *) b)->key))
return +1;
return 0;
}
#endif // INTEGRAL_CONSERVING_DISSIPATION
c4science · Help