Mapped grids

compressible/interface.py

@@ -419,7 +419,7 @@ def riemann_cgf(idir, ng,
                 lambdastar_r = ustar + cstar_r
 
                 if (pstar > p_r):
-                    # the wave if a shock -- find the shock speed
+                    # the wave is a shock -- find the shock speed
                     sigma = (lambda_r + lambdastar_r) / 2.0
 
                     if (sigma > 0.0):
@@ -713,7 +713,7 @@ def riemann_prim(idir, ng,
                 lambdastar_r = ustar + cstar_r
 
                 if (pstar > p_r):
-                    # the wave if a shock -- find the shock speed
+                    # the wave is a shock -- find the shock speed
                     sigma = (lambda_r + lambdastar_r) / 2.0
 
                     if (sigma > 0.0):
@@ -801,7 +801,7 @@ def riemann_prim(idir, ng,
     return q_int
 
 
-@njit(cache=True)
+# @njit(cache=True)
 def riemann_hllc(idir, ng,
                  idens, ixmom, iymom, iener, irhoX, nspec,
                  lower_solid, upper_solid,
@@ -857,6 +857,9 @@ def riemann_hllc(idir, ng,
             # primitive variable states
             rho_l = U_l[i, j, idens]
 
+            if (rho_l <= 0):
+                print(f"-ve density {rho_l} at i,j = ({i},{j})")
+
             # un = normal velocity; ut = transverse velocity
             if (idir == 1):
                 un_l = U_l[i, j, ixmom] / rho_l

compressible/problems/acoustic_pulse.py

@@ -41,18 +41,19 @@ def init_data(myd, rp):
     ymin = rp.get_param("mesh.ymin")
     ymax = rp.get_param("mesh.ymax")
 
-    xctr = 0.5*(xmin + xmax)
-    yctr = 0.5*(ymin + ymax)
+    xctr = 0.5 * (xmin + xmax)
+    yctr = 0.5 * (ymin + ymax)
 
     dist = np.sqrt((myd.grid.x2d - xctr)**2 +
                    (myd.grid.y2d - yctr)**2)
 
     dens[:, :] = rho0
     idx = dist <= 0.5
-    dens[idx] = rho0 + drho0*np.exp(-16*dist[idx]**2) * np.cos(np.pi*dist[idx])**6
+    dens[idx] = rho0 + drho0 * \
+        np.exp(-16 * dist[idx]**2) * np.cos(np.pi * dist[idx])**6
 
-    p = (dens/rho0)**gamma
-    ener[:, :] = p/(gamma - 1)
+    p = (dens / rho0)**gamma
+    ener[:, :] = p / (gamma - 1)
 
 
 def finalize():

compressible/problems/advect.py

@@ -38,22 +38,23 @@ 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)
+    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))
+    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
-    xmom[:, :] = dens[:, :]*u
-    ymom[:, :] = dens[:, :]*v
+    xmom[:, :] = dens[:, :] * u
+    ymom[:, :] = dens[:, :] * v
 
     # pressure is constant
     p = 1.0
-    ener[:, :] = p/(gamma - 1.0) + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]
+    ener[:, :] = p / (gamma - 1.0) + 0.5 * (xmom[:, :]
+                                            ** 2 + ymom[:, :]**2) / dens[:, :]
 
 
 def finalize():

compressible/problems/test.py

@@ -9,7 +9,7 @@ def init_data(my_data, rp):
 
     # make sure that we are passed a valid patch object
     if not isinstance(my_data, patch.CellCenterData2d):
-        print("ERROR: patch invalid in sedov.py")
+        print("ERROR: patch invalid in test.py")
         print(my_data.__class__)
         sys.exit()
 

compressible/simulation.py

@@ -10,18 +10,18 @@
 import compressible.derives as derives
 import compressible.unsplit_fluxes as flx
 import mesh.boundary as bnd
-from simulation_null import NullSimulation, grid_setup, bc_setup
+from simulation_null import NullSimulation, NullVariables, grid_setup, bc_setup
 import util.plot_tools as plot_tools
 import particles.particles as particles
 
 
-class Variables(object):
+class Variables(NullVariables):
     """
     a container class for easy access to the different compressible
     variable by an integer key
     """
-    def __init__(self, myd):
-        self.nvar = len(myd.names)
+
+    def initialize(self, myd):
 
         # conserved variables -- we set these when we initialize for
         # they match the CellCenterData2d object
@@ -58,19 +58,19 @@ def cons_to_prim(U, gamma, ivars, myg):
     q = myg.scratch_array(nvar=ivars.nq)
 
     q[:, :, ivars.irho] = U[:, :, ivars.idens]
-    q[:, :, ivars.iu] = U[:, :, ivars.ixmom]/U[:, :, ivars.idens]
-    q[:, :, ivars.iv] = U[:, :, ivars.iymom]/U[:, :, ivars.idens]
+    q[:, :, ivars.iu] = U[:, :, ivars.ixmom] / U[:, :, ivars.idens]
+    q[:, :, ivars.iv] = U[:, :, ivars.iymom] / U[:, :, ivars.idens]
 
     e = (U[:, :, ivars.iener] -
-         0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 +
-                                  q[:, :, ivars.iv]**2))/q[:, :, ivars.irho]
+         0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 +
+                                      q[:, :, ivars.iv]**2)) / q[:, :, ivars.irho]
 
     q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e)
 
     if ivars.naux > 0:
