crectify.py 11.2 KB
Newer Older
1
#!/usr/bin/python
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
2
3
4
from __future__ import print_function
from uns_simu import *
from csnapshot import *
5
import sys
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
6

7
8
9
10
try:
    import py_unstools # rectify swig
except ImportError:    
    print("WARNING !!!, failed to import module [py_unstools]",file=sys.stderr)
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
11
12
13
14

from multiprocessing import Lock
import time
import os
15
import glob
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
16
17
18
19
20
21

#
# ----
#
class CRectify:
    """
22
    Rectify UNS snapshots
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
23
    """
24
    __analysis=None
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
25
26
27
    __use_rho=False
    __dmin=0.
    __dmax=100.
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
28
29
30
31
    #
    # ----
    #
    # constructor
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
32
    def __init__(self,analysis=None,use_rho=False,dmin=0.,dmax=100.,verbose=False,verbose_debug=False):
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
33
34
35
36
37
38
39
        """
        Constructor of CRectify class

        - analysis : is a class object instantiate from CUnsAnalysis class
        """
        self.__vdebug=verbose_debug
        self.__verbose=verbose
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
40
41
42
        self.__use_rho=use_rho
        self.__dmin=dmin
        self.__dmax=dmax
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

        if analysis is not None:
            self.__analysis=analysis
            self.__smartAnalysisInit()

    #
    # __smartAnalysisInit
    #
    def __smartAnalysisInit(self):
        """
        start some initialisations
        """

        data=self.__analysis

        if not hasattr(data,'rectify_select'):
            print("\n\nYou must set a fied [rectify_select] in your <data> object, aborting...\n\n")
            sys.exit()

        data.rectify_select=data.rectify_select.replace(" ","") # remove blank

LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
64
        ### build CRECT Dir
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
65
66
67
68
69
70
71
72
        if hasattr(data,'rectify_dir'):
            self.__rectify_dir=data.rectify_dir
            #
        else: # default simdir simulation
            self.__rectify_dir=data.sim_info['dir']+"/ANALYSIS/rectify"

        print("RECTIFY DIR = ", self.__rectify_dir, data.sim_info['name'])
        self.simname=data.sim_info['name']
73
        #self.__com_file_base=self.__rectify_dir
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

        # lock process
        data.lock[data.lock_id].acquire()
        # build directory
        if (not os.path.isdir(self.__rectify_dir)) :
            try:
                print("Core ID ",data.core_id," create directory [%s]\n"%(self.__rectify_dir))
                os.makedirs(self.__rectify_dir)

            except OSError:
                print("Unable to create directory [%s]\n"%(self.__rectify_dir))
                data.lock[data.lock_id].release()
                sys.exit()

        # release process
        data.lock[data.lock_id].release()

        # check components vs cod
        ### COD Dir name
        if hasattr(data,'cod_dir'):
            pass
            #
        else: # default simdir simulation
97
98
            data.cod_dir=data.sim_info['dir']+"/ANALYSIS/cod"
            #data.cod_dir=data.sim_info['dir']+"/ANALYSIS"+self.__COD_DIR_NAME
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
        
        self.__cod_base_dir=data.cod_dir

        self.comp_cod=[] # list to store combo : component/cod to apply rectify

        ### 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
        self.__new_select=""
        for colon_s in data.rectify_select.split(":"): # colon separate analysis

            comp_L,cod_R=colon_s.split("#")    # hashtag separate comp_Left from cod_Right
            cod_R=cod_R.split("@")
            if len(cod_R)>1:
                rcut=float(cod_R[1])
                codf=cod_R[0]
            else:
115
                rcut=50.
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
                codf=cod_R[0]

            # check components exist?
            ok=True
            for comma_s in comp_L.split(","):
                xx=data.list_components.find(comma_s)  # find if selection exist
                if xx==-1:
                    ok=False
                    print("Rectify#Warning : In combo [%s], component <%s> from select <%s> does not exist...\n"\
                          %(colon_s,comma_s,comp_L))
                    break
            # components are ok, check cod file
            if ok: 
                cod_file=self.__cod_base_dir+"/"+self.simname+"."+codf+".cod"
                if (os.path.exists(cod_file)):
                    ok=True
                else:
                    ok=False
                    print("Rectify#Warning : In combo [%s], COD <%s> does not exist...\n"\
                          %(colon_s,cod_file))
            # both comp_L and cod_r exist, build combo
            if ok: 
                self.comp_cod.append([comp_L,codf,cod_file,rcut])
                

        if self.__vdebug:
            for comp_cod in self.comp_cod:
                print("Rectify combo: <%s> with [%s]\n"%(comp_cod[0],comp_cod[2]))

    # 
    # smartAnalysis
    #
    def smartAnalysis(self,analysis=None):
        """
150
        Main core function to compute RECTIFY on current snapshot stored in data_analysis
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
        """
        if analysis is None:
            data=self.__analysis
        else:
            data=analysis

        uns_snap=data.uns_snap # link to UNS object
        ok,time=uns_snap.getData("time")

        c=py_unstools.CRectify()

        # loop on combo comp/cod
        for comp_cod in self.comp_cod:
            select   = comp_cod[0]
            cod_sel  = comp_cod[1]
            cod_file = comp_cod[2]
            rcut     = comp_cod[3]
168
            print("select %s cod_sel %s cod_file %s rcut %s"%(select,cod_sel,cod_file,rcut))
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
169
170
            # pre computed rectify file
            pre_rect_file=self.__rectify_dir+"/"+self.simname+"."+select+"-"+cod_sel+".ev"
171
172
            print("rect_file =",pre_rect_file)

LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
173
174
175
            # lock file
            #data.lock[data.lock_id].acquire()
            exist_time,e_vectors=self.__checkTimeInPreRect(time,pre_rect_file)
