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

Author: zingale

compressible_fv4/fluxes.py

@@ -40,6 +40,9 @@ def flux_cons(ivars, idir, gamma, q):
 
 def fluxes(myd, rp, ivars, solid, tc):
 
+    alpha = 0.3
+    beta = 0.3
+
     myg = myd.grid
 
     gamma = rp.get_param("eos.gamma")
@@ -156,4 +159,38 @@ def fluxes(myd, rp, ivars, solid, tc):
                                 2*F_avg.v(n=n, buf=1) +
                                 F_avg.ip(-1, n=n, buf=1))
 
+        # artificial viscosity McCorquodale & Colella Eq. 35, 36
+        # first find face-centered div
+        lam = myg.scratch_array()
+
+        if idir == 1:
+            lam.v(buf=1)[:, :] = (q_bar.v(buf=1, n=ivars.iu) -
+                                  q_bar.ip(-1, buf=1, n=ivars.iu))/myg.dx + \
+                                  0.25*(q_bar.jp(1, buf=1, n=ivars.iv) -
+                                        q_bar.jp(-1, buf=1, n=ivars.iv) +
+                                        q_bar.ip_jp(-1, 1, buf=1, n=ivars.iv) -
+                                        q_bar.ip_jp(-1, -1, buf=1, n=ivars.iv))/myg.dy
+        else:
+            lam.v(buf=1)[:, :] = (q_bar.v(buf=1, n=ivars.iv) -
+                                  q_bar.jp(-1, buf=1, n=ivars.iv))/myg.dy + \
+                                  0.25*(q_bar.ip(1, buf=1, n=ivars.iu) -
+                                        q_bar.ip(-1, buf=1, n=ivars.iu) +
+                                        q_bar.ip_jp(1, -1, buf=1, n=ivars.iu) -
+                                        q_bar.ip_jp(-1, -1, buf=1, n=ivars.iu))/myg.dx
+
+        test = myg.scratch_array()
+        test.v(buf=1)[:, :] = (myg.dx*lam.v(buf=1))**2 / \
+                                (beta * gamma * q_bar.v(buf=1, n=ivars.ip) /
+                                 q_bar.v(buf=1, n=ivars.irho))
+
+        nu = myg.dx * lam * np.minimum(test, 1.0)
+        nu[lam >= 0.0] = 0.0
+
+        if idir == 1:
+            for n in range(ivars.nvar):
+                F_x.v(buf=1, n=n)[:, :] += alpha * nu.v(buf=1) * (U_avg.v(buf=1, n=n) - U_avg.ip(-1, buf=1, n=n))
+        else:
+            for n in range(ivars.nvar):
+                F_y.v(buf=1, n=n)[:, :] += alpha * nu.v(buf=1) * (U_avg.v(buf=1, n=n) - U_avg.jp(-1, buf=1, n=n))
+
     return F_x, F_y

docs/source/compressible_compare.rst

@@ -1,5 +1,102 @@
 Compressible solver comparisons
 ===============================
 
+We run various problems run with the different compressible solvers in pyro (standard Riemann, Runge-Kutta, fourth order).
+
+
 Kelvin-Helmholtz
-----------------
+^^^^^^^^^^^^^^^^
+The McNally Kelvin-Helmholtz problem sets up a heavier fluid moving in the negative x-direction sandwiched between regions of lighter fluid moving in the positive x-direction.
+
+The image below shows the KH problem initialized with McNally's test. It ran on a 128 x 128 grid, with gamma = 1.7, and ran until t = 2.0. This is run with:
+
+.. code-block:: none
+
+  ./pyro.py compressible kh inputs.kh kh.vbulk=0
+  ./pyro.py compressible_rk kh inputs.kh kh.vbulk=0
+  ./pyro.py compressible_fv4 kh inputs.kh kh.vbulk=0
+
+.. image:: ./solver_comparisons/kh.png
+   :align: center
+	     
+
+We vary the velocity in the positive y-direction (vbulk) to see how effective the solvers are at preserving the initial shape.
+
+
+Sedov
+^^^^^
+The Sedov problem ran on a 128 x 128 grid, with gamma = 1.4, and until t = 0.1, which can be run as:
+
+.. code-block:: none
+  
+   ./pyro.py compressible sedov inputs.sedov
+   ./pyro.py compressible_rk sedov inputs.sedov
+   ./pyro.py compressible_fv4 sedov inputs.sedov
+
+.. image:: ./solver_comparisons/sedov.png
+   :align: center
+	
+.. image:: ./solver_comparisons/sedov_rk.png
+   :align: center
+	   
+.. image:: ./solver_comparisons/sedov_fv4.png
+   :align: center
+	   
+Quad
+^^^^
+The quad problem ran on a 256 x 256 grid until t = 0.8, which can be run as:
+
+.. code-block:: none
+
+   ./pyro.py compressible quad inputs.quad
+   ./pyro.py compressible_rk quad inputs.quad
+   ./pyro.py compressible_fv4 quad inputs.quad
+
+.. image:: ./solver_comparisons/quad.png
+   :align: center
+	
+.. image:: ./solver_comparisons/quad_rk.png
+   :align: center
+
+.. image:: ./solver_comparisons/quad_fv4.png
+   :align: center
+	      
+   
+Bubble
+^^^^^^
+The bubble problem ran on a 128 x 256 grid until t = 3.0, which can be run as:
+
+.. code-block:: none
+
+   ./pyro.py compressible bubble inputs.bubble
+   ./pyro.py compressible_rk bubble inputs.bubble
+   ./pyro.py compressible_fv4 bubble inputs.bubble
+
+.. image:: ./solver_comparisons/bubble.png
+   :align: center
+
+.. image:: ./solver_comparisons/bubble_rk.png
+   :align: center
+
+.. image:: ./solver_comparisons/bubble_fv4.png
+   :align: center
+
+
+Rayleigh-Taylor
+^^^^^^^^^^^^^^^
+The Rayleigh-Taylor problem ran on a 64 x 192 grid until t = 3.0, which can be run as:
+
+.. code-block:: none
+
+   ./pyro.py compressible rt inputs.rt
+   ./pyro.py compressible_rk rt inputs.rt
+   ./pyro.py compressible_fv4 rt inputs.rt
+
+.. image:: ./solver_comparisons/rt.png
+   :align: center
+
+.. image:: ./solver_comparisons/rt_rk.png
+   :align: center
+		      
+.. image:: ./solver_comparisons/rt_fv4.png
+   :align: center