Add double Mach reflection problem to compressible solver as ramp. (#31)

* Add double Mach reflection problem to compressible solver as ramp.

* flake8 check.

* flake8 check.

* flake8 check.

Author: mojiastonybrook

compressible/BC.py

@@ -13,6 +13,9 @@
 import compressible.eos as eos
 from util import msg
 
+import math
+import numpy as np
+
 
 def user(bc_name, bc_edge, variable, ccdata):
     """
@@ -132,5 +135,117 @@ def user(bc_name, bc_edge, variable, ccdata):
         else:
             msg.fail("error: hse BC not supported for xlb or xrb")
 
+    elif bc_name == "ramp":
+        # Boundary conditions for double Mach reflection problem
+
+        gamma = ccdata.get_aux("gamma")
+
+        if bc_edge == "xlb":
+            # lower x boundary
+            # inflow condition with post shock setup
+
+            v = ccdata.get_var(variable)
+            i = myg.ilo - 1
+            if variable in ["density", "x-momentum", "y-momentum", "energy"]:
+                val = inflow_post_bc(variable, gamma)
+                while i >= 0:
+                    v[i, :] = val
+                    i = i - 1
+            else:
+                v[:, :] = 0.0   # no source term
+
+        elif bc_edge == "ylb":
+            # lower y boundary
+            # for x > 1./6., reflective boundary
+            # for x < 1./6., inflow with post shock setup
+
+            if variable in ["density", "x-momentum", "y-momentum", "energy"]:
+                v = ccdata.get_var(variable)
+                j = myg.jlo - 1
+                jj = 0
+                while j >= 0:
+                    xcen_l = myg.x < 1.0/6.0
+                    xcen_r = myg.x >= 1.0/6.0
+                    v[xcen_l, j] = inflow_post_bc(variable, gamma)
+
+                    if variable == "y-momentum":
+                        v[xcen_r, j] = -1.0*v[xcen_r, myg.jlo+jj]
+                    else:
+                        v[xcen_r, j] = v[xcen_r, myg.jlo+jj]
+                    j = j - 1
+                    jj = jj + 1
+            else:
+                v = ccdata.get_var(variable)
+                v[:, :] = 0.0   # no source term
+
+        elif bc_edge == "yrb":
+            # upper y boundary
+            # time-dependent boundary, the shockfront moves with a 10 mach velocity forming an angle
+            # to the x-axis of 30 degrees clockwise.
+            # x coordinate of the grid is used to judge whether the cell belongs to pure post shock area,
+            # the pure pre shock area or the mixed area.
+
+            if variable in ["density", "x-momentum", "y-momentum", "energy"]:
+                v = ccdata.get_var(variable)
+                for j in range(myg.jhi+1, myg.jhi+myg.ng+1):
+                    shockfront_up = 1.0/6.0 + (myg.y[j] + 0.5*myg.dy*math.sqrt(3))/math.tan(math.pi/3.0) \
+                                    + (10.0/math.sin(math.pi/3.0))*ccdata.t
+                    shockfront_down = 1.0/6.0 + (myg.y[j] - 0.5*myg.dy*math.sqrt(3))/math.tan(math.pi/3.0) \
+                                    + (10.0/math.sin(math.pi/3.0))*ccdata.t
+                    shockfront = np.array([shockfront_down, shockfront_up])
+                    for i in range(myg.ihi+myg.ng+1):
+                        v[i, j] = 0.0
+                        cx_down = myg.x[i] - 0.5*myg.dx*math.sqrt(3)
+                        cx_up = myg.x[i] + 0.5*myg.dx*math.sqrt(3)
+                        cx = np.array([cx_down, cx_up])
+
+                        for sf in shockfront:
+                            for x in cx:
+                                if x < sf:
+                                    v[i, j] = v[i, j] + 0.25*inflow_post_bc(variable, gamma)
+                                else:
+                                    v[i, j] = v[i, j] + 0.25*inflow_pre_bc(variable, gamma)
+            else:
+                v = ccdata.get_var(variable)
+                v[:, :] = 0.0   # no source term
+
     else:
         msg.fail("error: bc type %s not supported" % (bc_name))
+
+
+def inflow_post_bc(var, g):
+    # inflow boundary condition with post shock setup
+    r_l = 8.0
+    u_l = 7.1447096
+    v_l = -4.125
+    p_l = 116.5
+    if var == "density":
+        vl = r_l
+    elif var == "x-momentum":
+        vl = r_l*u_l
+    elif var == "y-momentum":
+        vl = r_l*v_l
+    elif var == "energy":
+        vl = p_l/(g - 1.0) + 0.5*r_l*(u_l*u_l + v_l*v_l)
+    else:
+        vl = 0.0
+    return vl
+
+
+def inflow_pre_bc(var, g):
+    # pre shock setup
+    r_r = 1.4
+    u_r = 0.0
+    v_r = 0.0
+    p_r = 1.0
+    if var == "density":
+        vl = r_r
+    elif var == "x-momentum":
+        vl = r_r*u_r
+    elif var == "y-momentum":
+        vl = r_r*v_r
+    elif var == "energy":
+        vl = p_r/(g - 1.0) + 0.5*r_r*(u_r*u_r + v_r*v_r)
+    else:
+        vl = 0.0
+    return vl

compressible/problems/__init__.py

@@ -1 +1,2 @@
-__all__ = ['acoustic_pulse', 'advect', 'bubble', 'hse', 'kh', 'logo', 'quad', 'rt', 'rt2', 'sedov', 'sod', 'test']
+__all__ = ['acoustic_pulse', 'advect', 'bubble', 'hse', 'kh', 'logo', 'quad', 'rt', 'rt2', 'sedov', 'sod',
+           'test', 'ramp']

compressible/problems/_ramp.defaults

@@ -0,0 +1,19 @@
+# these defaults comes from the third test problem in 
+# Woodward and Colella. Journal of Computational Physics, 54, 115-174, 1984
+#
+# Also, the numbers is consitent with the ones in Castro
+#
+
+
+[ramp]
+rhol = 8.0                 ; post-shock initial density
+ul = 7.1447096             ; post-shock initial x-velocity
+vl = -4.125                ; post-shock initial y-velocity
+pl = 116.5                 ; post-shock initial pressure
+ 
+rhor = 1.4                 ; pre-shock initial density
+ur = 0.0                   ; pre-shock initial x-velocity
+vr = 0.0                   ; pre-shock initial y-velocity
+pr = 1.0                   ; pre-shock initial pressure
+ 
+ 

compressible/problems/inputs.ramp

@@ -0,0 +1,40 @@
+# simple inputs files for the double Mach reflection problem.
+
+[driver]
+max_steps = 12500
+tmax = 0.25
+
+
+[compressible]
+limiter = 2
+cvisc = 0.1
+
+[io]
+basename = double_mach_reflection_
+dt_out = 0.1
+
+
+[mesh]
+nx = 1024
+ny = 256
+xmax = 4.0
+ymax = 1.0
+
+xlboundary = ramp
+xrboundary = outflow
+
+ylboundary = ramp
+yrboundary = ramp
+
+
+[ramp]
+rhol = 8.0
+ul = 7.1447096
+vl = -4.125
+pl = 116.5
+
+rhor = 1.4
+ur = 0.0
+vr = 0.0
+pr = 1.0
+

compressible/problems/ramp.py

@@ -0,0 +1,88 @@
+from __future__ import print_function
+
+import sys
+
+import numpy as np
+import math
+
+import mesh.patch as patch
+from util import msg
+
+
+def init_data(my_data, rp):
+    """ initialize the double Mach reflection problem """
+
+    msg.bold("initializing the double Mach reflection problem...")
+
+    # make sure that we are passed a valid patch object
+    if not isinstance(my_data, patch.CellCenterData2d):
+        print("ERROR: patch invalid in ramp.py")
+        print(my_data.__class__)
+        sys.exit()
+
+    # get the density, momenta, and energy as separate variables
+    dens = my_data.get_var("density")
+    xmom = my_data.get_var("x-momentum")
+    ymom = my_data.get_var("y-momentum")
+    ener = my_data.get_var("energy")
+
+    # initialize the components, remember, that ener here is
+    # rho*eint + 0.5*rho*v**2, where eint is the specific
+    # internal energy (erg/g)
+
+    r_l = rp.get_param("ramp.rhol")
+    u_l = rp.get_param("ramp.ul")
+    v_l = rp.get_param("ramp.vl")
+    p_l = rp.get_param("ramp.pl")
+
+    r_r = rp.get_param("ramp.rhor")
+    u_r = rp.get_param("ramp.ur")
+    v_r = rp.get_param("ramp.vr")
+    p_r = rp.get_param("ramp.pr")
+
+    gamma = rp.get_param("eos.gamma")
+
+    energy_l = p_l/(gamma - 1.0) + 0.5*r_l*(u_l*u_l + v_l*v_l)
+    energy_r = p_r/(gamma - 1.0) + 0.5*r_r*(u_r*u_r + v_r*v_r)
+
+    # there is probably an easier way to do this, but for now, we
+    # will just do an explicit loop.  Also, we really want to set
+    # the pressue and get the internal energy from that, and then
+    # compute the total energy (which is what we store).  For now
+    # we will just fake this
+
+    myg = my_data.grid
+    dens[:, :] = 1.4
+
+    for j in range(myg.jlo, myg.jhi+1):
+        cy_up = myg.y[j] + 0.5*myg.dy*math.sqrt(3)
+        cy_down = myg.y[j] - 0.5*myg.dy*math.sqrt(3)
+        cy = np.array([cy_down, cy_up])
+
+        for i in range(myg.ilo, myg.ihi+1):
+            dens[i, j] = 0.0
+            xmom[i, j] = 0.0
+            ymom[i, j] = 0.0
+            ener[i, j] = 0.0
+
+            sf_up = math.tan(math.pi/3.0)*(myg.x[i] + 0.5*myg.dx*math.sqrt(3)-1.0/6.0)
+            sf_down = math.tan(math.pi/3.0)*(myg.x[i] - 0.5*myg.dx*math.sqrt(3)-1.0/6.0)
+            sf = np.array([sf_down, sf_up])   # initial shock front
+
+            for y in cy:
+                for shockfront in sf:
+                    if y >= shockfront:
+                        dens[i, j] = dens[i, j] + 0.25*r_l
+                        xmom[i, j] = xmom[i, j] + 0.25*r_l*u_l
+                        ymom[i, j] = ymom[i, j] + 0.25*r_l*v_l
+                        ener[i, j] = ener[i, j] + 0.25*energy_l
+                    else:
+                        dens[i, j] = dens[i, j] + 0.25*r_r
+                        xmom[i, j] = xmom[i, j] + 0.25*r_r*u_r
+                        ymom[i, j] = ymom[i, j] + 0.25*r_r*v_r
+                        ener[i, j] = ener[i, j] + 0.25*energy_r
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """
+    pass

compressible/simulation.py

@@ -112,6 +112,7 @@ 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
 
         bc, bc_xodd, bc_yodd = bc_setup(self.rp)