Skip to content

GitLab

Commit e25e1480 authored by LAMBERT Jean-charles's avatar LAMBERT Jean-charles
Browse files

in progess...

parent a2de79dc
Hide whitespace changes
Inline Side-by-side
lib/projects/nemodep/crectify.cc
View file @ e25e1480
......@@ -34,8 +34,10 @@ using namespace std;
// PUBLIC CRectify
//------------------------------------------------------------------------------
//
CRectify::CRectify()
CRectify::CRectify(bool _verbose)
{
verbose=_verbose;
only_eigen_values=false;
init();
}
......@@ -216,7 +218,6 @@ void CRectify::processRho()
}
}
std::cerr << "Nbody = "<< nbody << " sel=" << ii << "\n";
// compute COD without shifting
if (cod_file.length()==0) { // we must compute COD
......@@ -266,7 +267,6 @@ void CRectify::findMoment()
}
}
DIVMS(w_qpole, w_qpole, w_sum);
std::cerr << " find moment, done\n";
}
//------------------------------------------------------------------------------
// computeVectors
......@@ -275,21 +275,26 @@ void CRectify::computeVectors()
{
// computing Eigen vectors according to moment w_qpole
eigenFrame(frame, w_qpole);
float tmp;
vectutils::dotvp(tmp,oldframe[0], frame[0]);
if (tmp < 0.0)
vectutils::mulvs(frame[0], frame[0], (float) -1.0);
vectutils::dotvp(tmp,oldframe[2], frame[2]);
if (tmp < 0.0)
vectutils::mulvs(frame[2], frame[2], (float )-1.0);
CROSSVP(frame[1], frame[2], frame[0]);
printvec("e_x:", frame[0]);
printvec("e_y:", frame[1]);
printvec("e_z:", frame[2]);
// update old frame vectors
for (int i = 0; i < NDIM; i++)
vectutils::setv(oldframe[i], frame[i]);
if (! only_eigen_values) { // does not compute only eigen values
float tmp;
vectutils::dotvp(tmp,oldframe[0], frame[0]);
if (tmp < 0.0)
vectutils::mulvs(frame[0], frame[0], (float) -1.0);
vectutils::dotvp(tmp,oldframe[2], frame[2]);
if (tmp < 0.0)
vectutils::mulvs(frame[2], frame[2], (float )-1.0);
CROSSVP(frame[1], frame[2], frame[0]);
if (verbose) {
printvec("e_x:", frame[0]);
printvec("e_y:", frame[1]);
printvec("e_z:", frame[2]);
}
// update old frame vectors
for (int i = 0; i < NDIM; i++)
vectutils::setv(oldframe[i], frame[i]);
}
}
//------------------------------------------------------------------------------
// snapTransform STATIC
......
lib/projects/nemodep/crectify.h
View file @ e25e1480
......@@ -23,13 +23,14 @@ http://carma.astro.umd.edu/nemo/man_html/snaprect.1.html
#define CRECTIFY_H
#include <string>
#include "cfalcon.h"
#include <cassert>
#define SINGLEPREC
#define real float
extern "C" {
#include <vectmath.h>
}
#include "csnaptools.h"
#undef real // prevent conflict with SWIG py_unsio.i
namespace uns_proj {
class CDataIndex {
......@@ -54,7 +55,7 @@ public:
class CRectify {
public:
CRectify();
CRectify(bool _verbose=false);
//------------------------------------------------------------------------------
// Rectify
// this method rectify position and vitesse for tge current time, according to
......@@ -76,6 +77,44 @@ public:
const bool _use_rho, const bool rectify=false,
std::string cod_file="",std::string rect_file="", const float radius=50.0, const float dmin=0.0, const float dmax=100.0);
//
// computeEigenVectors : (SWIG interface for python)
// compute eigen vectors on current snapshot
// return bool,array
//
// array is one dimension with 16 elements:
// t codx(3) codv(3) eigen_values(flatten matrix 3x3)
int computeEigenVectors(int n1, float * _pos , int n2, float * _vel,
int n3, float * _mass, int n4, float * _rho,
int n5, float * _rect_array,
const float _time,
const bool _use_rho, const bool _rectify=false,
std::string _cod_file="",
const float _radius=50.0,const float _dmin=0.0, const float _dmax=100.0) {
assert(n5==16);
only_eigen_values=true; // !!!
bool ok=rectify(n1/3,_time,_pos,_vel,_mass,_rho,_use_rho,_rectify,_cod_file,"",_radius,_dmin,_dmax);
// fill up returns array
int ii=0;
_rect_array[ii++] = time;
for (int i=0; i<6; i++) {
_rect_array[ii++] = codf[i];
}
for (int i=0; i<3; i++) {
_rect_array[ii++] = frame[i][0];
_rect_array[ii++] = frame[i][1];
_rect_array[ii++] = frame[i][2];
}
return ok;
}
void process();
//------------------------------------------------------------------------------
// snapTransform STATIC
......@@ -118,6 +157,8 @@ private:
float codf[6];
float frame[3][3];
std::vector <float> vpos,vvel,vmass,vrho;
bool verbose;
bool only_eigen_values; // toggle only eigen values computation
};
}
......
py/lib/CMakeLists.txt
View file @ e25e1480
......@@ -55,9 +55,11 @@ INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/swig
INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_SOURCE_DIR}/../src)
INCLUDE_DIRECTORIES(${NUMPY_INCLUDE_DIRS} ${CMAKE_CURRENT_SOURCE_DIR}/../../lib/utils/nemodep )
INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/../../lib/utils )
INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}/../../lib/projects/nemodep )
INCLUDE_DIRECTORIES(${DEHNEN}/falcON/utils/inc ${DEHNEN}/falcON/inc)
INCLUDE_DIRECTORIES(${UNSIOPATH}/inc)
add_definitions(${OPT} -DNO_CUDA -DfalcON_NEMO -DfalcON_SINGLE )
add_definitions(${OPT} -DNO_CUDA -DfalcON_NEMO -DfalcON_SINGLE -DNOBOOST)
# Set definitions
set (INTERFACE swig/py_unstools.i)
......@@ -80,10 +82,10 @@ ELSE (OSX) # Linux
SET(SOEXT "so")
ENDIF(OSX)
SWIG_ADD_MODULE(${MODULENAME} python ${INTERFACE} ${execpp_sources} ${UNSIOPATH}/lib/libJCLutils.${SOEXT} ${UNSIOPATH}/lib/libunsio.${SOEXT} ${UNSIOPATH}/lib/libnemo.a)
SWIG_ADD_MODULE(${MODULENAME} python ${INTERFACE} ${execpp_sources} ${UNSIOPATH}/lib/libJCLutils.${SOEXT} ${UNSIOPATH}/lib/libJCLprojects.${SOEXT} ${UNSIOPATH}/lib/libunsio.${SOEXT} ${UNSIOPATH}/lib/libnemo.a)
# ${UNSIOPATH}/lib/libJCLprojects.${SOEXT}
SWIG_LINK_LIBRARIES(${MODULENAME} ${PYTHON_LIBRARIES} ${UNSIOPATH}/lib/libJCLutils.${SOEXT} ${UNSIOPATH}/lib/libunsio.${SOEXT} ${UNSIOPATH}/lib/libWDutils.${SOEXT} ${UNSIOPATH}/lib/libfalcON.${SOEXT} ${UNSIOPATH}/lib/libnemo.a ${SQLITE3_LIB_PATH}/libsqlite3.${SOEXT} )
SWIG_LINK_LIBRARIES(${MODULENAME} ${PYTHON_LIBRARIES} ${UNSIOPATH}/lib/libJCLutils.${SOEXT} ${UNSIOPATH}/lib/libJCLprojects.${SOEXT} ${UNSIOPATH}/lib/libunsio.${SOEXT} ${UNSIOPATH}/lib/libWDutils.${SOEXT} ${UNSIOPATH}/lib/libfalcON.${SOEXT} ${UNSIOPATH}/lib/libnemo.a ${UNSIOPATH}/lib/libcpgplot.a ${UNSIOPATH}/lib/libpgplot.so ${SQLITE3_LIB_PATH}/libsqlite3.${SOEXT} X11 gfortran)
# ----------------------------------------------------------
# Install SETUP
# ----------------------------------------------------------
......
py/lib/swig/py_unstools.i
View file @ e25e1480
//-*- C -*-
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2008-2015
// Copyright Jean-Charles LAMBERT - 2008-2017
// Centre de donneeS Astrophysiques de Marseille (CeSAM)
// e-mail: Jean-Charles.Lambert@lam.fr
// address: Aix Marseille Universite, CNRS, LAM
......@@ -20,6 +20,8 @@
#define SWIG_FILE_WITH_INIT
#include "cfalcon.h"
#include "ctree.h"
#include "crectify.h"
#include "c2dplot.h"
%}
%include "numpy.i"
......@@ -45,7 +47,7 @@
%inline %{
typedef float real; // real is not known
typedef float real; // real is not known
%}
// Parse the original header file
......@@ -118,11 +120,59 @@
%template(CTreeF) jcltree::CTree<float>;
%template(CTreeD) jcltree::CTree<double>;
// static bool addGravity(const int nbody,
// const float * pos, const float * mass,
// float * acc, float * phi,
// const float eps,
// const float G=1.0,
// const float theta=0.6,
// const int kernel_type=1,
// const int ncrit=6);
//
// -- CRectify class
//
%apply (int DIM1 , float * INPLACE_ARRAY1) {
(int n1, float * _pos ),
(int n2, float * _vel ),
(int n3, float * _mass),
(int n4, float * _rho )
}
%apply (int DIM1 , float * ARGOUT_ARRAY1) {
(int n5, float * _rect_array )
}
%apply (int DIM1 , float * INPLACE_ARRAY1) {
(int nn1, float * ppos )
}
%include "crectify.h"
%extend uns_proj::CRectify {
}
//
// -- C2dplot class
//
%apply (int DIM1, float * IN_ARRAY1) {(int n1, float * pos),(int n2, float * weight),(int n3, float * hsml),(int n4, float * _range)}
%apply (int DIM1, double* IN_ARRAY1) {(int n1, double* pos),(int n2, double* weight),(int n3, double* hsml)}
%include "c2dplot.h"
%extend uns_proj::C2dplot {
bool compute_image(std::string _dev, const int _no_frame,int n1, T * pos , int n4, float * _range,
std::string _title,std::string _sel_comp,std::string _filename, float _timu,
bool _xy, bool _xz, bool _zy, bool _sview, int n2, T * weight,
const int _psort, int n3, T * hsml, const int _itf, const bool _wedge,
std::string _legend, const int _cmap) {
// float _range[3][2]
T *p_hsml = NULL;
if (n3 !=0) {
p_hsml=hsml;
}
T *p_weight = NULL;
if (n2 !=0) {
p_weight=weight;
}
$self->compute(_dev,_no_frame,n1/3,pos,(float (*)[2]) _range,_title,_sel_comp,_filename,
_timu,_xy,_xz,_zy,_sview, p_weight, _psort, p_hsml, _itf, _wedge, _legend,_cmap);
}
}
%template(C2dplotF) uns_proj::C2dplot<float >;
%template(C2dplotD) uns_proj::C2dplot<double>;
//
py/mains/g3cpu_parse.py
View file @ e25e1480
......@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
# cnd line
import sys, getopt,os.path
sys.path.append('/home/jcl/works/SVN/uns_projects/trunk/py/modules/')
#sys.path.append('/home/jcl/works/SVN/uns_projects/trunk/py/modules/')
from uns_simu import *
import argparse
......
py/mains/process_analysis.py
View file @ e25e1480
#!/usr/bin/env python
import numpy as np
import os,time
import sys
from multiprocessing import Process
......@@ -40,7 +39,7 @@ def commandLine():
# process, is the core function
def process(args):
try:
analysis=CUnsAnalysis(simname=args.simname, script=args.script,verbose_debug=True)
analysis=CUnsAnalysis(simname=args.simname, script=args.script,verbose_debug=args.verbose)
except Exception as x :
print (x.message)
else:
......
py/mains/script_analysis.txt
View file @ e25e1480
#!/usr/bin/env python
"""
This script is called by external python script through execfile() call
"""
from simulations.ccod import *
from simulations.ccom import *
from simulations.crectify import *
##
## COD
##
if 1:
data.lock_id=0 # lock file ID for COD analysis
data.cod_dir="ANALYSIS/"+data.simname+"/cod"
data.cod_select="halo: disk: stars: gas: bndry: bulge: halo_1: halo_2: halo,disk: halo,stars: disk,stars: halo,disk,stars"
if data.first:
cod=CCod("",analysis=data,verbose_debug=True) # initialyze COD analysis
cod.smartAnalysis(analysis=data)
##
## COD
## COM
##
data.lock_id=0 # lock file ID for COD analysis
data.cod_dir="ANALYSIS/cod"
data.cod_select="halo:disk:stars:gas:bndry:bulge:halo_1:halo_2:halo,disk:halo,stars:disk,stars:halo,disk,stars"
if 1:
data.lock_id=1
data.com_dir="ANALYSIS/"+data.simname+"/com"
data.com_select="halo: disk: stars: gas: disk,stars: halo,disk,stars"
if data.first:
cod=CCod("",analysis=data,verbose_debug=True) # initialyze COD analysis
if data.first:
com=CCom(analysis=data,verbose_debug=True) # initialyze COM analysis
cod.smartAnalysis(analysis=data)
com.smartAnalysis()
##
## RECTIFY
##
if 1:
data.lock_id=2
data.rectify_dir="ANALYSIS/"+data.simname+"/rectify"
data.rectify_select="stars#halo,disk,stars: stars#stars: gas#stars@10.0 " # syntax: "component#cod@rcut"
if data.first:
rectify=CRectify(analysis=data,verbose_debug=True) # initialyze RECTIFY analysis
rectify.smartAnalysis()
##
data.first=False #MANDATORY
......
py/modules/simulations/ccod.py
View file @ e25e1480
......@@ -557,7 +557,7 @@ class CCod:
start some initialisations
"""
### build COD Dir
### COD Dir name
if hasattr(data,'cod_dir'):
pass
#
......@@ -567,19 +567,7 @@ class CCod:
print("CoD DIR = ", data.cod_dir, data.sim_info['name'])
self.__cod_file_base=data.cod_dir
# lock process
data.lock[data.lock_id].acquire()
if (not os.path.isdir(data.cod_dir)) :
try:
print("Core ID ",data.core_id," create directory [%s]\n"%(data.cod_dir))
os.makedirs(data.cod_dir)
except OSError:
print("Unable to create directory [%s]\n"%(data.cod_dir))
data.lock[data.lock_id].release()
sys.exit()
# release process
data.lock[data.lock_id].release()
data.cod_select=data.cod_select.replace(" ", "")
### re build select component variable according to components existing at mid
### simulation time (pre-computed by cuns_analysis.py and set to data.list_components
......@@ -604,6 +592,69 @@ class CCod:
if self.__vdebug:
print("new select :",self.__new_select)
#
# build directories and links
#
# lock process
data.lock[data.lock_id].acquire()
if (not os.path.isdir(data.cod_dir)) : # build dir
try:
print("Core ID ",data.core_id," create directory [%s]\n"%(data.cod_dir))
os.makedirs(data.cod_dir)
except OSError:
print("Unable to create directory [%s]\n"%(data.cod_dir))
data.lock[data.lock_id].release()
sys.exit()
# make links
cod_basename=data.cod_dir+"/"+data.sim_info['name']
simname=data.sim_info['name']
if self.__vdebug:
print("Make links, cod_basename[%s]\n"%(cod_basename))
if self.__new_select.find('gas') ==-1 and \
self.__new_select.find('disk')!=-1 : # gas no and disk yes
mysrc=simname+".disk.cod"
mylink=cod_basename+".disk,stars.cod"
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
mysrc=simname+".halo,disk.cod"
mylink=cod_basename+".halo,disk,stars.cod"
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
if self.__new_select.find('gas') !=-1 and \
self.__new_select.find('disk')==-1 : # gas yes and disk no
mysrc=simname+".stars.cod"
mylink=cod_basename+".disk,stars.cod"
if self.__vdebug:
print ("Links : src[%s] dest[%s]\n"%(mysrc,mylink))
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
mysrc=simname+".halo,stars.cod"
mylink=cod_basename+".halo,disk,stars.cod"
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
if self.__new_select.find('halo_1') !=-1: # halo_1 yes
mysrc=simname+".halo_1.cod"
mylink=cod_basename+"-gal01.halo.cod"
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
if self.__new_select.find('halo_2') !=-1: # halo_2 yes
mysrc=simname+".halo_2.cod"
mylink=cod_basename+"-gal02.halo.cod"
if not os.path.islink(mylink):
os.symlink(mysrc,mylink)
# release process
data.lock[data.lock_id].release()
#
# smartAnalysis
#
......
py/modules/simulations/csnapshot.py
View file @ e25e1480
......@@ -7,6 +7,8 @@ from uns_simu import *
from cfalcon import *
import numpy as np
#
# class CSsnapshot
#
class CSnapshot:
"""Operations on UNS snapshots"""
......@@ -15,7 +17,7 @@ class CSnapshot:
__float32=True
simname=None
select=None
def __init__(self,simname,select="all",times="all",float32=True,verbose=False,verbose_debug=False):
select=select.encode('ascii')
times=times.encode('ascii')
......@@ -39,7 +41,9 @@ class CSnapshot:
raise RuntimeError("UNS not valid")
else:
None
#
# nextFrame
#
def nextFrame(self,bits=""):
"""
......@@ -58,17 +62,52 @@ class CSnapshot:
raise RuntimeError("UNS not valid")
else:
return self.__uns.nextFrame(bits)
#
# close
#
def close(self):
if not self.__uns.isValid():
raise RuntimeError("UNS not valid")
self.__uns.close()
#
# kill
#
def kill(self):
if not self.__uns.isValid():
raise RuntimeError("UNS not valid")
self.close()
#
# getInterfaceType
#
def getInterfaceType(self):
"""
return uns snapshot's Interface type
"""
if not self.__uns.isValid():
raise RuntimeError("UNS not valid")
return self.__uns.getInterfaceType()
#
# getFileStructure
#
def getFileStructure(self):
"""
return uns snapshot's File structure
"""
if not self.__uns.isValid():
raise RuntimeError("UNS not valid")
return self.__uns.getFileStructure()
#
# getFileName
#
def getFileName(self):
"""
return uns snapshot's File name
"""
if not self.__uns.isValid():
raise RuntimeError("UNS not valid")
return self.__uns.getFileName()
#
# getData
#
def getData(self,select,tag=None,data_type=np.float32):
"""
......@@ -105,7 +144,9 @@ class CSnapshot:
if self.__vdebug:
print ("ok=",ok," data=",data.size," ret_data=",ret_data.size)
return status,ret_data
#
# __getData
#
def __getData(self,comp_value,tag=None,data_type=np.float32):
if self.__vdebug:
print ("array=",tag,data_type)
......@@ -132,6 +173,8 @@ class CSnapshot:
ok,data=self.__uns.getArrayI(comp_value,tag)
return ok,2,data
#
# center
#
def center(self,pos,vel,data_weight,center=False):
"""
......@@ -154,7 +197,9 @@ class CSnapshot:
if vel is not None:
cxv[3:6]=self.__centerOnWeight(vel,data_weight,center)
return cxv
#
# __centerOnWeight
#
def __centerOnWeight(self,data,data_weight,center):
# reshape array in x,y,z arrays
data=np.reshape(data,(-1,3)) # data reshaped in a 2D array [nbody,3]
......@@ -173,6 +218,8 @@ class CSnapshot:
data[:,2] -= zcom
return xcom,ycom,zcom
#
# computeDensity
#
def computeDensity(self,pos,mass,K=32,N=1,method=0,ncrit=None):
"""
......
py/modules/simulations/ctree.py
View file @ e25e1480
......@@ -15,11 +15,13 @@ class CTree:
__time=0.0
# -----------------------------------------------------
#
def __init__(self,pos,vel=None,mass=None,fcells=0.9,rsize=0.4,time=0.0):
def __init__(self,pos,vel=None,mass=None,fcells=0.9,rsize=0.4,time=0.0,seed=None):
self.__pos=pos
self.__mass=mass
self.__vel=vel
self.__time=float(time)
# initialyse random generator number
np.random.seed(seed)
if (type(pos[0])==np.float32):
self.__tree=ut.CTreeF(pos,mass,fcells,rsize)
......
py/modules/simulations/cuns_analysis.py
View file @ e25e1480
......@@ -24,6 +24,15 @@ class CUnsAnalysis:
__analysis_script=None
__select="all"
__nb_lock=15 # #lock files
class CData:
"""
store information used during analysis process
"""
uns_snap=None # pointer to UNS snpshot
first=True
lock=[]
# constructor
def __init__(self,simname,script,dbname=None,verbose=False,verbose_debug=False):
"""
......@@ -49,6 +58,17 @@ class CUnsAnalysis:
else:
self.__detectComponents()
# create data analysis object
self.__data=self.CData()
self.__data.list_components=self.__new_select
self.__data.simname=self.simname
self.__data.sim_info=self.__r # unsio simulation info
self.__data.lock=[] # locking mechanism
for i in range