Plausible results for incompressible

Author: harpolea

advection/simulation.py

@@ -89,6 +89,8 @@ def evolve(self):
             self.particles.update_particles(u2d, v2d, self.dt,
                 self.rp.get_param("advection.limiter"))
 
+            self.particles.enforce_particle_boundaries()
+
         # increment the time
         self.cc_data.t += self.dt
         self.n += 1

incompressible/_defaults

@@ -4,3 +4,7 @@ proj_type = 2             ; what are we projecting? 1 includes -Gp term in U*
 
 [driver]
 cfl = 0.8
+
+[particles]
+do_particles = 1
+particle_generator = grid

incompressible/simulation.py

@@ -12,6 +12,7 @@
 
 from simulation_null import NullSimulation, grid_setup, bc_setup
 import multigrid.MG as MG
+import particles.particles as particles
 
 
 class Simulation(NullSimulation):
@@ -43,6 +44,9 @@ def initialize(self):
 
         self.cc_data = my_data
 
+        if self.rp.get_param("particles.do_particles") == 1:
+            self.particles = particles.Particles(self.cc_data, bc, self.rp)
+
         # now set the initial conditions for the problem
         problem = importlib.import_module("incompressible.problems.{}".format(self.problem_name))
         problem.init_data(self.cc_data, self.rp)
@@ -388,6 +392,10 @@ def evolve(self):
         self.cc_data.fill_BC("x-velocity")
         self.cc_data.fill_BC("y-velocity")
 
+        if self.particles is not None:
+            self.particles.update_particles(u, v, self.dt, limiter)
+            self.particles.enforce_particle_boundaries()
+
         # increment the time
         if not self.in_preevolve:
             self.cc_data.t += self.dt
@@ -437,6 +445,15 @@ def dovis(self):
 
             plt.colorbar(img, ax=ax)
 
+        ax = axes.flat[0]
+        if self.particles is not None:
+            particle_positions = self.particles.get_positions()
+
+            # plot particles
+            ax.scatter(particle_positions[:, 0], particle_positions[:, 1], s=5)
+            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))
 
         plt.pause(0.001)

particles/particles.py

@@ -185,10 +185,10 @@ def update_particles(self, u, v, dt, limiter=0):
                 v_vel -= y_frac*(1.0 + cy)*ldelta_vy.v()[x_idx, y_idx]
             else:
                 v_vel += y_frac*(1.0 - cy)*ldelta_vy.v()[x_idx, y_idx]
-
-            # transverse velocity
-            u_vel += y_frac / myg.dy * ldelta_uy.v()[x_idx, y_idx]
-            v_vel += x_frac / myg.dx * ldelta_vx.v()[x_idx, y_idx]
+            #
+            # # transverse velocity
+            u_vel += y_frac * ldelta_uy.v()[x_idx, y_idx]
+            v_vel += x_frac * ldelta_vx.v()[x_idx, y_idx]
 
             p.advect(u_vel, v_vel, dt)
 

pyro.py

@@ -96,9 +96,6 @@ def doit(solver_name, problem_name, param_file,
         # evolve for a single timestep
         sim.evolve()
 
-        if sim.particles is not None:
-            sim.particles.enforce_particle_boundaries()
-
         if verbose > 0:
             print("%5d %10.5f %10.5f" % (sim.n, sim.cc_data.t, sim.dt))