New solver for slotted-type advection problems

.pytest_cache/v/cache/lastfailed

@@ -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
+}

.pytest_cache/v/cache/nodeids

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

advection_nonuniform/__init__.py

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

advection_nonuniform/_defaults

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

advection_nonuniform/advective_fluxes.py

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

advection_nonuniform/problems/__init__.py

@@ -0,0 +1 @@
+__all__ = ['smooth', 'tophat', 'test', 'slotted']

advection_nonuniform/problems/_slotted.defaults

@@ -0,0 +1,3 @@
+[slotted]
+omega = 0.5     ;angular velocity
+offset = 0.25   ;offset of the slot center from domain center

advection_nonuniform/problems/_smooth.defaults

@@ -0,0 +1,2 @@
+[smooth]
+

advection_nonuniform/problems/_tophat.defaults

@@ -0,0 +1,2 @@
+[tophat]
+

advection_nonuniform/problems/inputs.slotted

@@ -0,0 +1,39 @@
+# simple inputs files for the unsplit CTU hydro scheme
+
+[driver]
+max_steps = 500
+tmax = 12.5664
+
+max_dt_change = 1.e33
+init_tstep_factor = 1.0
+
+
+[driver]
+cfl = 0.8
+
+[io]
+basename = slotted_
+dt_out = 0.2
+
+[mesh]
+nx = 64
+ny = 64
+xmax = 1.0
+ymax = 1.0
+
+xlboundary = periodic
+xrboundary = periodic
+
+ylboundary = periodic
+yrboundary = periodic
+
+
+[advection]
+u = 1.0
+v = 1.0
+
+limiter = 2
+
+[particles]
+do_particles = 0
+particle_generator = grid

advection_nonuniform/problems/inputs.smooth

@@ -0,0 +1,39 @@
+# simple inputs files for the unsplit CTU hydro scheme
+
+[driver]
+max_steps = 500
+tmax = 1.0
+
+max_dt_change = 1.e33
+init_tstep_factor = 1.0
+
+
+[driver]
+cfl = 0.8
+
+[io]
+basename = smooth_
+dt_out = 0.2
+
+[mesh]
+nx = 32
+ny = 32
+xmax = 1.0
+ymax = 1.0
+
+xlboundary = periodic
+xrboundary = periodic
+
+ylboundary = periodic
+yrboundary = periodic
+
+
+[advection]
+u = 1.0
+v = 1.0
+
+limiter = 2
+
+[particles]
+do_particles = 1
+particle_generator = grid

advection_nonuniform/problems/inputs.tophat

@@ -0,0 +1,34 @@
+# simple inputs files for the unsplit CTU hydro scheme
+
+[driver]
+max_steps = 500
+tmax = 1.0
+
+max_dt_change = 1.e33
+init_tstep_factor = 1.0
+
+cfl = 0.8
+
+
+[io]
+basename = tophat_
+n_out = 10
+
+[mesh]
+nx = 32
+ny = 32
+xmax = 1.0
+ymax = 1.0
+
+xlboundary = periodic
+xrboundary = periodic
+
+ylboundary = periodic
+yrboundary = periodic
+
+
+[advection]
+u = 1.0
+v = 1.0
+
+limiter = 2

advection_nonuniform/problems/slotted.py

@@ -0,0 +1,55 @@
+from __future__ import print_function
+import numpy as np
+
+import mesh.patch as patch
+from util import msg
+
+
+def init_data(my_data, rp):
+    """ initialize the slotted advection problem """
+    msg.bold("initializing the slotted advection problem...")
+
+    # make sure that we are passed a valid patch object
+    if not isinstance(my_data, patch.CellCenterData2d):
+        print(my_data.__class__)
+        msg.fail("ERROR: patch invalid in slotted.py")
+
+    offset = rp.get_param("slotted.offset")
+    omega = rp.get_param("slotted.omega")
+
+    myg = my_data.grid
+
+    xctr_dens = 0.5*(myg.xmin + myg.xmax)
+    yctr_dens = 0.5*(myg.ymin + myg.ymax) + offset
+
+    # setting initial condition for density
+    dens = my_data.get_var("density")
+    dens[:, :] = 0.0
+
+    R = 0.15
+    slot_width = 0.05
+
+    inside = (myg.x2d - xctr_dens)**2 + (myg.y2d - yctr_dens)**2 < R**2
+
+    slot_x = np.logical_and(myg.x2d > (xctr_dens - slot_width*0.5),
+                            myg.x2d < (xctr_dens + slot_width*0.5))
+    slot_y = np.logical_and(myg.y2d > (yctr_dens - R),
+                            myg.y2d < (yctr_dens))
+    slot = np.logical_and(slot_x, slot_y)
+
+    dens[inside] = 1.0
+    dens[slot] = 0.0
+
+    # setting initial condition for velocity
+    u = my_data.get_var("x-velocity")
+    v = my_data.get_var("y-velocity")
+
+    u[:, :] = omega*(myg.y2d - xctr_dens)
+    v[:, :] = -omega*(myg.x2d - (yctr_dens-offset))
+
+    print("extrema: ", np.amax(u), np.amin(u))
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """
+    pass

advection_nonuniform/problems/smooth.py

@@ -0,0 +1,41 @@
+from __future__ import print_function
+
+import mesh.patch as patch
+import numpy as np
+from util import msg
+
+
+def init_data(my_data, rp):
+    """ initialize the smooth advection problem """
+
+    msg.bold("initializing the smooth advection problem...")
+
+    # make sure that we are passed a valid patch object
+    if not isinstance(my_data, patch.CellCenterData2d):
+        print(my_data.__class__)
+        msg.fail("ERROR: patch invalid in slotted.py")
+
+    dens = my_data.get_var("density")
+
+    xmin = my_data.grid.xmin
+    xmax = my_data.grid.xmax
+
+    ymin = my_data.grid.ymin
+    ymax = my_data.grid.ymax
+
+    xctr = 0.5*(xmin + xmax)
+    yctr = 0.5*(ymin + ymax)
+
+    dens[:, :] = 1.0 + np.exp(-60.0*((my_data.grid.x2d-xctr)**2 +
+                                        (my_data.grid.y2d-yctr)**2))
+
+    u = my_data.get_var("x-velocity")
+    v = my_data.get_var("y-velocity")
+
+    u[:, :] = rp.get_param("advection.u")
+    v[:, :] = rp.get_param("advection.v")
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """
+    pass