Page MenuHomec4science
Files F77958152

pnbody.c
No OneTemporary
Actions

File Metadata

Created
Sat, Aug 17, 11:45

pnbody.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#ifdef PNBODY
#include <Python.h>
#include <numpy/arrayobject.h>
#include "allvars.h"
#include "proto.h"
double NextTime;
PyObject *ImageModule;
PyObject *ImageDict;
PyObject *ImageMaker;
PyObject *pyCreator;
PyObject *pyArgs;
PyObject *pyFile;
PyObject *pyNextTime;
PyArrayObject *Pos;
PyArrayObject *Vel;
PyArrayObject *Num;
PyArrayObject *Mass;
PyArrayObject *Tpe;
PyArrayObject *U;
PyArrayObject *Rho;
PyArrayObject *Rsp;
PyArrayObject *Metals;
void init_pnbody()
{
/* initialize python and numpy */
Py_Initialize();
import_array();
ImageModule = PyImport_ImportModule("Mkgmov");
/* create the object */
pyCreator = PyObject_GetAttrString(ImageModule, "Movie");
/* file */
pyFile = PyString_FromString("filmparam.py");
pyArgs = PyTuple_New(1);
PyTuple_SetItem(pyArgs, 0, pyFile);
//ImageMaker = PyObject_CallObject(pyCreator, NULL);
ImageMaker = PyObject_CallObject(pyCreator, pyArgs);
/* deallocate */
Py_DECREF(pyCreator);
Py_DECREF(pyArgs);
Py_DECREF(pyFile);
/* dictionary */
ImageDict =PyDict_New();
/* give some info */
PyObject_CallMethod(ImageMaker, "info",NULL);
}
/*! finalize pNbody
*
*/
void finalize_pnbody()
{
/* clean */
Py_DECREF(ImageDict);
Py_DECREF(ImageMaker);
Py_Finalize();
}
void compute_pnbody()
{
/****************************************/
/* check if it is time for the next image
/****************************************/
pyNextTime = PyObject_CallMethod(ImageMaker, "get_next_time",NULL);
if (!PyFloat_Check(pyNextTime))
{
printf("no more time left.\n");
Py_DECREF(pyNextTime);
return;
}
NextTime = PyFloat_AsDouble(pyNextTime);
printf(" Time=%g NextTime=%g\n",All.Time,NextTime);
if (NextTime<=All.Time)
{
/* set next time */
PyObject_CallMethod(ImageMaker, "set_next_time",NULL);
/***********************************/
/* create arrays
/***********************************/
/*
need to add :
- metals (gas and stars)
- tstar (stars)
- minit (stars) (bof)
- idp (stars) (bof)
- rho_stars (stars) (bof)
*/
int i,k;
npy_intp ld1[2],ld2[1],ldn[2];
FLOAT Ui,Rhoi,Rspi;
ld2[0]=NumPart;
ld2[1]=3;
ld1[0]=NumPart;
ldn[0]=NumPart;
ldn[1]=NELEMENTS;
Pos = (PyArrayObject *) PyArray_SimpleNew(2,ld2,NPY_FLOAT);
Vel = (PyArrayObject *) PyArray_SimpleNew(2,ld2,NPY_FLOAT);
Num = (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_INT);
Mass= (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_FLOAT);
Tpe = (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_INT);
U = (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_FLOAT);
Rho = (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_FLOAT);
Rsp = (PyArrayObject *) PyArray_SimpleNew(1,ld1,NPY_FLOAT);
Metals = (PyArrayObject *) PyArray_SimpleNew(2,ldn,NPY_FLOAT);
for (i=0;i<NumPart;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] ;
*(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] ;
*(int*) (Num->data + i*(Num->strides[0])) = P[i].ID;
*(int*) (Tpe->data + i*(Tpe->strides[0])) = P[i].Type;
*(float*)(Mass->data + i*(Mass->strides[0])) = P[i].Mass;
if (i<N_gas)
{
#ifdef ISOTHERM_EQS
Ui = P[i].Entropy ;
#else
Ui = dmax(All.MinEgySpec,SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].Density , GAMMA_MINUS1));
//Ui = dmax(All.MinEgySpec,SphP[i].Entropy / GAMMA_MINUS1 * pow(SphP[i].Density * a3inv, GAMMA_MINUS1));
#endif
Rhoi =SphP[i].Density ;
Rspi = SphP[i].Hsml ;
}
else
{
Ui = 0;
Rhoi=0;
Rspi=0;
}
*(float*)(U->data + i*( U->strides[0])) = Ui ;
*(float*)(Rho->data + i*(Rho->strides[0])) = Rhoi ;
*(float*)(Rsp->data + i*(Rsp->strides[0])) = Rspi ;
switch (P[i].Type)
{
case 0:
for(k=0;k<NELEMENTS;k++)
*(float*)(Metals->data + i*(Metals->strides[0]) + k*Metals->strides[1]) = SphP[i].Metal[k];
break;
case ST:
for(k=0;k<NELEMENTS;k++)
*(float*)(Metals->data + i*(Metals->strides[0]) + k*Metals->strides[1]) = StP[P[i].StPIdx].Metal[k];
break;
default:
for(k=0;k<NELEMENTS;k++)
*(float*)(Metals->data + i*(Metals->strides[0]) + k*Metals->strides[1]) = 0;
break;
}
}
/***********************************/
/* create image
/***********************************/
/* fill the dictionnary */
PyDict_Clear(ImageDict);
PyDict_SetItemString(ImageDict,"pos", (PyObject*)Pos);
PyDict_SetItemString(ImageDict,"vel", (PyObject*)Vel);
PyDict_SetItemString(ImageDict,"num", (PyObject*)Num);
PyDict_SetItemString(ImageDict,"mass", (PyObject*)Mass);
PyDict_SetItemString(ImageDict,"tpe", (PyObject*)Tpe);
PyDict_SetItemString(ImageDict,"u", (PyObject*)U);
PyDict_SetItemString(ImageDict,"rho", (PyObject*)Rho);
PyDict_SetItemString(ImageDict,"rsp", (PyObject*)Rsp);
PyDict_SetItemString(ImageDict,"metals", (PyObject*)Metals);
PyDict_SetItemString(ImageDict,"atime", PyFloat_FromDouble((double)All.Time));
/* call the dump method */
PyObject_CallMethod(ImageMaker, "dump","O",ImageDict);
Py_DECREF(Pos);
Py_DECREF(Vel);
Py_DECREF(Num);
Py_DECREF(Mass);
Py_DECREF(Tpe);
Py_DECREF(U);
Py_DECREF(Rho);
Py_DECREF(Rsp);
Py_DECREF(Metals);
}
Py_DECREF(pyNextTime);
}
#endif

Event Timeline

c4science · Help