Merge pull request #41 from harpolea/particles

Particles

Author: harpolea

_defaults

@@ -31,14 +31,15 @@ ymin = 0.0                ; domain minimum y-coordinate
 ymax = 1.0                ; domain maximum y-coordinate
 
 xlboundary = reflect     ; minimum x BC ('reflect', 'outflow', or 'periodic')
-xrboundary = reflect     ; maximum x BC ('reflect', 'outflow', or 'periodic') 
+xrboundary = reflect     ; maximum x BC ('reflect', 'outflow', or 'periodic')
 ylboundary = reflect     ; minimum y BC ('reflect', 'outflow', or 'periodic')
 yrboundary = reflect     ; maximum y BC ('reflect', 'outflow', or 'periodic')
 
 nx = 25                   ; number of zones in the x-direction
 ny = 25                   ; number of zones in the y-direction
 
 
-
-
-
+[particles]
+do_particles = 0          ; include particles? (1=yes, 0=no)
+n_particles = 100         ; number of particles
+particle_generator = random ; how do we generate particles? (random, grid)

advection/_defaults

@@ -9,5 +9,6 @@ v = 1.0      ; advective velocity in y-direction
 limiter = 2  ; limiter (0 = none, 1 = 2nd order, 2 = 4th order)
 
 
-
-
+[particles]
+do_particles = 0
+particle_generator = grid

advection/problems/inputs.smooth

@@ -34,3 +34,6 @@ v = 1.0
 
 limiter = 2
 
+[particles]
+do_particles = 1
+particle_generator = grid

advection/problems/inputs.tophat

@@ -32,4 +32,3 @@ u = 1.0
 v = 1.0
 
 limiter = 2
-

advection/simulation.py

@@ -5,6 +5,8 @@
 import advection.advective_fluxes as flx
 import mesh.patch as patch
 from simulation_null import NullSimulation, grid_setup, bc_setup
+import particles.particles as particles
+import util.plot_tools as plot_tools
 
 
 class Simulation(NullSimulation):
@@ -25,6 +27,11 @@ def initialize(self):
 
         self.cc_data = my_data
 
+        if self.rp.get_param("particles.do_particles") == 1:
+            n_particles = self.rp.get_param("particles.n_particles")
+            particle_generator = self.rp.get_param("particles.particle_generator")
+            self.particles = particles.Particles(self.cc_data, bc, n_particles, particle_generator)
+
         # now set the initial conditions for the problem
         problem = importlib.import_module("advection.problems.{}".format(self.problem_name))
         problem.init_data(self.cc_data, self.rp)
@@ -73,6 +80,16 @@ def evolve(self):
         dens.v()[:, :] = dens.v() + dtdx*(flux_x.v() - flux_x.ip(1)) + \
                                     dtdy*(flux_y.v() - flux_y.jp(1))
 
+        if self.particles is not None:
+            myg = self.cc_data.grid
+            u = self.rp.get_param("advection.u")
+            v = self.rp.get_param("advection.v")
+
+            u2d = myg.scratch_array() + u
+            v2d = myg.scratch_array() + v
+
+            self.particles.update_particles(self.dt, u2d, v2d)
+
         # increment the time
         self.cc_data.t += self.dt
         self.n += 1
@@ -87,16 +104,36 @@ def dovis(self):
 
         myg = self.cc_data.grid
 
-        plt.imshow(np.transpose(dens.v()),
+        _, axes, cbar_title = plot_tools.setup_axes(myg, 1)
+
+        # plot density
+        ax = axes[0]
+        img = ax.imshow(np.transpose(dens.v()),
                    interpolation="nearest", origin="lower",
                    extent=[myg.xmin, myg.xmax, myg.ymin, myg.ymax],
                    cmap=self.cm)
 
-        plt.xlabel("x")
-        plt.ylabel("y")
+        ax.set_xlabel("x")
+        ax.set_ylabel("y")
+
+        # needed for PDF rendering
+        cb = axes.cbar_axes[0].colorbar(img)
+        cb.solids.set_rasterized(True)
+        cb.solids.set_edgecolor("face")
+
         plt.title("density")
 
-        plt.colorbar()
+        if self.particles is not None:
+            particle_positions = self.particles.get_positions()
+
+            # dye particles
+            colors = self.particles.get_init_positions()[:, 0]
+
+            # plot particles
+            ax.scatter(particle_positions[:, 0],
+                particle_positions[:, 1], c=colors, cmap="Greys")
+            ax.set_xlim([myg.xmin, myg.xmax])
+            ax.set_ylim([myg.ymin, myg.ymax])
 
         plt.figtext(0.05, 0.0125, "t = {:10.5f}".format(self.cc_data.t))
 

advection/tests/test_advection.py

@@ -19,6 +19,7 @@ def setup_method(self):
 
         self.rp.params["mesh.nx"] = 8
         self.rp.params["mesh.ny"] = 8
+        self.rp.params["particles.do_particles"] = 0
 
         self.sim = sn.Simulation("advection", "test", self.rp)
         self.sim.initialize()

advection_fv4/simulation.py

@@ -6,6 +6,7 @@
 import mesh.array_indexer as ai
 import mesh.fv as fv
 from simulation_null import grid_setup, bc_setup
+import particles.particles as particles
 
 
 class Simulation(advection_rk.Simulation):
@@ -26,6 +27,11 @@ def initialize(self):
 
         self.cc_data = my_data
 
+        if self.rp.get_param("particles.do_particles") == 1:
+            n_particles = self.rp.get_param("particles.n_particles")
+            particle_generator = self.rp.get_param("particles.particle_generator")
+            self.particles = particles.Particles(self.cc_data, bc, n_particles, particle_generator)
+
         # now set the initial conditions for the problem
         problem = importlib.import_module("advection_fv4.problems.{}".format(self.problem_name))
         problem.init_data(self.cc_data, self.rp)

advection_rk/simulation.py

@@ -72,6 +72,15 @@ def evolve(self):
 
         rk.compute_final_update()
 
+        if self.particles is not None:
+            myg = self.cc_data.grid
+            u = self.rp.get_param("advection.u")
+            v = self.rp.get_param("advection.v")
+
+            u2d = myg.scratch_array() + u
+            v2d = myg.scratch_array() + v
+
+            self.particles.update_particles(self.dt, u2d, v2d)
         # increment the time
         myd.t += self.dt
         self.n += 1

advection_weno/simulation.py

@@ -72,6 +72,16 @@ def evolve(self):
 
         rk.compute_final_update()
 
+        if self.particles is not None:
+            myg = self.cc_data.grid
+            u = self.rp.get_param("advection.u")
+            v = self.rp.get_param("advection.v")
+
+            u2d = myg.scratch_array() + u
+            v2d = myg.scratch_array() + v
+
+            self.particles.update_particles(self.dt, u2d, v2d)
+
         # increment the time
         myd.t += self.dt
         self.n += 1

compressible/_defaults

@@ -21,3 +21,6 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = HLLC            ; HLLC or CGF
 
+[particles]
+do_particles = 0
+particle_generator = grid