Homec4science
Diffusion pNbody (master) Browse

pNbody/src/PyGadget/src/python_interface.0.c
021ca11fcae0master

Last Change
Raw File

python_interface.0.c
View Options

#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <numpy/libnumarray.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
#define TO_INT(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_INT) ,0) )
#define TO_DOUBLE(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_DOUBLE) ,0) )
#define TO_FLOAT(a) ( (PyArrayObject*) PyArray_CastToType(a, PyArray_DescrFromType(NPY_FLOAT) ,0) )
static int Init()
{
/* main.c */
RestartFlag = 0;
All.CPU_TreeConstruction = All.CPU_TreeWalk = All.CPU_Gravity = All.CPU_Potential = All.CPU_Domain =
All.CPU_Snapshot = All.CPU_Total = All.CPU_CommSum = All.CPU_Imbalance = All.CPU_Hydro =
All.CPU_HydCompWalk = All.CPU_HydCommSumm = All.CPU_HydImbalance =
All.CPU_EnsureNgb = All.CPU_Predict = All.CPU_TimeLine = All.CPU_PM = All.CPU_Peano = 0;
CPUThisRun = 0;
/* from init.c, after read ic */
int i, j;
double a3;
All.Time = All.TimeBegin;
All.Ti_Current = 0;
if(All.ComovingIntegrationOn)
{
All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
a3 = All.Time * All.Time * All.Time;
}
else
{
All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
a3 = 1;
}
set_softenings();
All.NumCurrentTiStep = 0; /* setup some counters */
All.SnapshotFileCount = 0;
if(RestartFlag == 2)
All.SnapshotFileCount = atoi(All.InitCondFile + strlen(All.InitCondFile) - 3) + 1;
All.TotNumOfForces = 0;
All.NumForcesSinceLastDomainDecomp = 0;
if(All.ComovingIntegrationOn)
if(All.PeriodicBoundariesOn == 1)
check_omega();
All.TimeLastStatistics = All.TimeBegin - All.TimeBetStatistics;
if(All.ComovingIntegrationOn) /* change to new velocity variable */
{
for(i = 0; i < NumPart; i++)
for(j = 0; j < 3; j++)
P[i].Vel[j] *= sqrt(All.Time) * All.Time;
}
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
for(j = 0; j < 3; j++)
P[i].GravAccel[j] = 0;
#ifdef PMGRID
for(j = 0; j < 3; j++)
P[i].GravPM[j] = 0;
#endif
P[i].Ti_endstep = 0;
P[i].Ti_begstep = 0;
P[i].OldAcc = 0;
P[i].GravCost = 1;
P[i].Potential = 0;
}
#ifdef PMGRID
All.PM_Ti_endstep = All.PM_Ti_begstep = 0;
#endif
#ifdef FLEXSTEPS
All.PresentMinStep = TIMEBASE;
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
P[i].FlexStepGrp = (int) (TIMEBASE * get_random_number(P[i].ID));
}
#endif
for(i = 0; i < N_gas; i++) /* initialize sph_properties */
{
for(j = 0; j < 3; j++)
{
SphP[i].VelPred[j] = P[i].Vel[j];
SphP[i].HydroAccel[j] = 0;
}
SphP[i].DtEntropy = 0;
if(RestartFlag == 0)
{
SphP[i].Hsml = 0;
SphP[i].Density = -1;
}
}
ngb_treeallocate(MAX_NGB);
force_treeallocate(All.TreeAllocFactor * All.MaxPart, All.MaxPart);
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
Flag_FullStep = 1; /* to ensure that Peano-Hilber order is done */
domain_Decomposition(); /* do initial domain decomposition (gives equal numbers of particles) */
ngb_treebuild(); /* will build tree */
setup_smoothinglengths();
TreeReconstructFlag = 1;
/* at this point, the entropy variable normally contains the
* internal energy, read in from the initial conditions file, unless the file
* explicitly signals that the initial conditions contain the entropy directly.
* Once the density has been computed, we can convert thermal energy to entropy.
*/
#ifndef ISOTHERM_EQS
if(header.flag_entropy_instead_u == 0)
for(i = 0; i < N_gas; i++)
SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
return 1;
}
static void Begrun1()
{
struct global_data_all_processes all;
if(ThisTask == 0)
{
printf("\nThis is Gadget, version `%s'.\n", GADGETVERSION);
printf("\nRunning on %d processors.\n", NTask);
}
read_parameter_file(ParameterFile); /* ... read in parameters for this run */
allocate_commbuffers(); /* ... allocate buffer-memory for particle
exchange during force computation */
set_units();
#if defined(PERIODIC) && (!defined(PMGRID) || defined(FORCETEST))
ewald_init();
#endif
open_outputfiles();
random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
gsl_rng_set(random_generator, 42); /* start-up seed */
#ifdef PMGRID
long_range_init();
#endif
All.TimeLastRestartFile = CPUThisRun;
if(RestartFlag == 0 || RestartFlag == 2)
{
set_random_numbers();
Init(); /* ... read in initial model */
}
else
{
all = All; /* save global variables. (will be read from restart file) */
restart(RestartFlag); /* ... read restart file. Note: This also resets
all variables in the struct `All'.
However, during the run, some variables in the parameter
file are allowed to be changed, if desired. These need to
copied in the way below.
Note: All.PartAllocFactor is treated in restart() separately.
*/
All.MinSizeTimestep = all.MinSizeTimestep;
All.MaxSizeTimestep = all.MaxSizeTimestep;
All.BufferSize = all.BufferSize;
All.BunchSizeForce = all.BunchSizeForce;
All.BunchSizeDensity = all.BunchSizeDensity;
All.BunchSizeHydro = all.BunchSizeHydro;
All.BunchSizeDomain = all.BunchSizeDomain;
All.TimeLimitCPU = all.TimeLimitCPU;
All.ResubmitOn = all.ResubmitOn;
All.TimeBetSnapshot = all.TimeBetSnapshot;
All.TimeBetStatistics = all.TimeBetStatistics;
All.CpuTimeBetRestartFile = all.CpuTimeBetRestartFile;
All.ErrTolIntAccuracy = all.ErrTolIntAccuracy;
All.MaxRMSDisplacementFac = all.MaxRMSDisplacementFac;
All.ErrTolForceAcc = all.ErrTolForceAcc;
All.TypeOfTimestepCriterion = all.TypeOfTimestepCriterion;
All.TypeOfOpeningCriterion = all.TypeOfOpeningCriterion;
All.NumFilesWrittenInParallel = all.NumFilesWrittenInParallel;
All.TreeDomainUpdateFrequency = all.TreeDomainUpdateFrequency;
All.SnapFormat = all.SnapFormat;
All.NumFilesPerSnapshot = all.NumFilesPerSnapshot;
All.MaxNumNgbDeviation = all.MaxNumNgbDeviation;
All.ArtBulkViscConst = all.ArtBulkViscConst;
All.OutputListOn = all.OutputListOn;
All.CourantFac = all.CourantFac;
All.OutputListLength = all.OutputListLength;
memcpy(All.OutputListTimes, all.OutputListTimes, sizeof(double) * All.OutputListLength);
strcpy(All.ResubmitCommand, all.ResubmitCommand);
strcpy(All.OutputListFilename, all.OutputListFilename);
strcpy(All.OutputDir, all.OutputDir);
strcpy(All.RestartFile, all.RestartFile);
strcpy(All.EnergyFile, all.EnergyFile);
strcpy(All.InfoFile, all.InfoFile);
strcpy(All.CpuFile, all.CpuFile);
strcpy(All.TimingsFile, all.TimingsFile);
strcpy(All.SnapshotFileBase, all.SnapshotFileBase);
if(All.TimeMax != all.TimeMax)
readjust_timebase(All.TimeMax, all.TimeMax);
}
}
static void Begrun2()
{
#ifdef PMGRID
long_range_init_regionsize();
#endif
if(All.ComovingIntegrationOn)
init_drift_table();
if(RestartFlag == 2)
All.Ti_nextoutput = find_next_outputtime(All.Ti_Current + 1);
else
All.Ti_nextoutput = find_next_outputtime(All.Ti_Current);
All.TimeLastRestartFile = CPUThisRun;
}
/************************************************************/
/* PYTHON INTERFACE */
/************************************************************/
static PyObject *gadget_Info(PyObject *self, PyObject *args, PyObject *kwds)
{
printf("I am proc %d among %d procs.\n",ThisTask,NTask);
return Py_BuildValue("i",1);
}
static PyObject *gadget_InitMPI(PyObject *self, PyObject *args, PyObject *kwds)
{
//MPI_Init(0, 0); /* this is done in mpi4py */
MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
MPI_Comm_size(MPI_COMM_WORLD, &NTask);
for(PTask = 0; NTask > (1 << PTask); PTask++);
return Py_BuildValue("i",1);
}
static PyObject * gadget_InitDefaultParameters(PyObject* self)
{
/* list of Gadget parameters */
All.ComovingIntegrationOn = 0;
All.PeriodicBoundariesOn = 0;
All.Omega0 = 0;
All.OmegaLambda = 0;
All.OmegaBaryon = 0;
All.HubbleParam = 0;
All.BoxSize = 0;
All.ErrTolTheta = 0.7;
All.TypeOfOpeningCriterion = 0;
All.ErrTolForceAcc = 0.005;
All.DesNumNgb = 33;
All.MaxNumNgbDeviation = 3;
All.PartAllocFactor = 2.0;
All.TreeAllocFactor = 2.0;
All.BufferSize = 30;
All.MinGasHsmlFractional = 0.25;
All.SofteningGas = 0.5;
All.SofteningHalo = 0.5;
All.SofteningDisk = 0.5;
All.SofteningBulge = 0.5;
All.SofteningStars = 0.5;
All.SofteningBndry = 0.5;
All.SofteningGasMaxPhys = 0.5;
All.SofteningHaloMaxPhys = 0.5;
All.SofteningDiskMaxPhys = 0.5;
All.SofteningBulgeMaxPhys = 0.5;
All.SofteningStarsMaxPhys = 0.5;
All.SofteningBndryMaxPhys = 0.5;
return Py_BuildValue("i",1);
}
static PyObject *gadget_LoadParticles(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *ntype,*pos,*vel,*mass,*num,*tpe;
static char *kwlist[] = {"npart", "pos","vel","mass","num","tpe", NULL};
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OOOOOO",kwlist,&ntype,&pos,&vel,&mass,&num,&tpe))
return Py_BuildValue("i",1);
/* check type */
if (!(PyArray_Check(pos)))
{
PyErr_SetString(PyExc_ValueError,"aruments 1 must be array.");
return NULL;
}
/* check type */
if (!(PyArray_Check(mass)))
{
PyErr_SetString(PyExc_ValueError,"aruments 2 must be array.");
return NULL;
}
/* check dimension */
if ( (pos->nd!=2))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 1 must be 2.");
return NULL;
}
/* check dimension */
if ( (mass->nd!=1))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 2 must be 1.");
return NULL;
}
/* check size */
if ( (pos->dimensions[1]!=3))
{
PyErr_SetString(PyExc_ValueError,"First size of argument must be 3.");
return NULL;
}
/* check size */
if ( (pos->dimensions[0]!=mass->dimensions[0]))
{
PyErr_SetString(PyExc_ValueError,"Size of argument 1 must be similar to argument 2.");
return NULL;
}
/* ensure double */
ntype = TO_INT(ntype);
pos = TO_FLOAT(pos);
vel = TO_FLOAT(vel);
mass = TO_FLOAT(mass);
num = TO_FLOAT(num);
tpe = TO_FLOAT(tpe);
/* count number of particles */
NtypeLocal[0] = *(int*) (ntype->data + 0*(ntype->strides[0]));
NtypeLocal[1] = *(int*) (ntype->data + 1*(ntype->strides[0]));
NtypeLocal[2] = *(int*) (ntype->data + 2*(ntype->strides[0]));
NtypeLocal[3] = *(int*) (ntype->data + 3*(ntype->strides[0]));
NtypeLocal[4] = *(int*) (ntype->data + 4*(ntype->strides[0]));
NtypeLocal[5] = *(int*) (ntype->data + 5*(ntype->strides[0]));
NumPart = 0;
N_gas = NtypeLocal[0];
for (i = 0; i < 6; i++)
NumPart += NtypeLocal[i];
MPI_Allreduce(&NumPart, &All.TotNumPart, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&N_gas, &All.TotN_gas, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask);
All.MaxPartSph = All.PartAllocFactor * (All.TotN_gas / NTask);
All.MinGasHsml = All.MinGasHsmlFractional * All.ForceSoftening[0];
All.BunchSizeDomain =
(All.BufferSize * 1024 * 1024) / (sizeof(struct particle_data) + sizeof(struct sph_particle_data) +
sizeof(peanokey));
if(All.BunchSizeDomain & 1)
All.BunchSizeDomain -= 1; /* make sure that All.BunchSizeDomain is even
--> 8-byte alignment of DomainKeyBuf for 64bit processors */
All.BunchSizeForce =
(All.BufferSize * 1024 * 1024) / (sizeof(struct gravdata_index) + 2 * sizeof(struct gravdata_in));
if(All.BunchSizeForce & 1)
All.BunchSizeForce -= 1; /* make sure that All.BunchSizeForce is an even number
--> 8-byte alignment for 64bit processors */
All.BunchSizeDensity =
(All.BufferSize * 1024 * 1024) / (2 * sizeof(struct densdata_in) + 2 * sizeof(struct densdata_out));
/*********************/
/* some allocation */
/*********************/
if(!(CommBuffer = malloc(bytes = All.BufferSize * 1024 * 1024)))
{
printf("failed to allocate memory for `CommBuffer' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(2);
}
Exportflag = malloc(NTask * sizeof(char));
GravDataIndexTable = (struct gravdata_index *) CommBuffer;
GravDataIn = (struct gravdata_in *) (GravDataIndexTable + All.BunchSizeForce);
GravDataGet = GravDataIn + All.BunchSizeForce;
GravDataOut = GravDataIn;
GravDataResult = GravDataGet;
DensDataIn = (struct densdata_in *) CommBuffer;
DensDataGet = DensDataIn + All.BunchSizeDensity;
DensDataResult = (struct densdata_out *) (DensDataGet + All.BunchSizeDensity);
DensDataPartialResult = DensDataResult + All.BunchSizeDensity;
/*********************/
/* create P */
/*********************/
if(!(P = malloc(bytes = All.MaxPart * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `P' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphP = malloc(bytes = All.MaxPartSph * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphP' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
endrun(1);
}
/*********************/
/* init P */
/*********************/
for (i = 0; i < pos->dimensions[0]; i++)
{
P[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
P[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
P[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
P[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
P[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
P[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
P[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
P[i].ID = *(unsigned int *) (num->data + i*(num->strides[0]));
P[i].Type = *(int *) (tpe->data + i*(tpe->strides[0])); /* this should be changed... */
//P[i].Active = 1;
}
/***************************************
* some inits *
/***************************************/
RestartFlag = 0;
Begrun();
/***************************************
* init ewald *
/***************************************/
//if (All.PeriodicBoundariesOn)
// ewald_init();
/***************************************
* domain decomposition construction *
/***************************************/
All.NumForcesSinceLastDomainDecomp = 1; /* a changer !!!! */
All.TreeDomainUpdateFrequency = 0; /* a changer !!!! */
allocate_commbuffers();
domain_Decomposition();
return Py_BuildValue("i",1);
}
/* definition of the method table */
static PyMethodDef gadgetMethods[] = {
{"Info", gadget_Info, METH_VARARGS,
"give some info"},
{"InitMPI", gadget_InitMPI, METH_VARARGS,
"Init MPI"},
{"InitDefaultParameters", gadget_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"LoadParticles", gadget_LoadParticles, METH_VARARGS,
"LoadParticles partilces"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initgadget(void)
{
(void) Py_InitModule("gadget", gadgetMethods);
import_array();
}
c4science · Help