-        for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux),
-                          range(ivars.irhox, ivars.irhox+ivars.naux)):
-            q[:, :, nq] = U[:, :, nu]/q[:, :, ivars.irho]
+        for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux),
+                          range(ivars.irhox, ivars.irhox + ivars.naux)):
+            q[:, :, nq] = U[:, :, nu] / q[:, :, ivars.irho]
 
     return q
 
@@ -81,18 +81,18 @@ def prim_to_cons(q, gamma, ivars, myg):
     U = myg.scratch_array(nvar=ivars.nvar)
 
     U[:, :, ivars.idens] = q[:, :, ivars.irho]
-    U[:, :, ivars.ixmom] = q[:, :, ivars.iu]*U[:, :, ivars.idens]
-    U[:, :, ivars.iymom] = q[:, :, ivars.iv]*U[:, :, ivars.idens]
+    U[:, :, ivars.ixmom] = q[:, :, ivars.iu] * U[:, :, ivars.idens]
+    U[:, :, ivars.iymom] = q[:, :, ivars.iv] * U[:, :, ivars.idens]
 
     rhoe = eos.rhoe(gamma, q[:, :, ivars.ip])
 
-    U[:, :, ivars.iener] = rhoe + 0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 +
-                                                           q[:, :, ivars.iv]**2)
+    U[:, :, ivars.iener] = rhoe + 0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 +
+                                                               q[:, :, ivars.iv]**2)
 
     if ivars.naux > 0:
-        for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux),
-                          range(ivars.irhox, ivars.irhox+ivars.naux)):
-            U[:, :, nu] = q[:, :, nq]*q[:, :, ivars.irho]
+        for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux),
+                          range(ivars.irhox, ivars.irhox + ivars.naux)):
+            U[:, :, nu] = q[:, :, nq] * q[:, :, ivars.irho]
 
     return U
 
@@ -113,7 +113,8 @@ def initialize(self, extra_vars=None, ng=4):
 
         # define solver specific boundary condition routines
         bnd.define_bc("hse", BC.user, is_solid=False)
-        bnd.define_bc("ramp", BC.user, is_solid=False)  # for double mach reflection problem
+        # for double mach reflection problem
+        bnd.define_bc("ramp", BC.user, is_solid=False)
 
         bc, bc_xodd, bc_yodd = bc_setup(self.rp)
 
