Homec4science
Diffusion Gear (master) Browse

Gear/src/python_interface.c
710addfceec3master

Last Change
Raw File

python_interface.c
View Options

#ifdef PY_INTERFACE
#include <Python.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <numpy/arrayobject.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 pyGadget, 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();
}
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()
{
if(RestartFlag == 0 || RestartFlag == 2)
Init(); /* ... read in initial model */
#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.InitCondFile ="ICs/cluster_littleendian.dat";
All.OutputDir ="cluster/";
All.EnergyFile ="energy.txt";
All.InfoFile ="info.txt";
All.TimingsFile ="timings.txt";
All.CpuFile ="cpu.txt";
All.RestartFile ="restart";
All.SnapshotFileBase ="snapshot";
All.OutputListFilename ="parameterfiles/outputs_lcdm_gas.txt";
*/
/* CPU time -limit */
All.TimeLimitCPU = 36000; /* = 10 hours */
All.ResubmitOn = 0;
//All.ResubmitCommand = "my-scriptfile";
All.ICFormat = 1;
All.SnapFormat = 1;
All.ComovingIntegrationOn = 0;
All.TypeOfTimestepCriterion = 0;
All.OutputListOn = 0;
All.PeriodicBoundariesOn = 0;
/* Caracteristics of run */
All.TimeBegin = 0.0; /*% Begin of the simulation (z=23)*/
All.TimeMax = 1.0;
All.Omega0 = 0;
All.OmegaLambda = 0;
All.OmegaBaryon = 0;
All.HubbleParam = 0;
All.BoxSize = 0;
/* Output frequency */
All.TimeBetSnapshot = 0.1;
All.TimeOfFirstSnapshot = 0.0; /*% 5 constant steps in log(a) */
All.CpuTimeBetRestartFile = 36000.0; /* here in seconds */
All.TimeBetStatistics = 0.05;
All.NumFilesPerSnapshot = 1;
All.NumFilesWrittenInParallel = 1;
/* Accuracy of time integration */
All.ErrTolIntAccuracy = 0.025;
All.MaxRMSDisplacementFac = 0.2;
All.CourantFac = 0.15;
All.MaxSizeTimestep = 0.03;
All.MinSizeTimestep = 0.0;
/* Tree algorithm, force accuracy, domain update frequency */
All.ErrTolTheta = 0.7;
All.TypeOfOpeningCriterion = 0;
All.ErrTolForceAcc = 0.005;
All.TreeDomainUpdateFrequency = 0.1;
/* Further parameters of SPH */
All.DesNumNgb = 50;
All.MaxNumNgbDeviation = 2;
All.ArtBulkViscConst = 0.8;
All.InitGasTemp = 0;
All.MinGasTemp = 0;
/* Memory allocation */
All.PartAllocFactor = 2.0;
All.TreeAllocFactor = 2.0;
All.BufferSize = 30;
/* System of units */
All.UnitLength_in_cm = 3.085678e21; /* 1.0 kpc */
All.UnitMass_in_g = 1.989e43; /* 1.0e10 solar masses */
All.UnitVelocity_in_cm_per_s = 1e5; /* 1 km/sec */
All.GravityConstantInternal = 0;
/* Softening lengths */
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_GetParameters()
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
/* CPU time -limit */
key = PyString_FromString("TimeLimitCPU");
value = PyFloat_FromDouble(All.TimeLimitCPU);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ResubmitOn");
value = PyFloat_FromDouble(All.ResubmitOn);
PyDict_SetItem(dict,key,value);
//All.ResubmitCommand
key = PyString_FromString("ICFormat");
value = PyInt_FromLong(All.ICFormat);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SnapFormat");
value = PyInt_FromLong(All.SnapFormat);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ComovingIntegrationOn");
value = PyInt_FromLong(All.ComovingIntegrationOn);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TypeOfTimestepCriterion");
value = PyInt_FromLong(All.TypeOfTimestepCriterion);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OutputListOn");
value = PyInt_FromLong(All.OutputListOn);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("PeriodicBoundariesOn");
value = PyInt_FromLong(All.PeriodicBoundariesOn);
PyDict_SetItem(dict,key,value);
/* Caracteristics of run */
key = PyString_FromString("TimeBegin");
value = PyFloat_FromDouble(All.TimeBegin);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeMax");
value = PyFloat_FromDouble(All.TimeMax);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("Omega0");
value = PyFloat_FromDouble(All.Omega0);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaLambda");
value = PyFloat_FromDouble(All.OmegaLambda);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("OmegaBaryon");
value = PyFloat_FromDouble(All.OmegaBaryon);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("HubbleParam");
value = PyFloat_FromDouble(All.HubbleParam);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("BoxSize");
value = PyFloat_FromDouble(All.BoxSize);
PyDict_SetItem(dict,key,value);
/* Output frequency */
key = PyString_FromString("TimeBetSnapshot");
value = PyFloat_FromDouble(All.TimeBetSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeOfFirstSnapshot");
value = PyFloat_FromDouble(All.TimeOfFirstSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("CpuTimeBetRestartFile");
value = PyFloat_FromDouble(All.CpuTimeBetRestartFile);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TimeBetStatistics");
value = PyFloat_FromDouble(All.TimeBetStatistics);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("NumFilesPerSnapshot");
value = PyInt_FromLong(All.NumFilesPerSnapshot);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("NumFilesWrittenInParallel");
value = PyInt_FromLong(All.NumFilesWrittenInParallel);
PyDict_SetItem(dict,key,value);
/* Accuracy of time integration */
key = PyString_FromString("ErrTolIntAccuracy");
value = PyFloat_FromDouble(All.ErrTolIntAccuracy);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxRMSDisplacementFac");
value = PyFloat_FromDouble(All.MaxRMSDisplacementFac);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("CourantFac");
value = PyFloat_FromDouble(All.CourantFac);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxSizeTimestep");
value = PyFloat_FromDouble(All.MaxSizeTimestep);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MinSizeTimestep");
value = PyFloat_FromDouble(All.MinSizeTimestep);
PyDict_SetItem(dict,key,value);
/* Tree algorithm, force accuracy, domain update frequency */
key = PyString_FromString("ErrTolTheta");
value = PyFloat_FromDouble(All.ErrTolTheta);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TypeOfOpeningCriterion");
value = PyInt_FromLong(All.TypeOfOpeningCriterion);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ErrTolForceAcc");
value = PyFloat_FromDouble(All.ErrTolForceAcc);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TreeDomainUpdateFrequency");
value = PyFloat_FromDouble(All.TreeDomainUpdateFrequency);
PyDict_SetItem(dict,key,value);
/* Further parameters of SPH */
key = PyString_FromString("DesNumNgb");
value = PyInt_FromLong(All.DesNumNgb);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MaxNumNgbDeviation");
value = PyInt_FromLong(All.MaxNumNgbDeviation);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("ArtBulkViscConst");
value = PyInt_FromLong(All.ArtBulkViscConst);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("InitGasTemp");
value = PyInt_FromLong(All.InitGasTemp);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("MinGasTemp");
value = PyInt_FromLong(All.MinGasTemp);
PyDict_SetItem(dict,key,value);
/* Memory allocation */
key = PyString_FromString("PartAllocFactor");
value = PyFloat_FromDouble(All.PartAllocFactor);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("TreeAllocFactor");
value = PyFloat_FromDouble(All.TreeAllocFactor);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("BufferSize");
value = PyInt_FromLong(All.BufferSize);
PyDict_SetItem(dict,key,value);
/* System of units */
key = PyString_FromString("UnitLength_in_cm");
value = PyFloat_FromDouble(All.UnitLength_in_cm);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitMass_in_g");
value = PyFloat_FromDouble(All.UnitMass_in_g);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("UnitVelocity_in_cm_per_s");
value = PyFloat_FromDouble(All.UnitVelocity_in_cm_per_s);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("GravityConstantInternal");
value = PyFloat_FromDouble(All.GravityConstantInternal);
PyDict_SetItem(dict,key,value);
/* Softening lengths */
key = PyString_FromString("MinGasHsmlFractional");
value = PyFloat_FromDouble(All.MinGasHsmlFractional);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningGas");
value = PyFloat_FromDouble(All.SofteningGas);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningHalo");
value = PyFloat_FromDouble(All.SofteningHalo);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningDisk");
value = PyFloat_FromDouble(All.SofteningDisk);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBulge");
value = PyFloat_FromDouble(All.SofteningBulge);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningStars");
value = PyFloat_FromDouble(All.SofteningStars);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBndry");
value = PyFloat_FromDouble(All.SofteningBndry);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningGasMaxPhys");
value = PyFloat_FromDouble(All.SofteningGasMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningHaloMaxPhys");
value = PyFloat_FromDouble(All.SofteningHaloMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningDiskMaxPhys");
value = PyFloat_FromDouble(All.SofteningDiskMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBulgeMaxPhys");
value = PyFloat_FromDouble(All.SofteningBulgeMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningStarsMaxPhys");
value = PyFloat_FromDouble(All.SofteningStarsMaxPhys);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("SofteningBndryMaxPhys");
value = PyFloat_FromDouble(All.SofteningBndryMaxPhys);
PyDict_SetItem(dict,key,value);
/*
key = PyString_FromString("OutputInfo");
value = PyFloat_FromDouble(All.OutputInfo);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("PeanoHilbertOrder");
value = PyFloat_FromDouble(All.PeanoHilbertOrder);
PyDict_SetItem(dict,key,value);
*/
return Py_BuildValue("O",dict);
}
static PyObject * gadget_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
/* CPU time -limit */
if(strcmp(PyString_AsString(key), "TimeLimitCPU")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeLimitCPU = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "ResubmitOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ResubmitOn = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "ICFormat")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ICFormat = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "SnapFormat")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SnapFormat = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ComovingIntegrationOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ComovingIntegrationOn = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "TypeOfTimestepCriterion")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TypeOfTimestepCriterion = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "OutputListOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OutputListOn = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "PeriodicBoundariesOn")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.PeriodicBoundariesOn = PyInt_AsLong(value);
}
/* Caracteristics of run */
if(strcmp(PyString_AsString(key), "TimeBegin")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBegin = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeMax")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeMax = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "Omega0")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.Omega0 = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaLambda")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaLambda = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "OmegaBaryon")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.OmegaBaryon = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "HubbleParam")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.HubbleParam = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "BoxSize")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.BoxSize = PyFloat_AsDouble(value);
}
/* Output frequency */
if(strcmp(PyString_AsString(key), "TimeBetSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBetSnapshot = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeOfFirstSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeOfFirstSnapshot = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "CpuTimeBetRestartFile")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.CpuTimeBetRestartFile = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TimeBetStatistics")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TimeBetStatistics = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "NumFilesPerSnapshot")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.NumFilesPerSnapshot = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "NumFilesWrittenInParallel")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.NumFilesWrittenInParallel = PyInt_AsLong(value);
}
/* Accuracy of time integration */
if(strcmp(PyString_AsString(key), "ErrTolIntAccuracy")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolIntAccuracy = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MaxRMSDisplacementFac")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxRMSDisplacementFac = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "CourantFac")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.CourantFac = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MaxSizeTimestep")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxSizeTimestep = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "MinSizeTimestep")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinSizeTimestep = PyFloat_AsDouble(value);
}
/* Tree algorithm, force accuracy, domain update frequency */
if(strcmp(PyString_AsString(key), "ErrTolTheta")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolTheta = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TypeOfOpeningCriterion")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TypeOfOpeningCriterion = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ErrTolForceAcc")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ErrTolForceAcc = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TreeDomainUpdateFrequency")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TreeDomainUpdateFrequency = PyFloat_AsDouble(value);
}
/* Further parameters of SPH */
if(strcmp(PyString_AsString(key), "DesNumNgb")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.DesNumNgb = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "MaxNumNgbDeviation")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MaxNumNgbDeviation = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "ArtBulkViscConst")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.ArtBulkViscConst = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "InitGasTemp")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.InitGasTemp = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "MinGasTemp")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinGasTemp = PyInt_AsLong(value);
}
/* Memory allocation */
if(strcmp(PyString_AsString(key), "PartAllocFactor")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.PartAllocFactor = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "TreeAllocFactor")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.TreeAllocFactor = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "BufferSize")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.BufferSize = PyInt_AsLong(value);
}
/* System of units */
if(strcmp(PyString_AsString(key), "UnitLength_in_cm")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitLength_in_cm = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitMass_in_g")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitMass_in_g = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "UnitVelocity_in_cm_per_s")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.UnitVelocity_in_cm_per_s = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "GravityConstantInternal")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.GravityConstantInternal = PyFloat_AsDouble(value);
}
/* Softening lengths */
if(strcmp(PyString_AsString(key), "MinGasHsmlFractional")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.MinGasHsmlFractional = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningGas")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningGas = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningHalo")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningHalo = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningDisk")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningDisk = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBulge")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBulge = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningStars")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningStars = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBndry")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBndry = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningGasMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningGasMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningHaloMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningHaloMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningDiskMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningDiskMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBulgeMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBulgeMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningStarsMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningStarsMaxPhys = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "SofteningBndryMaxPhys")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.SofteningBndryMaxPhys = PyFloat_AsDouble(value);
}
}
}
return Py_BuildValue("i",1);
}
static PyObject *
gadget_check_parser(PyObject *self, PyObject *args, PyObject *keywds)
{
int voltage;
char *state = "a stiff";
char *action = "voom";
char *type = "Norwegian Blue";
static char *kwlist[] = {"voltage", "state", "action", "type", NULL};
if (!PyArg_ParseTupleAndKeywords(args, keywds, "i|sss", kwlist,
&voltage, &state, &action, &type))
return NULL;
printf("-- This parrot wouldn't %s if you put %i Volts through it.\n",
action, voltage);
printf("-- Lovely plumage, the %s -- It's %s!\n", type, state);
Py_INCREF(Py_None);
return Py_None;
}
static PyObject *gadget_Free(PyObject *self, PyObject *args, PyObject *kwds)
{
free_memory();
ngb_treefree();
force_treefree();
free(Exportflag );
free(DomainStartList);
free(DomainEndList);
free(TopNodes);
free(DomainWork);
free(DomainCount);
free(DomainCountSph);
free(DomainTask);
free(DomainNodeIndex);
free(DomainTreeNodeLen);
free(DomainHmax);
free(DomainMoment);
free(CommBuffer);
gsl_rng_free(random_generator);
return Py_BuildValue("i",1);
}
static PyObject *gadget_LoadParticles(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *ntype=Py_None;
PyArrayObject *pos=Py_None;
PyArrayObject *vel=Py_None;
PyArrayObject *mass=Py_None;
PyArrayObject *num=Py_None;
PyArrayObject *tpe=Py_None;
PyArrayObject *u=Py_None;
PyArrayObject *rho=Py_None;
static char *kwlist[] = {"npart", "pos","vel","mass","num","tpe","u","rho", NULL};
if (! PyArg_ParseTupleAndKeywords(args,kwds, "OOOOOO|OO",kwlist,&ntype,&pos,&vel,&mass,&num,&tpe,&u,&rho ))
return NULL;
/* check type */
if (!(PyArray_Check(pos)))
{
PyErr_SetString(PyExc_ValueError,"aruments 2 must be array.");
return NULL;
}
/* check type */
if (!(PyArray_Check(mass)))
{
PyErr_SetString(PyExc_ValueError,"aruments 3 must be array.");
return NULL;
}
/* check dimension */
if ( (pos->nd!=2))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 2 must be 2.");
return NULL;
}
/* check dimension */
if ( (mass->nd!=1))
{
PyErr_SetString(PyExc_ValueError,"Dimension of argument 3 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 2 must be similar to argument 3.");
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);
/* optional variables */
if (u!=Py_None)
u = TO_FLOAT(u);
if (rho!=Py_None)
rho = TO_FLOAT(rho);
/***************************************
* some inits *
/***************************************/
RestartFlag = 0;
Begrun1();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPart = 0;
N_gas = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPart += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPart!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"Numpart != pos->dimensions[0].");
return NULL;
}
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);
/*********************/
/* allocate 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 < NumPart; 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]));
//P[i].Active = 1;
}
/*********************/
/* init SphP */
/*********************/
if (u!=Py_None)
for (i = 0; i < NumPart; i++)
{
SphP[i].Entropy = *(float *) (u->data + i*(u->strides[0]));
//#ifndef ISOTHERM_EQS
// SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
//#endif
}
if (rho!=Py_None)
for (i = 0; i < NumPart; i++)
{
SphP[i].Density = *(float *) (rho->data + i*(rho->strides[0]));
}
/***************************************
* END LOAD PARTICLES *
/***************************************/
/***************************************
* finish inits *
/***************************************/
Begrun2();
/***************************************
* free memory *
/***************************************/
/* free the memory allocated,
because these vectors where casted
and their memory is now handeled by the C part */
Py_DECREF(ntype);
Py_DECREF(pos);
Py_DECREF(vel);
Py_DECREF(mass);
Py_DECREF(num);
Py_DECREF(tpe);
/* optional variables */
if (u!=Py_None)
Py_DECREF(u);
if (rho!=Py_None)
Py_DECREF(rho);
return Py_BuildValue("i",1);
}
static PyObject *gadget_LoadParticlesQ(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);
/***************************************
* some inits *
/***************************************/
allocate_commbuffersQ();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPartQ = 0;
N_gasQ = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPartQ += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPartQ!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"NumpartQ != pos->dimensions[0].");
return NULL;
}
MPI_Allreduce(&NumPartQ, &All.TotNumPartQ, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&N_gasQ, &All.TotN_gasQ, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
All.MaxPartQ = All.PartAllocFactor * (All.TotNumPartQ / NTask);
All.MaxPartSphQ = All.PartAllocFactor * (All.TotN_gasQ / NTask);
/*********************/
/* allocate Q */
/*********************/
if(!(Q = malloc(bytes = All.MaxPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = All.MaxPartSphQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
endrun(1);
}
/*********************/
/* init P */
/*********************/
for (i = 0; i < NumPartQ; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
Q[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
Q[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
Q[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
Q[i].ID = *(unsigned int *) (num->data + i*(num->strides[0]));
Q[i].Type = *(int *) (tpe->data + i*(tpe->strides[0]));
//Q[i].Active = 1;
}
/***************************************
* END LOAD PARTILES *
/***************************************/
domain_DecompositionQ();
/***************************************
* finish inits *
/***************************************/
return Py_BuildValue("i",1);
}
static PyObject *gadget_AllPotential(PyObject *self)
{
compute_potential();
return Py_BuildValue("i",1);
}
static PyObject * gadget_GetAllPotential(PyObject* self)
{
PyArrayObject *pot;
npy_intp ld[1];
int i;
ld[0] = NumPart;
pot = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < pot->dimensions[0]; i++)
{
*(float *) (pot->data + i*(pot->strides[0])) = P[i].Potential;
}
return PyArray_Return(pot);
}
static PyObject *gadget_AllAcceleration(PyObject *self)
{
NumForceUpdate = NumPart;
gravity_tree();
return Py_BuildValue("i",1);
}
static PyObject * gadget_GetAllAcceleration(PyObject* self)
{
PyArrayObject *acc;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
acc = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < acc->dimensions[0]; i++)
{
*(float *) (acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = P[i].GravAccel[0];
*(float *) (acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = P[i].GravAccel[1];
*(float *) (acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = P[i].GravAccel[2];
}
return PyArray_Return(acc);
}
static PyObject *gadget_GetAllDensities(PyObject* self)
{
PyArrayObject *rho;
npy_intp ld[1];
int i;
ld[0] = N_gas;
rho = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < rho->dimensions[0]; i++)
{
*(float *) (rho->data + i*(rho->strides[0])) = SphP[i].Density;
}
return PyArray_Return(rho);
}
static PyObject *gadget_GetAllHsml(PyObject* self)
{
PyArrayObject *hsml;
npy_intp ld[1];
int i;
ld[0] = N_gas;
hsml = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < hsml->dimensions[0]; i++)
{
*(float *) (hsml->data + i*(hsml->strides[0])) = SphP[i].Hsml;
}
return PyArray_Return(hsml);
}
static PyObject *gadget_GetAllPositions(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = P[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[i].Pos[2];
}
return PyArray_Return(pos);
}
static PyObject *gadget_GetAllVelocities(PyObject* self)
{
PyArrayObject *vel;
npy_intp ld[2];
int i;
ld[0] = NumPart;
ld[1] = 3;
vel = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < vel->dimensions[0]; i++)
{
*(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]) = P[i].Vel[0];
*(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]) = P[i].Vel[1];
*(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]) = P[i].Vel[2];
}
return PyArray_Return(vel);
}
static PyObject *gadget_GetAllMasses(PyObject* self)
{
PyArrayObject *mass;
npy_intp ld[1];
int i;
ld[0] = NumPart;
mass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < mass->dimensions[0]; i++)
{
*(float *) (mass->data + i*(mass->strides[0])) = P[i].Mass;
}
return PyArray_Return(mass);
}
static PyObject *gadget_GetAllEntropy(PyObject* self)
{
PyArrayObject *entropy;
npy_intp ld[1];
int i;
ld[0] = NumPart;
entropy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < entropy->dimensions[0]; i++)
{
*(float *) (entropy->data + i*(entropy->strides[0])) = SphP[i].Entropy;
}
return PyArray_Return(entropy);
}
static PyObject *gadget_GetAllEnergySpec(PyObject* self)
{
PyArrayObject *energy;
npy_intp ld[1];
int i;
double a3inv;
ld[0] = NumPart;
energy = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
if(All.ComovingIntegrationOn)
{
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
a3inv = 1;
for (i = 0; i < energy->dimensions[0]; i++)
{
#ifdef ISOTHERM_EQS
*(float *) (energy->data + i*(energy->strides[0])) = SphP[i].Entropy;
#else
a3inv = 1.;
*(float *) (energy->data + i*(energy->strides[0])) = dmax(All.MinEgySpec,SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].Density * a3inv, GAMMA_MINUS1));
#endif
}
return PyArray_Return(energy);
}
static PyObject *gadget_GetAllID(PyObject* self)
{
PyArrayObject *id;
npy_intp ld[1];
int i;
ld[0] = NumPart;
id = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < id->dimensions[0]; i++)
{
*(float *) (id->data + i*(id->strides[0])) = P[i].ID;
}
return PyArray_Return(id);
}
static PyObject *gadget_GetAllTypes(PyObject* self)
{
PyArrayObject *type;
npy_intp ld[1];
int i;
ld[0] = NumPart;
type = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < type->dimensions[0]; i++)
{
*(int *) (type->data + i*(type->strides[0])) = P[i].Type;
}
return PyArray_Return(type);
}
static PyObject *gadget_GetAllPositionsQ(PyObject* self)
{
PyArrayObject *pos;
npy_intp ld[2];
int i;
ld[0] = NumPartQ;
ld[1] = 3;
pos = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = Q[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = Q[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = Q[i].Pos[2];
}
return PyArray_Return(pos);
}
static PyObject *gadget_GetAllVelocitiesQ(PyObject* self)
{
PyArrayObject *vel;
npy_intp ld[2];
int i;
ld[0] = NumPartQ;
ld[1] = 3;
vel = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < vel->dimensions[0]; i++)
{
*(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]) = Q[i].Vel[0];
*(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]) = Q[i].Vel[1];
*(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]) = Q[i].Vel[2];
}
return PyArray_Return(vel);
}
static PyObject *gadget_GetAllMassesQ(PyObject* self)
{
PyArrayObject *mass;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
mass = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
for (i = 0; i < mass->dimensions[0]; i++)
{
*(float *) (mass->data + i*(mass->strides[0])) = Q[i].Mass;
}
return PyArray_Return(mass);
}
static PyObject *gadget_GetAllIDQ(PyObject* self)
{
PyArrayObject *id;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
id = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < id->dimensions[0]; i++)
{
*(float *) (id->data + i*(id->strides[0])) = Q[i].ID;
}
return PyArray_Return(id);
}
static PyObject *gadget_GetAllTypesQ(PyObject* self)
{
PyArrayObject *type;
npy_intp ld[1];
int i;
ld[0] = NumPartQ;
type = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_INT);
for (i = 0; i < type->dimensions[0]; i++)
{
*(int *) (type->data + i*(type->strides[0])) = Q[i].Type;
}
return PyArray_Return(type);
}
static PyObject *gadget_GetPos(PyObject *self, PyObject *args, PyObject *kwds)
{
int i,j;
size_t bytes;
PyArrayObject *pos;
if (! PyArg_ParseTuple(args, "O",&pos))
return PyString_FromString("error : GetPos");
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]) = P[i].Pos[0];
*(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]) = P[i].Pos[1];
*(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]) = P[i].Pos[2];
}
//return PyArray_Return(Py_None);
return Py_BuildValue("i",1);
}
static PyObject * gadget_Potential(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *pot;
int i;
npy_intp ld[1];
int input_dimension;
size_t bytes;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
pos = TO_FLOAT(pos);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
pot = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
NumPartQ = pos->dimensions[0];
All.ForceSofteningQ = eps;
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Type = 0;
Q[i].Mass = 0;
Q[i].Potential = 0;
}
compute_potential_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(pot->data + i*(pot->strides[0])) = Q[i].Potential;
}
free(Q);
free(SphQ);
return PyArray_Return(pot);
}
static PyObject * gadget_Acceleration(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *acc;
int i;
int input_dimension;
size_t bytes;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
pos = TO_FLOAT(pos);
/* create a NumPy object */
acc = (PyArrayObject *) PyArray_SimpleNew(pos->nd,pos->dimensions,PyArray_FLOAT);
NumPartQ = pos->dimensions[0];
All.ForceSofteningQ = eps;
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
Q[i].Type = 0;
Q[i].Mass = 0;
Q[i].GravAccel[0] = 0;
Q[i].GravAccel[1] = 0;
Q[i].GravAccel[2] = 0;
}
gravity_tree_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = Q[i].GravAccel[0];
*(float *)(acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = Q[i].GravAccel[1];
*(float *)(acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = Q[i].GravAccel[2];
}
free(Q);
free(SphQ);
return PyArray_Return(acc);
}
static PyObject * gadget_InitHsml(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml;
if (! PyArg_ParseTuple(args, "OO",&pos,&hsml))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vden,*vhsml;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vden = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
vhsml = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
setup_smoothinglengths_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vhsml->data + i*(vhsml->strides[0])) = SphQ[i].Hsml;
*(float *)(vden->data + i*(vden->strides[0])) = SphQ[i].Density;
}
free(Q);
free(SphQ);
return Py_BuildValue("OO",vden,vhsml);
}
static PyObject * gadget_Density(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml;
if (! PyArg_ParseTuple(args, "OO",&pos,&hsml))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vden,*vhsml;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vden = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
vhsml = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
density_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vhsml->data + i*(vhsml->strides[0])) = SphQ[i].Hsml;
*(float *)(vden->data + i*(vden->strides[0])) = SphQ[i].Density;
}
free(Q);
free(SphQ);
return Py_BuildValue("OO",vden,vhsml);
}
static PyObject * gadget_SphEvaluate(PyObject* self, PyObject *args)
{
PyArrayObject *pos,*hsml,*obs;
if (! PyArg_ParseTuple(args, "OOO",&pos,&hsml,&obs))
return PyString_FromString("error");
int i;
int input_dimension;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->dimensions[0] != hsml->dimensions[0])
PyErr_SetString(PyExc_ValueError,"pos and hsml must have the same dimension.");
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
pos = TO_FLOAT(pos);
hsml = TO_FLOAT(hsml);
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=pos->dimensions[0];
vobs = (PyArrayObject *) PyArray_SimpleNew(1,pos->dimensions,pos->descr->type_num);
NumPartQ = pos->dimensions[0];
N_gasQ = NumPartQ;
All.Ti_Current=1; /* need to flag active particles */
if(!(Q = malloc(bytes = NumPartQ * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `Q' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
if(!(SphQ = malloc(bytes = NumPartQ * sizeof(struct sph_particle_data))))
{
printf("failed to allocate memory for `SphQ' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
Q[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
Q[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
Q[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
SphQ[i].Hsml = *(float *) (hsml->data + i*(hsml->strides[0]));
}
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph_sub();
for (i = 0; i < pos->dimensions[0]; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphQ[i].Observable;
}
free(Q);
free(SphQ);
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_Ngbs(PyObject* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *poss;
int i,j,n,nn;
int input_dimension;
size_t bytes;
int startnode,numngb;
int phase=0;
FLOAT searchcenter[3];
double dx,dy,dz,r2,eps2;
input_dimension =pos->nd;
if (input_dimension != 1)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 1");
pos = TO_FLOAT(pos);
eps2 = eps*eps;
searchcenter[0] = (FLOAT)*(float *) (pos->data + 0*(pos->strides[0]));
searchcenter[1] = (FLOAT)*(float *) (pos->data + 1*(pos->strides[0]));
searchcenter[2] = (FLOAT)*(float *) (pos->data + 2*(pos->strides[0]));
startnode = All.MaxPart;
/* ici, il faut faire une fct qui fonctionne en //, cf hydra --> Exportflag */
numngb = ngb_treefind_pairs(&searchcenter[0], (FLOAT)eps, phase, &startnode);
nn=0;
for(n = 0;n < numngb; n++)
{
j = Ngblist[n];
dx = searchcenter[0] - P[j].Pos[0];
dy = searchcenter[1] - P[j].Pos[1];
dz = searchcenter[2] - P[j].Pos[2];
r2 = dx * dx + dy * dy + dz * dz;
if (r2<=eps2)
{
printf("%d r=%g\n",nn,sqrt(r2));
nn++;
}
}
return PyArray_Return(pos);
}
static PyObject * gadget_SphEvaluateOrigAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=NumPart;
vobs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
/* flag all particles as active */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph_orig();
for (i = 0; i < NumPart; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphP[i].Observable;
}
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_SphEvaluateAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
int ld[1];
PyArrayObject *vobs;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
/* create a NumPy object */
ld[0]=NumPart;
vobs = (PyArrayObject *) PyArray_SimpleNew(1,ld,PyArray_FLOAT);
/* flag all particles as active */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph();
for (i = 0; i < NumPart; i++)
{
*(float *)(vobs->data + i*(vobs->strides[0])) = SphP[i].Observable;
}
return Py_BuildValue("O",vobs);
}
static PyObject * gadget_SphEvaluateGradientAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
All.Ti_Current=1; /* need to flag active particles */
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
sph();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject * gadget_DensityEvaluateAll(PyObject* self, PyObject *args)
{
int i;
/* flag all particles as active */
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
density(0);
return Py_BuildValue("i",1);
}
static PyObject * gadget_DensityEvaluateGradientAll(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
density_sph_gradient();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject * gadget_EvaluateThermalConductivity(PyObject* self, PyObject *args)
{
PyArrayObject *obs;
if (! PyArg_ParseTuple(args, "O",&obs))
return PyString_FromString("error");
int i;
size_t bytes;
npy_intp ld[2];
PyArrayObject *grad;
if (obs->nd != 1)
PyErr_SetString(PyExc_ValueError,"dimension of obs must be 1.");
if (obs->dimensions[0] != NumPart)
PyErr_SetString(PyExc_ValueError,"The size of obs must be NumPart.");
obs = TO_FLOAT(obs);
//All.Ti_Current=1; /* need to flag active particles */
for (i = 0; i < NumPart; i++)
P[i].Ti_endstep == All.Ti_Current;
/* now, give observable value for P */
for (i = 0; i < NumPart; i++)
{
SphP[i].Observable = *(float *) (obs->data + i*(obs->strides[0]));
}
density_sph_gradient();
// equivalent of SphP[i].GradEnergyInt[0] is now in
//SphP[i].GradObservable[0]
//SphP[i].GradObservable[1]
//SphP[i].GradObservable[2]
/* now, compute thermal conductivity */
sph_thermal_conductivity();
/* create a NumPy object */
ld[0] = NumPart;
ld[1] = 3;
grad = (PyArrayObject *) PyArray_SimpleNew(2,ld,PyArray_FLOAT);
for (i = 0; i < NumPart; i++)
{
*(float *) (grad->data + i*(grad->strides[0]) + 0*grad->strides[1]) = SphP[i].GradObservable[0];
*(float *) (grad->data + i*(grad->strides[0]) + 1*grad->strides[1]) = SphP[i].GradObservable[1];
*(float *) (grad->data + i*(grad->strides[0]) + 2*grad->strides[1]) = SphP[i].GradObservable[2];
}
return PyArray_Return(grad);
}
static PyObject *gadget_LoadParticles2(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);
/***************************************
* some inits *
/***************************************/
RestartFlag = 0;
Begrun1();
/***************************************
* LOAD PARTILES *
/***************************************/
NumPart = 0;
N_gas = *(int*) (ntype->data + 0*(ntype->strides[0]));
for (i = 0; i < 6; i++)
NumPart += *(int*) (ntype->data + i*(ntype->strides[0]));
if (NumPart!=pos->dimensions[0])
{
PyErr_SetString(PyExc_ValueError,"Numpart != pos->dimensions[0].");
return NULL;
}
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);
/*********************/
/* allocate 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 */
/*********************/
float * fpt;
for (i = 0; i < NumPart; 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].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]);
fpt = (float *) (pos->data + i*(pos->strides[0]) + 0*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]));
//P[i].Active = 1;
}
/***************************************
* END LOAD PARTILES *
/***************************************/
/***************************************
* finish inits *
/***************************************/
Begrun2();
return Py_BuildValue("i",1);
}
#ifdef TESSEL
static PyObject *gadget_ConstructDelaunay(PyObject *self, PyObject *args, PyObject *kwds)
{
ConstructDelaunay();
return Py_BuildValue("i",1);
}
static PyObject *gadget_ComputeVoronoi(PyObject *self, PyObject *args, PyObject *kwds)
{
ComputeVoronoi();
return Py_BuildValue("i",1);
}
#endif /*TESSEL*/
/* definition of the method table */
static PyMethodDef gadgetMethods[] = {
{"Info", (PyCFunction)gadget_Info, METH_VARARGS,
"give some info"},
{"InitMPI", (PyCFunction)gadget_InitMPI, METH_VARARGS,
"Init MPI"},
{"InitDefaultParameters", (PyCFunction)gadget_InitDefaultParameters, METH_VARARGS,
"Init default parameters"},
{"GetParameters", (PyCFunction)gadget_GetParameters, METH_VARARGS,
"get gadget parameters"},
{"SetParameters", (PyCFunction)gadget_SetParameters, METH_VARARGS,
"Set gadget parameters"},
{"check_parser", (PyCFunction)gadget_check_parser, METH_VARARGS|METH_KEYWORDS,
"check the parser"},
{"LoadParticles", (PyCFunction)gadget_LoadParticles, METH_VARARGS|METH_KEYWORDS,
"LoadParticles partilces"},
{"LoadParticlesQ", (PyCFunction)gadget_LoadParticlesQ, METH_VARARGS,
"LoadParticles partilces Q"},
{"LoadParticles2", (PyCFunction)gadget_LoadParticles2, METH_VARARGS,
"LoadParticles partilces"},
{"AllPotential", (PyCFunction)gadget_AllPotential, METH_VARARGS,
"Computes the potential for each particle"},
{"AllAcceleration", (PyCFunction)gadget_AllAcceleration, METH_VARARGS,
"Computes the gravitational acceleration for each particle"},
{"GetAllAcceleration", (PyCFunction)gadget_GetAllAcceleration, METH_VARARGS,
"get the gravitational acceleration for each particle"},
{"GetAllPotential", (PyCFunction)gadget_GetAllPotential, METH_VARARGS,
"get the potential for each particle"},
{"GetAllDensities", (PyCFunction)gadget_GetAllDensities, METH_VARARGS,
"get the densities for each particle"},
{"GetAllHsml", (PyCFunction)gadget_GetAllHsml, METH_VARARGS,
"get the sph smoothing length for each particle"},
{"GetAllPositions", (PyCFunction)gadget_GetAllPositions, METH_VARARGS,
"get the position for each particle"},
{"GetAllVelocities", (PyCFunction)gadget_GetAllVelocities, METH_VARARGS,
"get the velocities for each particle"},
{"GetAllMasses", (PyCFunction)gadget_GetAllMasses, METH_VARARGS,
"get the mass for each particle"},
{"GetAllEnergySpec", (PyCFunction)gadget_GetAllEnergySpec, METH_VARARGS,
"get the specific energy for each particle"},
{"GetAllEntropy", (PyCFunction)gadget_GetAllEntropy, METH_VARARGS,
"get the entropy for each particle"},
{"GetAllID", (PyCFunction)gadget_GetAllID, METH_VARARGS,
"get the ID for each particle"},
{"GetAllTypes", (PyCFunction)gadget_GetAllTypes, METH_VARARGS,
"get the type for each particle"},
{"GetPos", (PyCFunction)gadget_GetPos, METH_VARARGS,
"get the position for each particle (no memory overhead)"},
{"Potential", (PyCFunction)gadget_Potential, METH_VARARGS,
"get the potential for a givent sample of points"},
{"Acceleration", (PyCFunction)gadget_Acceleration, METH_VARARGS,
"get the acceleration for a givent sample of points"},
{"SphEvaluateOrigAll", (PyCFunction)gadget_SphEvaluateOrigAll, METH_VARARGS,
"run the original sph routine."},
{"SphEvaluateAll", (PyCFunction)gadget_SphEvaluateAll, METH_VARARGS,
"compute mean value of a given field based on the sph convolution for all points."},
{"SphEvaluateGradientAll", (PyCFunction)gadget_SphEvaluateGradientAll, METH_VARARGS,
"compute the gradient of a given field based on the sph convolution for all points."},
{"DensityEvaluateAll", (PyCFunction)gadget_DensityEvaluateAll, METH_VARARGS,
"simply run the density function"},
{"DensityEvaluateGradientAll", (PyCFunction)gadget_DensityEvaluateGradientAll, METH_VARARGS,
"run the modified density function and compute a gradient from the given value"},
{"EvaluateThermalConductivity", (PyCFunction)gadget_EvaluateThermalConductivity, METH_VARARGS,
"evaluate the thermal coductivity for each particle"},
{"SphEvaluate", (PyCFunction)gadget_SphEvaluate, METH_VARARGS,
"compute mean value based on the sph convolution for a given number of points"},
{"InitHsml", (PyCFunction)gadget_InitHsml, METH_VARARGS,
"Init hsml based on the three for a given number of points"},
{"Density", (PyCFunction)gadget_Density, METH_VARARGS,
"compute Density based on the three for a given number of points"},
{"Ngbs", (PyCFunction)gadget_Ngbs, METH_VARARGS,
"return the position of the neighbors for a given point"},
{"GetAllPositionsQ", (PyCFunction)gadget_GetAllPositionsQ, METH_VARARGS,
"get the position for each particle Q"},
{"GetAllVelocitiesQ", (PyCFunction)gadget_GetAllVelocitiesQ, METH_VARARGS,
"get the velocities for each particle Q"},
{"GetAllMassesQ", (PyCFunction)gadget_GetAllMassesQ, METH_VARARGS,
"get the mass for each particle Q"},
{"GetAllIDQ", (PyCFunction)gadget_GetAllIDQ, METH_VARARGS,
"get the ID for each particle Q"},
{"GetAllTypesQ", (PyCFunction)gadget_GetAllTypesQ, METH_VARARGS,
"get the type for each particle Q"},
{"Free", (PyCFunction)gadget_Free, METH_VARARGS,
"release memory"},
#ifdef TESSEL
{"ConstructDelaunay", (PyCFunction)gadget_ConstructDelaunay, METH_VARARGS,
"Construct the Delaunay tesselation"},
{"ComputeVoronoi", (PyCFunction)gadget_ComputeVoronoi, METH_VARARGS,
"Compute the Voronoi tesselation"},
{"GetAllDelaunayTriangles", (PyCFunction)gadget_GetAllDelaunayTriangles, METH_VARARGS,
"Get all the Delaunay Triangles"},
{"GetAllvPoints", gadget_GetAllvPoints, METH_VARARGS,
"Get voronoi points"},
{"GetAllvDensities", gadget_GetAllvDensities, METH_VARARGS,
"Get voronoi density of all points"},
{"GetAllvVolumes", gadget_GetAllvVolumes, METH_VARARGS,
"Get voronoi volumes of all points"},
{"GetvPointsForOnePoint", gadget_GetvPointsForOnePoint, METH_VARARGS,
"Get voronoi points for a given point"},
#endif
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initgadget(void)
{
(void) Py_InitModule("gadget", gadgetMethods);
import_array();
}
#endif /*PY_INTERFACE*/
c4science · Help