Merge branch 'master' of ssh://github.com/zingale/pyro2

Author: zingale

CONTRIBUTING.md

@@ -0,0 +1,31 @@
+# Contributing
+
+Contributions are welcomed from anyone, including posting issues or
+submitting pull requests to the [pyro github](https://github.com/zingale/pyro).
+
+## Issues
+
+Creating an issue on github is a good way to request new features,
+file a bug report, or notify us of any difficulties that arise using
+pyro.
+
+To request support using pyro, please create an issue on the pyro
+github and the developers will be happy to assist you.
+
+If you are reporting a bug, please indicate any information necessary
+to reproduce the bug including your version of python.
+
+## Pull Requests
+
+*Any contributions that have the potential to change answers should be
+done via pull requests.* A pull request should be generated from your
+fork of pyro and target the `master` branch.
+
+The unit and regression tests will run automatically once the PR is
+submitted, and then one of the pyro developers will review the PR and
+if needed, suggest modifications prior to merging the PR.
+
+If there are a number of small commits making up the PR, we may wish
+to squash commits upon merge to have a clean history.  *Please ensure
+that your PR title and first post are descriptive, since these will be
+used for a squashed commit message.*

DEVELOPMENT

@@ -11,7 +11,7 @@ finite-volume framework.  The code is mainly written in python and is
 designed with simplicity in mind.  The algorithms are written to
 encourage experimentation and allow for self-learning of these code
 methods.
- 
+
 The latest version of pyro is always available at:
 
 https://github.com/zingale/pyro2
@@ -154,6 +154,8 @@ pyro provides the following solvers (all in 2-d):
     variable-coefficient Poisson equation (which inherits from the
     constant-coefficient solver).
 
+  - `swe`: a solver for the shallow water equations.
+
 
 ## Working with data
 
@@ -175,6 +177,12 @@ with their data.
 
   - `analysis/`
 
+      * `dam_compare.py`: this takes an output file from the
+        shallow water dam break problem and plots a slice through the domain
+        together with the analytic solution (calculated in the script).
+
+         usage: `./dam_compare.py file`
+
       * `gauss_diffusion_compare.py`: this is for the diffusion solver's
         Gaussian diffusion problem. It takes a sequence of output
         files as arguments, computes the angle-average, and the plots

README.md

@@ -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)

_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/_defaults

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

advection/problems/inputs.smooth

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

advection/problems/inputs.tophat

@@ -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/simulation.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/tests/test_advection.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)