@@ -182,10 +183,10 @@ def method_compute_timestep(self):
         u, v, cs = self.cc_data.get_var(["velocity", "soundspeed"])
 
         # the timestep is min(dx/(|u| + cs), dy/(|v| + cs))
-        xtmp = self.cc_data.grid.dx/(abs(u) + cs)
-        ytmp = self.cc_data.grid.dy/(abs(v) + cs)
+        xtmp = self.cc_data.grid.dx / (abs(u) + cs)
+        ytmp = self.cc_data.grid.dy / (abs(v) + cs)
 
-        self.dt = cfl*float(min(xtmp.min(), ytmp.min()))
+        self.dt = cfl * float(min(xtmp.min(), ytmp.min()))
 
     def evolve(self):
         """
@@ -211,19 +212,19 @@ def evolve(self):
         old_ymom = ymom.copy()
 
         # conservative update
-        dtdx = self.dt/myg.dx
-        dtdy = self.dt/myg.dy
+        dtdx = self.dt / myg.dx
+        dtdy = self.dt / myg.dy
 
         for n in range(self.ivars.nvar):
             var = self.cc_data.get_var_by_index(n)
 
             var.v()[:, :] += \
-                dtdx*(Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \
-                dtdy*(Flux_y.v(n=n) - Flux_y.jp(1, n=n))
+                dtdx * (Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \
+                dtdy * (Flux_y.v(n=n) - Flux_y.jp(1, n=n))
 
         # gravitational source terms
-        ymom[:, :] += 0.5*self.dt*(dens[:, :] + old_dens[:, :])*grav
-        ener[:, :] += 0.5*self.dt*(ymom[:, :] + old_ymom[:, :])*grav
+        ymom[:, :] += 0.5 * self.dt * (dens[:, :] + old_dens[:, :]) * grav
+        ener[:, :] += 0.5 * self.dt * (ymom[:, :] + old_ymom[:, :]) * grav
 
         if self.particles is not None:
             self.particles.update_particles(self.dt)
@@ -257,7 +258,7 @@ def dovis(self):
         u = q[:, :, ivars.iu]
         v = q[:, :, ivars.iv]
         p = q[:, :, ivars.ip]
-        e = eos.rhoe(gamma, p)/rho
+        e = eos.rhoe(gamma, p) / rho
 
         magvel = np.sqrt(u**2 + v**2)
 
@@ -297,7 +298,7 @@ def dovis(self):
 
             # plot particles
             ax.scatter(particle_positions[:, 0],
-                particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys")
+                       particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys")
             ax.set_xlim([myg.xmin, myg.xmax])
             ax.set_ylim([myg.ymin, myg.ymax])
 

compressible_mapped/__init__.py

@@ -0,0 +1,8 @@
+"""A method-of-lines compressible hydrodynamics solver.  Piecewise
+constant reconstruction is done in space and a Runge-Kutta time
+integration is used to advance the solutiion.
+
+"""
+__all__ = ["simulation"]
+
+from .simulation import Simulation

compressible_mapped/_defaults

@@ -0,0 +1,24 @@
+[driver]
+cfl = 0.8
+
+
+[eos]
+gamma = 1.4    ; pres = rho ener (gamma - 1)
+
+
+[compressible]
+use_flattening = 1        ; apply flattening at shocks (1)
+
+z0 = 0.75                 ; flattening z0 parameter
+z1 = 0.85                 ; flattening z1 parameter
+delta = 0.33              ; flattening delta parameter
+
+cvisc = 0.1               ; artifical viscosity coefficient
+
+limiter = 2               ; limiter (0 = none, 1 = 2nd order, 2 = 4th order)
+
+temporal_method = RK4     ; integration method (see mesh/integration.py)
+
+grav = 0.0                ; gravitational acceleration (in y-direction)
+
+riemann = HLLC            ; HLLC or CGF