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

many fixes

parent 59c5047a
...@@ -168,10 +168,15 @@ void CBar::save(std::string out, const float timu, const bool mvcod) ...@@ -168,10 +168,15 @@ void CBar::save(std::string out, const float timu, const bool mvcod)
// save data // save data
uns::CunsOut * unsout = new uns::CunsOut(out,"nemo",false); uns::CunsOut * unsout = new uns::CunsOut(out,"nemo",false);
unsout->snapshot->setData("time",timu); unsout->snapshot->setData("time",timu);
if (mass)
unsout->snapshot->setData("mass",nbody,mass,false); unsout->snapshot->setData("mass",nbody,mass,false);
if (pos)
unsout->snapshot->setData("pos" ,nbody,pos ,false); unsout->snapshot->setData("pos" ,nbody,pos ,false);
if (vel)
unsout->snapshot->setData("vel" ,nbody,vel ,false); unsout->snapshot->setData("vel" ,nbody,vel ,false);
if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false); unsout->snapshot->setData("rho" ,nbody,rho ,false);
if (hsml)
unsout->snapshot->setData("hsml",nbody,hsml,false); unsout->snapshot->setData("hsml",nbody,hsml,false);
if (id) if (id)
unsout->snapshot->setData("id" ,nbody,id ,false); unsout->snapshot->setData("id" ,nbody,id ,false);
......
...@@ -163,6 +163,7 @@ void CRectify::processRho() ...@@ -163,6 +163,7 @@ void CRectify::processRho()
{ {
if (!rho) { if (!rho) {
std::cerr << "Computing rho........\n";
// Instantiate a density object // Instantiate a density object
density = new CDensity(nbody,&pos[0],&mass[0]); density = new CDensity(nbody,&pos[0],&mass[0]);
density->compute(0,32,1,8); // estimate density density->compute(0,32,1,8); // estimate density
......
...@@ -97,7 +97,7 @@ public: ...@@ -97,7 +97,7 @@ public:
assert(n5==16); assert(n5==16);
only_eigen_values=true; // !!! only_eigen_values=true; // !!!
std::cerr << "Use rho="<<_use_rho << " dmin=" << _dmin << " dmax=" << _dmax << "\n";
bool ok=rectify(n1/3,_time,_pos,_vel,_mass,_rho,_use_rho,_rectify,_cod_file,"",_radius,_dmin,_dmax); bool ok=rectify(n1/3,_time,_pos,_vel,_mass,_rho,_use_rho,_rectify,_cod_file,"",_radius,_dmin,_dmax);
// fill up returns array // fill up returns array
......
...@@ -35,6 +35,7 @@ def commandLine(): ...@@ -35,6 +35,7 @@ def commandLine():
#parser.add_argument('--no-plot', help='Disable ploting',dest='plot',action='store_false',default=True) #parser.add_argument('--no-plot', help='Disable ploting',dest='plot',action='store_false',default=True)
parser.add_argument('--dbname',help='UNS database file name', default=dbname) parser.add_argument('--dbname',help='UNS database file name', default=dbname)
parser.add_argument('--verbose',help='verbose mode on', dest="verbose", action="store_true", default=False) parser.add_argument('--verbose',help='verbose mode on', dest="verbose", action="store_true", default=False)
parser.add_argument('--saveanalysis',help='save merging time in ANALYSIS directory', dest="store_analysis", action="store_true", default=False)
#parser.add_argument('--no-verbose',help='verbose mode off', dest="verbose", action="store_false", default=True) #parser.add_argument('--no-verbose',help='verbose mode off', dest="verbose", action="store_false", default=True)
# parse # parse
...@@ -51,7 +52,7 @@ def commandLine(): ...@@ -51,7 +52,7 @@ def commandLine():
def process(args): def process(args):
cod = CCod(simname=args.simname,verbose_debug=args.verbose,dbname=args.dbname) cod = CCod(simname=args.simname,verbose_debug=args.verbose,dbname=args.dbname)
try: try:
mt=cod.computeMergingTime(halo_1=args.cod1,halo_2=args.cod2,pngfile=args.pngfile,txtfile=args.txtfile,dmax=args.dmax,seger=args.seger,plot=args.plot) mt=cod.computeMergingTime(halo_1=args.cod1,halo_2=args.cod2,pngfile=args.pngfile,txtfile=args.txtfile,dmax=args.dmax,seger=args.seger,plot=args.plot,store_analysis=args.store_analysis)
if mt<0.0: if mt<0.0:
mt=0.0 mt=0.0
print("%f\n"%(mt)) print("%f\n"%(mt))
......
#!/usr/bin/env python #!/usr/bin/env python
from __future__ import print_function
import os,time import os,time
import sys import sys
from multiprocessing import Process from multiprocessing import Process
...@@ -6,6 +7,8 @@ import Queue ...@@ -6,6 +7,8 @@ import Queue
import multiprocessing import multiprocessing
import numpy as np import numpy as np
import argparse import argparse
import matplotlib
matplotlib.use('Agg')
#sys.path=['/home/jcl/works/GIT/uns_projects/py/modules/','/home/jcl/works/GIT/uns_projects/py/modules/simulations']+sys.path #sys.path=['/home/jcl/works/GIT/uns_projects/py/modules/','/home/jcl/works/GIT/uns_projects/py/modules/simulations']+sys.path
from simulations.cuns_analysis import * from simulations.cuns_analysis import *
...@@ -28,6 +31,8 @@ def commandLine(): ...@@ -28,6 +31,8 @@ def commandLine():
parser.add_argument('--dbname',help='UNS database file name', default=dbname) parser.add_argument('--dbname',help='UNS database file name', default=dbname)
parser.add_argument('--verbose',help='verbose mode', default=False) parser.add_argument('--verbose',help='verbose mode', default=False)
print ("Matplotlib backend Using:",matplotlib.get_backend(),file=sys.stderr)
# parse # parse
args = parser.parse_args() args = parser.parse_args()
......
...@@ -48,6 +48,7 @@ def process(args): ...@@ -48,6 +48,7 @@ def process(args):
ok,t=s.getData("time") ok,t=s.getData("time")
if ok: if ok:
print("Latest time [%f]\n"%(t),file=sys.stderr) print("Latest time [%f]\n"%(t),file=sys.stderr)
print("%f"%(t))
if args.save: if args.save:
filename=args.filename filename=args.filename
if args.filename is None: #simulation dir if args.filename is None: #simulation dir
...@@ -62,7 +63,7 @@ def process(args): ...@@ -62,7 +63,7 @@ def process(args):
f.write("last_time %f\n"%(t)) f.write("last_time %f\n"%(t))
f.close() f.close()
except: except:
print("\n\nUnable to save in <%s>, aborting....\n"%(filename)) print("\n\nUnable to save in <%s>, aborting....\n"%(filename),file=sys.stderr)
sys.exit() sys.exit()
# ----------------------------------------------------- # -----------------------------------------------------
......
...@@ -46,6 +46,7 @@ class C2dplot: ...@@ -46,6 +46,7 @@ class C2dplot:
self.__verbose=verbose self.__verbose=verbose
if analysis is not None: if analysis is not None:
matplotlib.use('Agg')
self.__analysis=analysis self.__analysis=analysis
self.__smartAnalysisInit() self.__smartAnalysisInit()
...@@ -104,6 +105,8 @@ class C2dplot: ...@@ -104,6 +105,8 @@ class C2dplot:
dir_movie=dirmovie dir_movie=dirmovie
else: else:
dir_movie=r['dir']+'/ANALYSIS/movie2' # movie directory dir_movie=r['dir']+'/ANALYSIS/movie2' # movie directory
if self.__vdebug:
print("Dir movie =",dir_movie)
for d in ["contour","std"]: for d in ["contour","std"]:
for s in glob.glob(dir_movie+"/"+d+"/*"): for s in glob.glob(dir_movie+"/"+d+"/*"):
dir_frame=os.path.basename(s) dir_frame=os.path.basename(s)
...@@ -487,8 +490,13 @@ class C2dplot: ...@@ -487,8 +490,13 @@ class C2dplot:
else: else:
outfile=out+"."+"%05d"%cpt+".jpg" outfile=out+"."+"%05d"%cpt+".jpg"
if self.__vdebug: if self.__vdebug:
print (">> ",outfile) print (">>>>> ",outfile)
try:
plt.savefig(outfile, bbox_inches=0,dpi=fig.dpi) plt.savefig(outfile, bbox_inches=0,dpi=fig.dpi)
except:
if not os.path.isfile(outfile):
print("\nUnable to save figure [%s], aborting...\n"%(outfile),file=sys.stderr)
sys.exit()
plt.close(fig) plt.close(fig)
......
...@@ -179,7 +179,7 @@ class CCod: ...@@ -179,7 +179,7 @@ class CCod:
test_snap.close() test_snap.close()
self.__is_multiple_halo = self.__getHaloParticlesNumber() # self.__is_multiple_halo = self.__getHaloParticlesNumber() #
self.__startProcess(select,ncores,cod_file) self.__startProcess(select,ncores,cod_file)
self.computeMergingTime() self.computeMergingTime(outdir=None,store_analysis=True)
# #
# ---- # ----
...@@ -400,7 +400,8 @@ class CCod: ...@@ -400,7 +400,8 @@ class CCod:
if self.__vdebug: if self.__vdebug:
print (sline,ind,type(self.__halo_part),file=sys.stderr) print (sline,ind,type(self.__halo_part),file=sys.stderr)
self.__halo_part[ind-1] = int(sline[1]) self.__halo_part[ind-1] = int(sline[1])
else:
print("\n\n!!! File [%s] does not exist !!!!\n\n"%(model_param),file=sys.stderr)
self.__halo_part=self.__halo_part[np.where(self.__halo_part>0)] self.__halo_part=self.__halo_part[np.where(self.__halo_part>0)]
if self.__vdebug: if self.__vdebug:
print (self.__halo_part,self.__halo_part.size,file=sys.stderr) print (self.__halo_part,self.__halo_part.size,file=sys.stderr)
...@@ -411,7 +412,7 @@ class CCod: ...@@ -411,7 +412,7 @@ class CCod:
# #
# ---- # ----
# Compute mering time between two halos galaxy # Compute mering time between two halos galaxy
def computeMergingTime(self,halo_1=None,halo_2=None,simname=None,outdir="./",txtfile="merging_time.txt",pngfile="merging_time.png",dmax=1.0,seger=False, plot=True): def computeMergingTime(self,halo_1=None,halo_2=None,simname=None,outdir="./",txtfile="merging_time.txt",pngfile="merging_time.png",dmax=1.0,seger=False, plot=True,store_analysis=False):
merging_time=-1 merging_time=-1
...@@ -422,6 +423,7 @@ class CCod: ...@@ -422,6 +423,7 @@ class CCod:
if not seger: if not seger:
c_halo_1=self.__cod_file_base+"/"+self.simname+".halo_1.cod" c_halo_1=self.__cod_file_base+"/"+self.simname+".halo_1.cod"
c_halo_2=self.__cod_file_base+"/"+self.simname+".halo_2.cod" c_halo_2=self.__cod_file_base+"/"+self.simname+".halo_2.cod"
if outdir is None:
outdir=self.__cod_file_base outdir=self.__cod_file_base
simname=self.simname simname=self.simname
else: # seger else: # seger
...@@ -470,6 +472,12 @@ class CCod: ...@@ -470,6 +472,12 @@ class CCod:
f.write(str(merging_time)+"\n") f.write(str(merging_time)+"\n")
f.close() f.close()
if store_analysis and self.__r is not None:
m_txt=self.__r["dir"]+"/ANALYSIS/merging_time.txt"
f=open(m_txt,"w")
f.write("merging_time "+str(merging_time)+"\n")
f.close()
if plot: if plot:
# plot # plot
# first plot: plot distance between 2 halos on whole simulation # first plot: plot distance between 2 halos on whole simulation
......
...@@ -473,7 +473,7 @@ class CInert: ...@@ -473,7 +473,7 @@ class CInert:
self.__fdesc.close() self.__fdesc.close()
# #
# __checkTimeExist # __saveNewTime
# #
def __saveNewTime(self,time,t_file): def __saveNewTime(self,time,t_file):
""" """
......
...@@ -22,11 +22,14 @@ class CRectify: ...@@ -22,11 +22,14 @@ class CRectify:
Rectify UNS snapshots Rectify UNS snapshots
""" """
__analysis=None __analysis=None
__use_rho=False
__dmin=0.
__dmax=100.
# #
# ---- # ----
# #
# constructor # constructor
def __init__(self,analysis=None,verbose=False,verbose_debug=False): def __init__(self,analysis=None,use_rho=False,dmin=0.,dmax=100.,verbose=False,verbose_debug=False):
""" """
Constructor of CRectify class Constructor of CRectify class
...@@ -34,6 +37,9 @@ class CRectify: ...@@ -34,6 +37,9 @@ class CRectify:
""" """
self.__vdebug=verbose_debug self.__vdebug=verbose_debug
self.__verbose=verbose self.__verbose=verbose
self.__use_rho=use_rho
self.__dmin=dmin
self.__dmax=dmax
if analysis is not None: if analysis is not None:
self.__analysis=analysis self.__analysis=analysis
...@@ -174,12 +180,14 @@ class CRectify: ...@@ -174,12 +180,14 @@ class CRectify:
ok1,pos = uns_snap.getData(select,"pos" ) ok1,pos = uns_snap.getData(select,"pos" )
ok2,vel = uns_snap.getData(select,"vel" ) ok2,vel = uns_snap.getData(select,"vel" )
ok3,mass = uns_snap.getData(select,"mass") ok3,mass = uns_snap.getData(select,"mass")
ok4,rho = uns_snap.getData(select,"rho")
#print("ok1 ok2 ok3 ",ok1,ok2,ok3) #print("ok1 ok2 ok3 ",ok1,ok2,ok3)
if pos.size>1 and vel.size>1 and mass.size>1 : if pos.size>1 and vel.size>1 and mass.size>1 :
if self.__vdebug: if self.__vdebug:
print("Rectifing time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file)) print("Rectifing time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file))
try: try:
ok,eigen_v=c.computeEigenVectors(pos,vel,mass,mass,16,time.item(),False,False,str(cod_file),rcut) print(self.__use_rho,self.__dmin, self.__dmax)
ok,eigen_v=c.computeEigenVectors(pos,vel,mass,rho,16,time.item(),self.__use_rho,False,str(cod_file),rcut,self.__dmin, self.__dmax)
print("Time [%f] eigen:\n"%(time),eigen_v) print("Time [%f] eigen:\n"%(time),eigen_v)
data.lock[data.lock_id].acquire() data.lock[data.lock_id].acquire()
self.__saveEigenVectors(pre_rect_file,eigen_v) self.__saveEigenVectors(pre_rect_file,eigen_v)
...@@ -264,7 +272,7 @@ class CRectify: ...@@ -264,7 +272,7 @@ class CRectify:
"] in UNS database..." "] in UNS database..."
raise Exception(message) raise Exception(message)
if ev_in is None: if ev_in is None:
simdir=self.__r['dir']+"/ANALYSIS/rectify/" simdir=self.__r['dir']+"/ANALYSIS/rectify2/"
for ev_in in glob.glob(simdir+"*.ev"): for ev_in in glob.glob(simdir+"*.ev"):
rect_out=os.path.basename(ev_in) rect_out=os.path.basename(ev_in)
rect_out=simdir+rect_out.split(".ev")[0]+".rect" rect_out=simdir+rect_out.split(".ev")[0]+".rect"
......
...@@ -114,10 +114,10 @@ int main(int argc, char ** argv ) ...@@ -114,10 +114,10 @@ int main(int argc, char ** argv )
assert(ok==true); assert(ok==true);
// get RHO if exist // get RHO if exist
ok=uns->snapshot->getData("rho",&cnbody,&rho); ok=uns->snapshot->getData(select_c, "rho",&cnbody,&rho);
// get HSML if exist // get HSML if exist
ok=uns->snapshot->getData("hsml",&cnbody,&hsml); ok=uns->snapshot->getData(select_c,"hsml",&cnbody,&hsml);
std::cerr << "nbody=" << nbody << " time="<<time <<"\n"; std::cerr << "nbody=" << nbody << " time="<<time <<"\n";
...@@ -154,6 +154,8 @@ int main(int argc, char ** argv ) ...@@ -154,6 +154,8 @@ int main(int argc, char ** argv )
unsout->snapshot->setData("time",time); unsout->snapshot->setData("time",time);
unsout->snapshot->setData("pos" ,nbody,pos,false); unsout->snapshot->setData("pos" ,nbody,pos,false);
unsout->snapshot->setData("mass" ,nbody,mass,false); unsout->snapshot->setData("mass" ,nbody,mass,false);
if (vel)
unsout->snapshot->setData("vel" ,nbody,vel,false);
if (rho) if (rho)
unsout->snapshot->setData("rho" ,nbody,rho ,false); unsout->snapshot->setData("rho" ,nbody,rho ,false);
if (hsml) if (hsml)
......
...@@ -111,6 +111,21 @@ template <class T> void displayInfo(bool display,int maxlines, std::string comp, ...@@ -111,6 +111,21 @@ template <class T> void displayInfo(bool display,int maxlines, std::string comp,
if (ok && display) { if (ok && display) {
displayFormat(maxlines,"metal[1] = ",metal,1,nbody, 3); displayFormat(maxlines,"metal[1] = ",metal,1,nbody, 3);
} }
// RAMSES hydro
if (comp=="gas") {
int nvarh=0;
ok = uns->getData("nvarh",&nvarh);
std::cerr << "nvarh=" << nvarh << "\n";
for (int i=0; i<std::min(nvarh,20); i++) {
T * hydro;
stringstream s_int;
s_int << i;
ok = uns->snapshot->getData("hydro",s_int.str(),&nbody,&hydro);
if (ok && display) {
displayFormat(maxlines,"hydro["+s_int.str()+"] = ",hydro,1,nbody, 3);
}
}
}
//} //}
//if (comp == "stars") { //if (comp == "stars") {
T * age;//, * metal; T * age;//, * metal;
...@@ -204,6 +219,7 @@ template <class T> void start() ...@@ -204,6 +219,7 @@ template <class T> void start()
std::cout << "File name : "<<uns->snapshot->getFileName()<<"\n"; std::cout << "File name : "<<uns->snapshot->getFileName()<<"\n";
std::cout << "File type : "<<stype<<"\n"; std::cout << "File type : "<<stype<<"\n";
int nbody; T time; int nbody; T time;
//std::vector<T> vrect=uns->snapshot->getData("","/PartType1/Coordinates");
// get the input number of bodies according to the selection // get the input number of bodies according to the selection
uns->snapshot->getData("nsel",&nbody); uns->snapshot->getData("nsel",&nbody);
// get the simulation time // get the simulation time
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment