Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77958152
pnbody.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 17, 11:45
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 11:45 (2 d)
Engine
blob
Format
Raw Data
Handle
19956338
Attached To
rGEAR Gear
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
Log In to Comment