Particle post-processing code is now in slimPost.py

dgpy/scripts/particle.py

-
-"""Second-generation Louvain-la-Neuve Ice-ocean model - User interface
-
-    Contact: jonathan.lambrechts@uclouvain.be
-    webpage: www.climate.be/slim
-
-    SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.
-
-    particle is a post processing tools usable with results from slim (for simulations, see slim.py file) allowing for particle(s) tracking.
-
-    To generate the post processing results, run the script my_postpro_script.py:
-        - On a single core: rundgpy my_postpro_script.py
-
-    
-"""
-import dgpy
-import numpy as np
-import matplotlib.pyplot as plt
-import time
-import os
-from dgpy.scripts import slim_private
-
-
-class Particle_tracker:
-    """ predict the path of particles from the velocity field predicted by SLIM"""
-    def __init__(self, mesh_file_name, path_input, bathy, index_start, n_iter, n_period, dt, diffusivity = None) :
-        """keyword arguments:
-        
-        * mesh_file_name            
-            path to the mesh file (.msh format)
-        * path_input           
-            path to the output files provided by slim (sw2d folder ex: path_input/sw2d/eta/eta.idx)
-        * bathy         
-            bathymetric file [in meters, positive] (.msh, .idx or .nc format)
-        * diffusivity        
-            diffusivity of the model (.msh, .idx or .nc format) (default: 0)
-        * index_start
-            index of the first iteration 
-        * n_iter
-            number of iteration
-        * n_period
-            number of times the loop of hydrodynamic solution will be run 
-        * dt
-            time step 
-        """
-        self._mesh = mesh_file_name
-        self._path_input = path_input
-        self._index_start = index_start
-        self._n_iter = n_iter
-        self._nP = n_period
-        self._dt = dt
-        self._tot_iter = self._nP * self._n_iter
-        self._groups = dgpy.dgGroupCollection(self._mesh)
-        self._bath = dgpy.dgDofContainer(self._groups, 1)
-        slim_private._load(self._bath, bathy)
-        self._diffDof = dgpy.dgDofContainer(self._groups,1)
-        if diffusivity is None: 
-            self._alpha = dgpy.functionConstant(0)
-            self._diffDof.interpolate(self._alpha)
-        else: 
-            slim_private._load(self._diffDof, diffusivity)
-        self._particles = dgpy.particleArray()
-        self._parts = []
-          
-    def addParticleAtPoint(self, x, y, number_of_particles, time_step, locationId, status=0, particleId=-1):
-        """Seed particle(s) at location (x,y) and time step time_step
-        
-        keyword arguments:
-        
-        * x
-            x coordinates where the particle is seeded 
-        * y
-            y coordinates where the particle is seeded
-        * number_of_particles
-            number of particles seeded at this location
-        *time_step
-            index of the time step at which the particle is seeded
-        *status
-            status of the particle (dead or alive, default: 0)
-        *LocationId
-            Identification of the source
-        *particleId
-            Identification of the particle (default: each particle is identified depending on the source) 
-        """
-        self._parts.append([x, y, number_of_particles, time_step, status, locationId, particleId])
-        
-            
-    def loop(self, path_output):
-        """ Compute the particle(s) position at each time step and store it in a .pos file
-        
-        keyword arguments:
-        
-        * path_output
-            name of the .pos files which will contain the position of the particle(s)
-        """
-        self._path_output = path_output
-        if(not os.path.exists(path_output)):
-            try : os.mkdir(path_output)
-            except: pass
-      
-        self._hydro_sol_dof = dgpy.dgDofContainer(self._groups,3)
-        self._dof_eta = dgpy.dgDofContainer(self._groups,1)
-        self._dof_u = dgpy.dgDofContainer(self._groups,1)
-        self._dof_v = dgpy.dgDofContainer(self._groups,1)
-        self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%self._index_start,{0:0})
-        self._dof_eta.scatter()
-        self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:0})
-        self._dof_u.scatter()
-        self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:1})
-        self._dof_v.scatter()
-        self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
-        self._hydro_sol_dof.interpolate(self._hydrosolF)
-        
-        self._particleTracker = dgpy.dgParticleTracker2D(self._groups, self._particles, self._bath, self._diffDof, self._tot_iter, self._dt, self._path_output)
-        
-        startcpu = time.clock()
-        connectivityMatrix = dgpy.fullMatrixDouble(1,1)
-        for part in self._parts: 
-            if part[3]==0:
-                self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
-        self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( 0 ), 'pos', 0)
-        
-        for i in range(self._index_start, self._index_start+self._tot_iter+1):
-            for part in self._parts: 
-                if i>self._index_start and part[3]==(i-self._index_start):
-                    self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
-            print ('LPT iteration:', i, 'of',self._index_start+self._tot_iter+1,'CPU time:', time.clock()-startcpu)
-            _iter_read = self._index_start +(i-self._index_start)%self._nP
-            self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%_iter_read,{0:0})
-            self._dof_eta.scatter()
-            self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:0})
-            self._dof_u.scatter()
-            self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:1})
-            self._dof_v.scatter()
-            self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
-            self._hydro_sol_dof.interpolate(self._hydrosolF)
-            
-            self._particleTracker.particleMove(self._hydro_sol_dof, connectivityMatrix, i)
-            self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( i ), 'pos', 0)
-                

