Merge pull request #50 from taataam/master

New solver for slotted-type advection problems

Author: harpolea

.coveragerc

@@ -0,0 +1,6 @@
+{
+  "compressible/tests/test_compressible.py": true,
+  "compressible/tests/test_eos.py": true,
+  "compressible_rk/tests/test_compressible_rk.py": true,
+  "swe/tests/test_swe.py": true
+}

.flake8

@@ -0,0 +1,3 @@
+[
+  "advection_nonuniform/tests/test_advection_nonuniform.py::TestSimulation::()::test_initializationst"
+]

.gitignore

@@ -0,0 +1,7 @@
+"""The pyro advection solver.  This implements a second-order,
+unsplit method for linear advection based on the Colella 1990 paper.
+
+"""
+
+__all__ = ['simulation']
+from .simulation import *

.pylintrc

@@ -0,0 +1,14 @@
+[driver]
+cfl       = 0.8           ; advective CFL number
+
+
+[advection]
+u = 1.0      ; advective velocity in x-direction
+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

.pytest_cache/v/cache/lastfailed

@@ -0,0 +1,126 @@
+import mesh.reconstruction as reconstruction
+import numpy as np
+
+
+def unsplit_fluxes(my_data, rp, dt, scalar_name):
+    """
+    Construct the fluxes through the interfaces for the linear advection
+    equation:
+
+    .. math::
+
+       a_t  + u a_x  + v a_y  = 0
+
+    We use a second-order (piecewise linear) unsplit Godunov method
+    (following Colella 1990).
+
+    In the pure advection case, there is no Riemann problem we need to
+    solve -- we just simply do upwinding.  So there is only one 'state'
+    at each interface, and the zone the information comes from depends
+    on the sign of the velocity.
+
+    Our convection is that the fluxes are going to be defined on the
+    left edge of the computational zones::
+
+        |             |             |             |
+        |             |             |             |
+       -+------+------+------+------+------+------+--
+        |     i-1     |      i      |     i+1     |
+
+                 a_l,i  a_r,i   a_l,i+1
+
+
+    a_r,i and a_l,i+1 are computed using the information in
+    zone i,j.
+
+    Parameters
+    ----------
+    my_data : CellCenterData2d object
+        The data object containing the grid and advective scalar that
+        we are advecting.
+    rp : RuntimeParameters object
+        The runtime parameters for the simulation
+    dt : float
+        The timestep we are advancing through.
+    scalar_name : str
+        The name of the variable contained in my_data that we are
+        advecting
+
+    Returns
+    -------
+    out : ndarray, ndarray
+        The fluxes on the x- and y-interfaces
+
+    """
+
+    myg = my_data.grid
+
+    a = my_data.get_var(scalar_name)
+
+    u = my_data.get_var("x-velocity")
+    v = my_data.get_var("y-velocity")
+
+    cx = u*dt/myg.dx
+    cy = v*dt/myg.dy
+
+    #--------------------------------------------------------------------------
+    # monotonized central differences
+    #--------------------------------------------------------------------------
+
+    limiter = rp.get_param("advection.limiter")
+
+    ldelta_ax = reconstruction.limit(a, myg, 1, limiter)
+    ldelta_ay = reconstruction.limit(a, myg, 2, limiter)
+
+    # upwind in x-direction
+    a_x = myg.scratch_array()
+    shift_x = my_data.get_var("x-shift").astype(int)
+
+    for index, vel in np.ndenumerate(u.v(buf=1)):
+        if vel < 0:
+            a_x.v(buf=1)[index] = a.ip(shift_x.v(buf=1)[index], buf=1)[index] - \
+                                  0.5*(1.0 + cx.v(buf=1)[index]) * \
+                                  ldelta_ax.ip(shift_x.v(buf=1)[index], buf=1)[index]
+        else:
+            a_x.v(buf=1)[index] = a.ip(shift_x.v(buf=1)[index], buf=1)[index] + \
+                                  0.5*(1.0 - cx.v(buf=1)[index]) * \
+                                  ldelta_ax.ip(shift_x.v(buf=1)[index], buf=1)[index]
+
+    # upwind in y-direction
+    a_y = myg.scratch_array()
+    shift_y = my_data.get_var("y-shift").astype(int)
+
+    for index, vel in np.ndenumerate(v.v(buf=1)):
+        if vel < 0:
+            a_y.v(buf=1)[index] = a.jp(shift_y.v(buf=1)[index], buf=1)[index] - \
+                                  0.5*(1.0 + cy.v(buf=1)[index]) * \
+                                  ldelta_ay.jp(shift_y.v(buf=1)[index], buf=1)[index]
+        else:
+            a_y.v(buf=1)[index] = a.jp(shift_y.v(buf=1)[index], buf=1)[index] + \
+                                  0.5*(1.0 - cy.v(buf=1)[index]) * \
+                                  ldelta_ay.jp(shift_y.v(buf=1)[index], buf=1)[index]
+
+    # compute the transverse flux differences.  The flux is just (u a)
+    # HOTF
+    F_xt = u*a_x
+    F_yt = v*a_y
+
+    F_x = myg.scratch_array()
+    F_y = myg.scratch_array()
+
+    # the zone where we grab the transverse flux derivative from
+    # depends on the sign of the advective velocity
+    dtdx2 = 0.5*dt/myg.dx
+    dtdy2 = 0.5*dt/myg.dy
+
+    for index, vel in np.ndenumerate(u.v(buf=1)):
+        F_x.v(buf=1)[index] = vel*(a_x.v(buf=1)[index] -
+                              dtdy2*(F_yt.ip_jp(shift_x.v(buf=1)[index], 1, buf=1)[index] -
+                                     F_yt.ip(shift_x.v(buf=1)[index], buf=1)[index]))
+
+    for index, vel in np.ndenumerate(v.v(buf=1)):
+        F_y.v(buf=1)[index] = vel*(a_y.v(buf=1)[index] -
+                              dtdx2*(F_xt.ip_jp(1, shift_y.v(buf=1)[index], buf=1)[index] -
+                                     F_xt.jp(shift_y.v(buf=1)[index], buf=1)[index]))
+
+    return F_x, F_y