Implemented midpoint method for particle update, added figure to docs

harpolea · harpolea · commit 9e85c6ba2a1c · 2018-07-06T15:28:56.000+01:00
diff --git a/compressible/simulation.py b/compressible/simulation.py
@@ -143,9 +143,7 @@ def initialize(self, extra_vars=None, ng=4):
         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)
+            self.particles = particles.Particles(self.cc_data, bc, self.rp)
 
         # some auxillary data that we'll need to fill GC in, but isn't
         # really part of the main solution
diff --git a/docs/source/particles_basics.rst b/docs/source/particles_basics.rst
@@ -127,3 +127,10 @@ x-position, we can plot them on the figure axis ``ax`` using the following code:
 
       ax.set_xlim([myg.xmin, myg.xmax])
       ax.set_ylim([myg.ymin, myg.ymax])
+
+Applying this to the Kelvin-Helmholtz problem with the ``compressible`` solver,
+we can produce a plot such as the one below, where the particles have been
+plotted on top of the fluid density.
+
+.. image:: particles_kh_compressible.png
+   :align: center
diff --git a/docs/source/particles_kh_compressible.png b/docs/source/particles_kh_compressible.png
diff --git a/particles/particles.py b/particles/particles.py
@@ -42,6 +42,48 @@ def update(self, u, v, dt):
         self.x += u * dt
         self.y += v * dt
 
+    def interpolate_velocity(self, myg, u, v):
+        """
+        Interpolate the x- and y-velocities defined on grid myg to the
+        particle's position.
+
+        Parameters
+        ----------
+        myg : Grid2d
+            grid which the velocities are defined on
+        u : ArrayIndexer
+            x-velocity
+        v : ArrayIndexer
+            y_velocity
+        """
+
+        # find what cell it lives in
+        x_idx = (self.x - myg.xmin) / myg.dx - 0.5
+        y_idx = (self.y - myg.ymin) / myg.dy - 0.5
+
+        x_frac = x_idx % 1
+        y_frac = y_idx % 1
+
+        # get the index of the particle's closest bottom left cell
+        # we'll add one as going to use buf'd grids -
+        # this will catch the cases where the particle is on the edges
+        # of the grid.
+        x_idx = int(x_idx) + 1
+        y_idx = int(y_idx) + 1
+
+        # interpolate velocity
+        u_vel = (1-x_frac)*(1-y_frac)*u.v(buf=1)[x_idx, y_idx] + \
+                x_frac*(1-y_frac)*u.v(buf=1)[x_idx+1, y_idx] + \
+                (1-x_frac)*y_frac*u.v(buf=1)[x_idx, y_idx+1] + \
+                x_frac*y_frac*u.v(buf=1)[x_idx+1, y_idx+1]
+
+        v_vel = (1-x_frac)*(1-y_frac)*v.v(buf=1)[x_idx, y_idx] + \
+                x_frac*(1-y_frac)*v.v(buf=1)[x_idx+1, y_idx] + \
+                (1-x_frac)*y_frac*v.v(buf=1)[x_idx, y_idx+1] + \
+                x_frac*y_frac*v.v(buf=1)[x_idx+1, y_idx+1]
+
+        return u_vel, v_vel
+
 
 class Particles(object):
     """
@@ -194,31 +236,20 @@ def update_particles(self, dt, u=None, v=None):
         elif v is None:
             v = self.sim_data.get_var("y-velocity")
 
-        for _, p in self.particles.items():
-            # find what cell it lives in
-            x_idx = (p.x - myg.xmin) / myg.dx - 0.5
-            y_idx = (p.y - myg.ymin) / myg.dy - 0.5
+        for k, p in self.particles.items():
+
+            # save old position and predict location at dt/2
+            initial_position = p.pos()
 
-            x_frac = x_idx % 1
-            y_frac = y_idx % 1
+            u_vel, v_vel = p.interpolate_velocity(myg, u, v)
 
-            # get the index of the particle's closest bottom left cell
-            # we'll add one as going to use buf'd grids -
-            # this will catch the cases where the particle is on the edges
-            # of the grid.
-            x_idx = int(x_idx) + 1
-            y_idx = int(y_idx) + 1
+            p.update(u_vel, v_vel, 0.5*dt)
 
-            # interpolate velocity
-            u_vel = (1-x_frac)*(1-y_frac)*u.v(buf=1)[x_idx, y_idx] + \
-                    x_frac*(1-y_frac)*u.v(buf=1)[x_idx+1, y_idx] + \
-                    (1-x_frac)*y_frac*u.v(buf=1)[x_idx, y_idx+1] + \
-                    x_frac*y_frac*u.v(buf=1)[x_idx+1, y_idx+1]
+            # find velocity at dt/2
+            u_vel, v_vel = p.interpolate_velocity(myg, u, v)
 
-            v_vel = (1-x_frac)*(1-y_frac)*v.v(buf=1)[x_idx, y_idx] + \
-                    x_frac*(1-y_frac)*v.v(buf=1)[x_idx+1, y_idx] + \
-                    (1-x_frac)*y_frac*v.v(buf=1)[x_idx, y_idx+1] + \
-                    x_frac*y_frac*v.v(buf=1)[x_idx+1, y_idx+1]
+            # update to final time using original position and vel at dt/2
+            p.x, p.y = initial_position
 
             p.update(u_vel, v_vel, dt)
 
@@ -234,6 +265,7 @@ def enforce_particle_boundaries(self):
         boundaries.
         """
         old_particles = self.particles.copy()
