the exact Riemann solver

Author: zingale

riemann_exact.py

@@ -0,0 +1,162 @@
+# solve the Riemann problem for a gamma-law gas
+
+import numpy as np
+import scipy.optimize as optimize
+
+class State(object):
+    """ a simple object to hold a primitive variable state """
+
+    def __init__(self, p=1.0, u=0.0, rho=1.0):
+        self.p = p
+        self.u = u
+        self.rho = rho
+
+    def __str__(self):
+        return "rho: {}; u: {}; p: {}".format(self.rho, self.u, self.p)
+
+class RiemannProblem(object):
+    """ a class to define a Riemann problem.  It takes a left
+        and right state.  Note: we assume a constant gamma """
+
+    def __init__(self, left_state, right_state, gamma=1.4):
+        self.left = left_state
+        self.right = right_state
+        self.gamma = gamma
+
+        self.ustar = None
+        self.pstar = None
+
+    def u_hugoniot(self, p, side):
+        """define the Hugoniot curve, u(p)."""
+
+        if side == "left":
+            state = self.left
+            s = 1.0
+        elif side == "right":
+            state = self.right
+            s = -1.0
+
+        c = np.sqrt(self.gamma*state.p/state.rho)
+
+        if p < state.p:
+            # rarefaction
+            u = state.u + s*(2.0*c/(self.gamma-1.0))* \
+                (1.0 - (p/state.p)**((self.gamma-1.0)/(2.0*self.gamma)))
+        else:
+            # shock
+            beta = (self.gamma+1.0)/(self.gamma-1.0)
+            u = state.u + s*(2.0*c/np.sqrt(2.0*self.gamma*(self.gamma-1.0)))* \
+                (1.0 - p/state.p)/np.sqrt(1.0 + beta*p/state.p)
+
+        return u
+
+    def find_star_state(self, p_min=0.001, p_max=1000.0):
+        """ root find the Hugoniot curve to find ustar, pstar """
+
+        # we need to root-find on
+        self.pstar = optimize.brentq(
+            lambda p: self.u_hugoniot(p, "left") - self.u_hugoniot(p, "right"),
+            p_min, p_max)
+        self.ustar = self.u_hugoniot(self.pstar, "left")
+
+
+    def shock_solution(self, sgn, state):
+        """return the interface solution considering a shock"""
+
+        p_ratio = self.pstar/state.p
+        c = np.sqrt(self.gamma*state.p/state.rho)
+
+        # Toro, eq. 4.52 / 4.59
+        S = state.u + sgn*c*np.sqrt(0.5*(self.gamma + 1.0)/self.gamma*p_ratio +
+                                    0.5*(self.gamma - 1.0)/self.gamma)
+
+        # are we to the left or right of the shock?
+        if (self.ustar < 0 and S < 0) or (self.ustar > 0 and S > 0):
+            # R/L region
+            solution = state
+        else:
+            # * region -- get rhostar from Toro, eq. 4.50 / 4.57
+            gam_fac = (self.gamma - 1.0)/(self.gamma + 1.0)
+            rhostar = state.rho * (p_ratio + gam_fac)/(gam_fac * p_ratio + 1.0)
+            solution = State(rho=rhostar, u=self.ustar, p=self.pstar)
+
+        return solution
+
+    def rarefaction_solution(self, sgn, state):
+        """return the interface solution considering a rarefaction wave"""
+
+        # find the speed of the head and tail of the rarefaction fan
+
+        # isentropic (Toro eq. 4.54 / 4.61)
+        p_ratio = self.pstar/state.p
+        c = np.sqrt(self.gamma*state.p/state.rho)
+        cstar = c*p_ratio**((self.gamma-1.0)/(2*self.gamma))
+
+        lambda_head = state.u + sgn*c
+        lambda_tail = self.ustar + sgn*cstar
+
+        if (sgn > 0 and lambda_head < 0) or (sgn < 0 and lambda_head > 0):
+            # R/L region
+            solution = state
+
+        elif (sgn > 0 and lambda_tail > 0) or (sgn < 0 and lambda_tail < 0):
+            # * region, we use the isentropic density (Toro 4.53 / 4.60)
+            solution = State(rho = state.rho*p_ratio**(1.0/self.gamma),
+                             u = self.ustar, p = self.pstar)
+
+        else:
+            # we are in the fan -- Toro 4.56 / 4.63
+            rho = state.rho * (2/(self.gamma + 1.0) -
+                               sgn*gam_fac*state.u/c)**(2.0/(self.gamma-1.0))
+            u = 2.0/(self.gamma + 1.0) * ( -sgn*c + 0.5*(self.gamma - 1.0)*state.u)
+            p = state.p * (2/(self.gamma + 1.0) -
+                           sgn*gam_fac*state.u/c)**(2.0*self.gamma/(self.gamma-1.0))
+            solution = State(rho=rho, u=u, p=p)
+
+        return solution
+
+    def sample_solution(self):
+        """given the star state (ustar, pstar), find the state on the interface"""
+
+        gam_fac = (self.gamma - 1.0)/(self.gamma + 1.0)
+
+        if self.ustar < 0:
+            # we are in the R* or R region
+            state = self.right
+            sgn = 1.0
+        else:
+            # we are in the L* or L region
+            state = self.left
+            sgn = -1.0
+
+        # is the non-contact wave a shock or rarefaction?
+        if self.pstar > state.p:
+            # compression! we are a shock
+            solution = self.shock_solution(sgn, state)
+
+        else:
+            # rarefaction
+            solution = self.rarefaction_solution(sgn, state)
+
+        return solution
+
+def cons_flux(state, v):
+    """ given an interface state, return the conservative flux"""
+    flux = np.zeros((v.nvar), dtype=np.float64)
+
+    flux[v.urho] = state.rho * state.u
+    flux[v.umx] = flux[v.urho] * state.u + state.p
+    flux[v.uener] = (0.5 * state.rho * state.u**2 +
+                     state.p/(v.gamma - 1.0) + state.p) * state.u
+    return flux
+
+if __name__ == "__main__":
+
+    q_l = State(rho=1.0, u=0.0, p=1.0)
+    q_r = State(rho=0.125, u=0.0, p=0.1)
+
+    rp = RiemannProblem(q_l, q_r, gamma=1.4)
+
+    rp.find_star_state()
+    q_int = rp.sample_solution()
+    print(q_int)