dgpy/scripts/slimPost.py

@@ -16,6 +16,8 @@ import dgpy
 import numpy as np
 import matplotlib.pyplot as plt
 import time
+import os
+from dgpy.scripts import slim_private
 
 class Fourier_series:
     """ compute the fourier series"""
@@ -270,4 +272,122 @@ class Fourier_series:
               PhaseMat.set(iNode, 0, Phase[k])
               k=k+1
         AmpliDof.exportMsh(range_file)
-        PhaseDof.exportMsh(phase_file)
\ No newline at end of file
+        PhaseDof.exportMsh(phase_file)
+        
+class Particle_tracker:
+    """ predict the path of particles from the velocity field predicted by SLIM"""
+    def __init__(self, mesh_file_name, path_input, bathy, index_start, n_iter, n_period, dt, diffusivity = None) :
+        """keyword arguments:
+        
+        * mesh_file_name            
+            path to the mesh file (.msh format)
+        * path_input           
+            path to the output files provided by slim (sw2d folder ex: path_input/sw2d/eta/eta.idx)
+        * bathy         
+            bathymetric file [in meters, positive] (.msh, .idx or .nc format)
+        * diffusivity        
+            diffusivity of the model (.msh, .idx or .nc format) (default: 0)
+        * index_start
+            index of the first iteration 
+        * n_iter
+            number of iteration
+        * n_period
+            number of times the loop of hydrodynamic solution will be run 
+        * dt
+            time step 
+        """
+        self._mesh = mesh_file_name
+        self._path_input = path_input
+        self._index_start = index_start
+        self._n_iter = n_iter
+        self._nP = n_period
+        self._dt = dt
+        self._tot_iter = self._nP * self._n_iter
+        self._groups = dgpy.dgGroupCollection(self._mesh)
+        self._bath = dgpy.dgDofContainer(self._groups, 1)
+        slim_private._load(self._bath, bathy)
+        self._diffDof = dgpy.dgDofContainer(self._groups,1)
+        if diffusivity is None: 
+            self._alpha = dgpy.functionConstant(0)
+            self._diffDof.interpolate(self._alpha)
+        else: 
+            slim_private._load(self._diffDof, diffusivity)
+        self._particles = dgpy.particleArray()
+        self._parts = []
+          
+    def addParticleAtPoint(self, x, y, number_of_particles, time_step, locationId, status=0, particleId=-1):
+        """Seed particle(s) at location (x,y) and time step time_step
+        
+        keyword arguments:
+        
+        * x
+            x coordinates where the particle is seeded 
+        * y
+            y coordinates where the particle is seeded
+        * number_of_particles
+            number of particles seeded at this location
+        *time_step
+            index of the time step at which the particle is seeded
+        *status
+            status of the particle (dead or alive, default: 0)
+        *LocationId
+            Identification of the source
+        *particleId
+            Identification of the particle (default: each particle is identified depending on the source) 
+        """
+        self._parts.append([x, y, number_of_particles, time_step, status, locationId, particleId])
+        
+            
+    def loop(self, path_output):
+        """ Compute the particle(s) position at each time step and store it in a .pos file
+        
+        keyword arguments:
+        
+        * path_output
+            name of the .pos files which will contain the position of the particle(s)
+        """
+        self._path_output = path_output
+        if(not os.path.exists(path_output)):
+            try : os.mkdir(path_output)
+            except: pass
+      
+        self._hydro_sol_dof = dgpy.dgDofContainer(self._groups,3)
+        self._dof_eta = dgpy.dgDofContainer(self._groups,1)
+        self._dof_u = dgpy.dgDofContainer(self._groups,1)
+        self._dof_v = dgpy.dgDofContainer(self._groups,1)
+        self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%self._index_start,{0:0})
+        self._dof_eta.scatter()
+        self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:0})
+        self._dof_u.scatter()
+        self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:1})
+        self._dof_v.scatter()
+        self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
+        self._hydro_sol_dof.interpolate(self._hydrosolF)
+        
+        self._particleTracker = dgpy.dgParticleTracker2D(self._groups, self._particles, self._bath, self._diffDof, self._tot_iter, self._dt, self._path_output)
+        
+        startcpu = time.clock()
+        connectivityMatrix = dgpy.fullMatrixDouble(1,1)
+        for part in self._parts: 
+            if part[3]==0:
+                self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
+        self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( 0 ), 'pos', 0)
+        
+        for i in range(self._index_start, self._index_start+self._tot_iter+1):
+            for part in self._parts: 
+                if i>self._index_start and part[3]==(i-self._index_start):
+                    self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
+            print ('LPT iteration:', i, 'of',self._index_start+self._tot_iter+1,'CPU time:', time.clock()-startcpu)
+            _iter_read = self._index_start +(i-self._index_start)%self._nP
+            self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%_iter_read,{0:0})
+            self._dof_eta.scatter()
+            self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:0})
+            self._dof_u.scatter()
+            self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:1})
+            self._dof_v.scatter()
+            self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
+            self._hydro_sol_dof.interpolate(self._hydrosolF)
+            
+            self._particleTracker.particleMove(self._hydro_sol_dof, connectivityMatrix, i)
+            self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( i ), 'pos', 0)
+