+        self.particles = dict()
 
         myg = self.sim_data.grid
 
@@ -248,7 +280,6 @@ def enforce_particle_boundaries(self):
             # -x boundary
             if p.x < myg.xmin:
                 if xlb in ["outflow", "neumann"]:
-                    del(self.particles[k])
                     continue
                 elif xlb == "periodic":
                     p.x = myg.xmax + p.x - myg.xmin
@@ -260,7 +291,6 @@ def enforce_particle_boundaries(self):
             # +x boundary
             if p.x > myg.xmax:
                 if xrb in ["outflow", "neumann"]:
-                    del(self.particles[k])
                     continue
                 elif xrb == "periodic":
                     p.x = myg.xmin + p.x - myg.xmax
@@ -272,7 +302,6 @@ def enforce_particle_boundaries(self):
             # -y boundary
             if p.y < myg.ymin:
                 if ylb in ["outflow", "neumann"]:
-                    del(self.particles[k])
                     continue
                 elif ylb == "periodic":
                     p.y = myg.ymax + p.y - myg.ymin
@@ -284,7 +313,6 @@ def enforce_particle_boundaries(self):
             # +y boundary
             if p.y > myg.ymax:
                 if yrb in ["outflow", "neumann"]:
-                    del(self.particles[k])
                     continue
                 elif yrb == "periodic":
                     p.y = myg.ymin + p.y - myg.ymax
@@ -293,6 +321,10 @@ def enforce_particle_boundaries(self):
                 else:
                     msg.fail("ERROR: yrb = %s invalid BC for particles" % (yrb))
 
