New solver for slotted-type advection problems where velocity field is constant but nonuniform

Author: cheginit

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,130 @@
+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 = myg.scratch_array()
+    cy = myg.scratch_array()
+    
+    cx.v(buf=1)[:, :] = u.v(buf=1)*dt/myg.dx
+    cy.v(buf=1)[:, :] = v.v(buf=1)*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)
+    
+    # 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)):
+        # upwind
+        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]
+
+    # 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)):
+        # upwind
+        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