176
            print("exist_time =",exist_time)
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
177
178
179
180
181
182
            #data.lock[data.lock_id].release()

            if not exist_time :
                ok1,pos  = uns_snap.getData(select,"pos" )
                ok2,vel  = uns_snap.getData(select,"vel" )
                ok3,mass = uns_snap.getData(select,"mass")
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
183
                ok4,rho = uns_snap.getData(select,"rho")
184
                #print("ok1 ok2 ok3 ",ok1,ok2,ok3)
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
185
186
187
                if pos.size>1 and vel.size>1 and mass.size>1 :
                    if self.__vdebug:
                        print("Rectifing time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file))
188
                    try:
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
189
190
                        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)
191
192
193
194
195
196
197
                        print("Time [%f] eigen:\n"%(time),eigen_v)
                        data.lock[data.lock_id].acquire()
                        self.__saveEigenVectors(pre_rect_file,eigen_v)
                        data.lock[data.lock_id].release()
                    except:
                        print("UNABLE TO COMPUTE EIGENVECTORS for time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file))
                        
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
                

        
    # 
    # checkTimeInPreRect
    #
    def __checkTimeInPreRect(self,time,eigen_file):
        """
        check in pre-computed Eigen Vectors stored in file if time has been already computed
        
        *IN*
        time    : current time
        egein_file : file with eigen vectors values

        *OUT*
        boolean : True if time exist, False otherwise
        numpy   : numppy array[16] => time[1] codx[3] codv[3] eigen_vectors[9]
        array
        """

        if os.path.isfile(eigen_file):
            try:
                f=open(eigen_file,"r")
                while True:
                    tce=map(np.float64,(f.readline().split()))
                    if (len(tce)>0):
                        atime=tce[0]
                        if (atime-0.001)<time and (atime+0.001)>time:
                            f.close()
                            return True,tce
                    else:
                        f.close()
                        return False,np.array([])

            except EOFError:
                pass

        return False,np.array([])

    # 
    # saveEigenVectors
    #
    def __saveEigenVectors(self,pre_rect_file,eigen_v):
        """
        """

        try:
            print("SAVING to pre_rect_file[%s]\n"%(pre_rect_file))
            f=open(pre_rect_file,"a")
            eigen_v.tofile(f,sep=" ",format="%e")
            f.write("\n")
            f.close()
        except:
            print("\n\nUnable to save pre_rect_file[%s] ... aborting...\n"%(pre_rect_file))
            sys.exit()
        
    # 
    # buildRectFile
    #
    def buildRectFile(self,simname=None,ev_in=None,rect_out=None):
        """
        Build rect file from eigen vectors file
        """

        if simname is not None:
            self.simname = simname
264
            self.__sql3 = UnsSimu(simname,verbose=self.__verbose)
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
265
266
267
268
269
270
271
272
273
            self.__r = self.__sql3.getInfo() # return None if does not exist

            if self.__vdebug:
                self.__sql3.printInfo(simname)
            self.__slist = self.__sql3.getSnapshotList()
            if (self.__slist is None):
                message="In CLASS <"+self.__class__.__name__+"> :  Unknown simulation ["+simname+\
                        "] in UNS database..."
                raise Exception(message)
274
            if ev_in is None:
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
275
              simdir=self.__r['dir']+"/ANALYSIS/rectify2/"
276
277
278
279
280
              for ev_in in glob.glob(simdir+"*.ev"):
                rect_out=os.path.basename(ev_in)
                rect_out=simdir+rect_out.split(".ev")[0]+".rect"
                print(">> %s [%s]\n"%(ev_in,rect_out))
                self.__convertEv2Rect(ev_in,rect_out)
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
281
282
283
284
285
286
287
288
289
290
291
292
293
294
        
        else :
            self.__convertEv2Rect(ev_in,rect_out)
        

    # 
    # __convertEv2Rect
    #
    def __convertEv2Rect(self,ev_in=None,rect_out=None):
        if ev_in is not None and os.path.isfile(ev_in):
            f=None
            try:
                if rect_out is not None:
                    f=open(rect_out,"w")
295
            except (IOError,EOFError):
LAMBERT Jean-charles's avatar
LAMBERT Jean-charles committed
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
                print("Unable to open file [%s] for writing, aborting....\n"%(rect_out))
                sys.exit()
            a=np.loadtxt(ev_in) # load file
            ii=a[:,0].argsort()
            ref0=np.array([1.,0.,0.])
            ref1=np.array([0.,1.,0.])
            ref2=np.array([0.,0.,1.])
            for id in ii:
                time=a[id,0]
                cx=a[id,1:4]
                cv=a[id,4:7]
                ev0=a[id,7:10]
                ev1=a[id,10:13]
                ev2=a[id,13:16]
                tmp=np.dot(ref0,ev0)
                if tmp<0.0:
                    ev0 = ev0*-1.
                tmp=np.dot(ref2,ev2)
                if tmp<0.0:
                    ev2 = ev2*-1.
                ev1=np.cross(ev2,ev0)
                ref0=ev0
                ref1=ev1
                ref2=ev2
                if f is not None:
                    f.write("%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n"\
                            %(a[id,0],cx[0],cx[1],cx[2],cv[0],cv[1],cv[2],ev0[0],\
                              ev0[1],ev0[2],ev1[0],ev1[1],ev1[2],ev2[0],ev2[1],ev2[2]))
                else:
                    print("%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e"\
                            %(a[id,0],cx[0],cx[1],cx[2],cv[0],cv[1],cv[2],ev0[0],\
                              ev0[1],ev0[2],ev1[0],ev1[1],ev1[2],ev2[0],ev2[1],ev2[2]))