+            self.particles[k] = p
+
+        self.n_particles = len(self.particles)
+
     def get_positions(self):
         """
         Return an array of current particle positions.
diff --git a/particles/tests/test_particles.py b/particles/tests/test_particles.py
@@ -41,6 +41,10 @@ def setup_test(n_particles=50, extra_rp_params=None):
 
     rp.params["mesh.nx"] = 8
     rp.params["mesh.ny"] = 8
+    rp.params["mesh.xmin"] = 0
+    rp.params["mesh.xmax"] = 1
+    rp.params["mesh.ymin"] = 0
+    rp.params["mesh.ymax"] = 1
     rp.params["particles.do_particles"] = 1
     n_particles = n_particles
 
@@ -134,11 +138,7 @@ def test_particles_advect():
     extra_rp_params = {"mesh.xlboundary": "periodic",
                        "mesh.xrboundary": "periodic",
                        "mesh.ylboundary": "periodic",
-                       "mesh.yrboundary": "periodic",
-                       "mesh.xmin": 0,
-                       "mesh.xmax": 1,
-                       "mesh.ymin": 0,
-                       "mesh.ymax": 1}
+                       "mesh.yrboundary": "periodic"}
 
     myd, bc, n_particles = setup_test(extra_rp_params=extra_rp_params)
 
@@ -164,11 +164,11 @@ def test_particles_advect():
     # now move constant speed to the right
     u[:, :] = 1
 
-    ps.update_particles(1, u, v)
+    ps.update_particles(0.1, u, v)
 
     positions = set([(p.x, p.y) for p in ps.particles.values()])
 
-    correct_positions = set([((x+1) % 1, y) for x in xs for y in xs])
+    correct_positions = set([((x+0.1) % 1, y) for x in xs for y in xs])
 
     assert positions == correct_positions, "sets are not the same"
 
@@ -177,23 +177,23 @@ def test_particles_advect():
     v[:, :] = 1
 
     ps = particles.Particles(myd, bc, n_particles, "grid")
-    ps.update_particles(1, u, v)
+    ps.update_particles(0.1, u, v)
 
     positions = set([(p.x, p.y) for p in ps.particles.values()])
 
-    correct_positions = set([(x, (y+1) % 1) for x in xs for y in xs])
+    correct_positions = set([(x, (y+0.1) % 1) for x in xs for y in xs])
 
     assert positions == correct_positions, "sets are not the same"
 
     # constant speed right + up
     u[:, :] = 1
 
     ps = particles.Particles(myd, bc, n_particles, "grid")
-    ps.update_particles(1, u, v)
+    ps.update_particles(0.1, u, v)
 
     positions = set([(p.x, p.y) for p in ps.particles.values()])
 
-    correct_positions = set([((x+1) % 1, (y+1) % 1) for x in xs for y in xs])
+    correct_positions = set([((x+0.1) % 1, (y+0.1) % 1) for x in xs for y in xs])
 
     assert positions == correct_positions, "sets are not the same"
 
@@ -206,17 +206,13 @@ def test_reflect_bcs():
     extra_rp_params = {"mesh.xlboundary": "reflect-even",
                        "mesh.xrboundary": "reflect-even",
                        "mesh.ylboundary": "reflect-even",
-                       "mesh.yrboundary": "reflect-even",
-                       "mesh.xmin": 0,
-                       "mesh.xmax": 1,
-                       "mesh.ymin": 0,
-                       "mesh.ymax": 1}
+                       "mesh.yrboundary": "reflect-even"}
 
     myd, bc, _ = setup_test(extra_rp_params=extra_rp_params)
 
     # create an array of particles at the edge of the domain.
-    init_particle_positions = [[0.5, 0.2], [0.5, 0.8],
-                               [0.2, 0.5], [0.8, 0.5]]
+    init_particle_positions = [[0.5, 0.03], [0.5, 0.96],
+                               [0.04, 0.5], [0.97, 0.5]]
 
     ps = particles.Particles(myd, bc, 4, "array", init_particle_positions)
 
@@ -232,15 +228,15 @@ def test_reflect_bcs():
     u[(x > y) & (x > (1-y))] = 1
     v[(x > y) & (x < (1-y))] = -1
 
-    ps.update_particles(0.3, u, v)
+    ps.update_particles(0.1, u, v)
 
     # extract positions by their initial positions so they're in the right order
     # (as particles are stored in a dictionary they may not be returned in the
     # same order as we first inserted them.)
     positions = [ps.particles[(x, y)].pos() for (x, y) in init_particle_positions]
 
-    correct_positions = [[0.5, 0.1], [0.5, 0.9],
-                         [0.1, 0.5], [0.9, 0.5]]
+    correct_positions = [[0.5, 0.07], [0.5, 0.94],
+                         [0.06, 0.5], [0.93, 0.5]]
 
     np.testing.assert_array_almost_equal(positions, correct_positions)
 
@@ -254,20 +250,14 @@ def test_outflow_bcs():
     extra_rp_params = {"mesh.xlboundary": "outflow",
                        "mesh.xrboundary": "outflow",
                        "mesh.ylboundary": "outflow",
-                       "mesh.yrboundary": "outflow",
-                       "mesh.xmin": 0,
-                       "mesh.xmax": 1,
-                       "mesh.ymin": 0,
-                       "mesh.ymax": 1,
-                       "mesh.nx": 100,
-                       "mesh.ny": 100}
+                       "mesh.yrboundary": "outflow"}
 
     myd, bc, _ = setup_test(extra_rp_params=extra_rp_params)
 
     # create an array of particles with some at the edge of the domain.
-    init_particle_positions = [[0.5, 0.2], [0.5, 0.8],
-                               [0.2, 0.5], [0.8, 0.5],
-                               [0.5, 0.4], [0.6, 0.5]]
+    init_particle_positions = [[0.5, 0.03], [0.5, 0.96],
+                               [0.04, 0.5], [0.97, 0.5],
+                               [0.5, 0.2], [0.8, 0.5]]
 
     ps = particles.Particles(myd, bc, 4, "array", init_particle_positions)
 
@@ -283,7 +273,7 @@ def test_outflow_bcs():
     u[(x > y) & (x > (1-y))] = 1
     v[(x > y) & (x < (1-y))] = -1
 
-    ps.update_particles(0.3, u, v)
+    ps.update_particles(0.1, u, v)
 
     assert len(ps.particles) == 2, "All but two of the particles should have flowed out of the domain"
 
