add spherical geometry support in compressible (#204)

this extends the unsplit / CTU compressible solver to work in a 2D spherical geometry (r, theta) modeled with
phi axisymmetry.

Author: zhichen3

pyro/compressible/BC.py

@@ -53,7 +53,9 @@ def user(bc_name, bc_edge, variable, ccdata):
 
             # we will take the density to be constant, the velocity to
             # be outflow, and the pressure to be in HSE
-            if variable in ["density", "x-momentum", "y-momentum", "ymom_src", "E_src", "fuel", "ash"]:
+            if variable in ["density", "x-momentum", "y-momentum",
+                            "dens_src", "xmom_src", "ymom_src", "E_src",
+                            "fuel", "ash"]:
                 v = ccdata.get_var(variable)
                 j = myg.jlo-1
                 while j >= 0:
@@ -99,7 +101,9 @@ def user(bc_name, bc_edge, variable, ccdata):
 
             # we will take the density to be constant, the velocity to
             # be outflow, and the pressure to be in HSE
-            if variable in ["density", "x-momentum", "y-momentum", "ymom_src", "E_src", "fuel", "ash"]:
+            if variable in ["density", "x-momentum", "y-momentum",
+                            "dens_src", "xmom_src", "ymom_src", "E_src",
+                            "fuel", "ash"]:
                 v = ccdata.get_var(variable)
                 for j in range(myg.jhi+1, myg.jhi+myg.ng+1):
                     v[:, j] = v[:, myg.jhi]

pyro/compressible/interface.py

@@ -67,7 +67,7 @@ def states(idir, ng, dx, dt,
         Are we predicting to the edges in the x-direction (1) or y-direction (2)?
     ng : int
         The number of ghost cells
-    dx : float
+    dx : ndarray
         The cell spacing
     dt : float
         The timestep
@@ -133,13 +133,13 @@ def states(idir, ng, dx, dt,
                                  q[irho] / cs, 0.0, 0.5 / (cs * cs)]
                 lvec[1, :ns] = [1.0, 0.0,
                                  0.0, -1.0 / (cs * cs)]
-                lvec[2, :ns] = [0.0, 0.0,             1.0, 0.0]
+                lvec[2, :ns] = [0.0, 0.0, 1.0, 0.0]
                 lvec[3, :ns] = [0.0, 0.5 *
                                  q[irho] / cs,  0.0, 0.5 / (cs * cs)]
 
                 rvec[0, :ns] = [1.0, -cs / q[irho], 0.0, cs * cs]
-                rvec[1, :ns] = [1.0, 0.0,       0.0, 0.0]
-                rvec[2, :ns] = [0.0, 0.0,       1.0, 0.0]
+                rvec[1, :ns] = [1.0, 0.0, 0.0, 0.0]
+                rvec[2, :ns] = [0.0, 0.0, 1.0, 0.0]
                 rvec[3, :ns] = [1.0, cs / q[irho],  0.0, cs * cs]
 
                 # now the species -- they only have a 1 in their corresponding slot
@@ -174,29 +174,30 @@ def states(idir, ng, dx, dt,
             if idir == 1:
                 # this is one the right face of the current zone,
                 # so the fastest moving eigenvalue is e_val[3] = u + c
-                factor = 0.5 * (1.0 - dtdx * max(e_val[3], 0.0))
+                factor = 0.5 * (1.0 - dtdx[i, j] * max(e_val[3], 0.0))
                 q_l[i + 1, j, :] = q + factor * dq
 
                 # left face of the current zone, so the fastest moving
                 # eigenvalue is e_val[3] = u - c
-                factor = 0.5 * (1.0 + dtdx * min(e_val[0], 0.0))
+                factor = 0.5 * (1.0 + dtdx[i, j] * min(e_val[0], 0.0))
                 q_r[i,  j, :] = q - factor * dq
 
             else:
 
-                factor = 0.5 * (1.0 - dtdx * max(e_val[3], 0.0))
+                factor = 0.5 * (1.0 - dtdx[i, j] * max(e_val[3], 0.0))
                 q_l[i, j + 1, :] = q + factor * dq
 
-                factor = 0.5 * (1.0 + dtdx * min(e_val[0], 0.0))
+                factor = 0.5 * (1.0 + dtdx[i, j] * min(e_val[0], 0.0))
                 q_r[i, j, :] = q - factor * dq
 
             # compute the Vhat functions
             for m in range(nvar):
                 asum = np.dot(lvec[m, :], dq)
 
-                betal[m] = dtdx4 * (e_val[3] - e_val[m]) * \
+                # Should we change to max(e_val[3], 0.0) and min(e_val[0], 0.0)?
+                betal[m] = dtdx4[i, j] * (e_val[3] - e_val[m]) * \
                     (np.copysign(1.0, e_val[m]) + 1.0) * asum
-                betar[m] = dtdx4 * (e_val[0] - e_val[m]) * \
+                betar[m] = dtdx4[i, j] * (e_val[0] - e_val[m]) * \
                     (1.0 - np.copysign(1.0, e_val[m])) * asum
 
             # construct the states

pyro/compressible/problems/advect.py

@@ -38,16 +38,38 @@ def init_data(my_data, rp):
     ymin = rp.get_param("mesh.ymin")
     ymax = rp.get_param("mesh.ymax")
 
-    xctr = 0.5*(xmin + xmax)
-    yctr = 0.5*(ymin + ymax)
+    myg = my_data.grid
 
-    # this is identical to the advection/smooth problem
-    dens[:, :] = 1.0 + np.exp(-60.0*((my_data.grid.x2d-xctr)**2 +
-                                    (my_data.grid.y2d-yctr)**2))
+    if myg.coord_type == 0:
+        xctr = 0.5*(xmin + xmax)
+        yctr = 0.5*(ymin + ymax)
+
+        # this is identical to the advection/smooth problem
+        dens[:, :] = 1.0 + np.exp(-60.0*((my_data.grid.x2d-xctr)**2 +
+                                         (my_data.grid.y2d-yctr)**2))
+
+        # velocity is diagonal
+        u = 1.0
+        v = 1.0
+
+    else:
+        x = myg.scratch_array()
+        y = myg.scratch_array()
+
+        xctr = 0.5*(xmin + xmax) * np.sin((ymin + ymax) * 0.25)
+        yctr = 0.5*(xmin + xmax) * np.cos((ymin + ymax) * 0.25)
+
+        x[:, :] = myg.x2d.v(buf=myg.ng) * np.sin(myg.y2d.v(buf=myg.ng))
+        y[:, :] = myg.x2d.v(buf=myg.ng) * np.cos(myg.y2d.v(buf=myg.ng))
+
+        # this is identical to the advection/smooth problem
+        dens[:, :] = 1.0 + np.exp(-120.0*((x-xctr)**2 +
+                                         (y-yctr)**2))
+
+        # velocity in theta direction.
+        u = 0.0
+        v = 1.0
 
-    # velocity is diagonal
-    u = 1.0
-    v = 1.0
     xmom[:, :] = dens[:, :]*u
     ymom[:, :] = dens[:, :]*v
 

pyro/compressible/problems/inputs.advect.128

@@ -12,7 +12,7 @@ cvisc = 0.1
 
 
 [io]
-basename = advect_
+basename = advect_128_
 
 
 [eos]

pyro/compressible/problems/inputs.advect.256

@@ -12,7 +12,7 @@ cvisc = 0.1
 
 
 [io]
-basename = advect_
+basename = advect_256_
 
 
 [eos]

pyro/compressible/problems/inputs.advect.64

@@ -12,7 +12,7 @@ cvisc = 0.1
 
 
 [io]
-basename = advect_
+basename = advect_64_
 
 
 [eos]

pyro/compressible/problems/inputs.advect.spherical.128

@@ -0,0 +1,39 @@
+[driver]
+max_steps = 500
+tmax = 1.0
+
+init_tstep_factor = 1.0
+fix_dt = 0.0025
+
+
+[compressible]
+limiter = 0
+cvisc = 0.1
+riemann = CGF
+
+[io]
+basename = spherical_advect_128_
+
+
+[eos]
+gamma = 1.4
+
+
+[mesh]
+grid_type = SphericalPolar
+nx = 128
+ny = 128
+xmin = 1.0
+xmax = 2.0
+ymin = 0.523
+ymax = 2.617
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = outflow
+yrboundary = outflow
+
+
+[vis]
+dovis = 1

pyro/compressible/problems/inputs.advect.spherical.256

@@ -0,0 +1,39 @@
+[driver]
+max_steps = 1000
+tmax = 1.0
+
+init_tstep_factor = 1.0
+fix_dt = 0.00125
+
+
+[compressible]
+limiter = 0
+cvisc = 0.1
+riemann = CGF
+
+[io]
+basename = spherical_advect_256_
+
+
+[eos]
+gamma = 1.4
+
+
+[mesh]
+grid_type = SphericalPolar
+nx = 256
+ny = 256
+xmin = 1.0
+xmax = 2.0
+ymin = 0.523
+ymax = 2.617
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = outflow
+yrboundary = outflow
+
+
+[vis]
+dovis = 1

pyro/compressible/problems/inputs.advect.spherical.64

@@ -0,0 +1,39 @@
+[driver]
+max_steps = 500
+tmax = 1.0
+
+init_tstep_factor = 1.0
+fix_dt = 0.005
+
+
+[compressible]
+limiter = 0
+cvisc = 0.1
+riemann = CGF
+
+[io]
+basename = spherical_advect_64_
+
+
+[eos]
+gamma = 1.4
+
+
+[mesh]
+grid_type = SphericalPolar
+nx = 64
+ny = 64
+xmin = 1.0
+xmax = 2.0
+ymin = 0.523
+ymax = 2.617
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = outflow
+yrboundary = outflow
+
+
+[vis]
+dovis = 1

pyro/compressible/problems/inputs.sedov.spherical

@@ -0,0 +1,41 @@
+[driver]
+max_steps = 5000
+tmax = 0.1
+
+
+[compressible]
+limiter = 2
+cvisc = 0.1
+riemann = CGF
+
+[io]
+basename = sedov_spherical_
+dt_out = 0.0125
+
+
+[eos]
+gamma = 1.4
+
+
+[mesh]
+grid_type = SphericalPolar
+nx = 128
+ny = 128
+xmin = 0.1
+xmax = 1.0
+ymin = 0.785
+ymax = 2.355
+
+xlboundary = reflect-odd
+xrboundary = outflow
+
+ylboundary = outflow
+yrboundary = outflow
+
+
+[sedov]
+r_init = 0.13
+
+
+[vis]
+dovis = 1