In [1]:
from sympy import init_session
init_session()
In [2]:
from sympy.abc import rho, tau, alpha
rho, tau, c, h, p = symbols("rho tau c h p", real=True, positive=True)
re = symbols(r"(\rho{}e)", real=True, positive=True)
ge = symbols(r"\gamma_e", real=True, positive=True)
alpha, u = symbols("alpha u", real=True)
The two routines below simplify the analysis of the eigensystem by simultaneously finding the set of orthonormal left and right eigenvectors of a matrix from the primitive form of the hydro equations
In [3]:
plus = u + c
zero = u
minus = u - c
class Eigenvector(object):
def __init__(self, name, ev, r, l=None):
self.name = name
if name == "minus":
self.d = 0
elif name == "zero":
self.d = 1
elif name == "plus":
self.d = 2
else:
self.d = None
self.ev = ev
self.l = l
self.r = r
def __lt__(self, other):
return self.d < other.d
def __str__(self):
return "{} wave, r = {}, l = {}".format(self.eigenvalue, self.r, self.l)
def eigensystem(A, suba=None, subb=None):
# get the left and right eigenvectors that diagonalize the system.
# it is best to use sympy diagonalize() for this purpose than getting
# the left and right eigenvectors independently.
e = []
R, D = A.diagonalize()
# the columns of R are the right eigenvectors and the diagonal
# element of D is the corresponding eigenvalues
for n in range(A.shape[0]):
r = R.col(n)
ev = D[n,n]
#print("here", r, ev)
if suba is not None and subb is not None:
ev = ev.subs(suba, subb)
# which eigenvalue are we?
if simplify(ev - minus) == 0:
name = "minus"
elif simplify(ev - plus) == 0:
name = "plus"
elif simplify(ev - zero) == 0:
name = "zero"
else:
return None
# normalize the right eigenvector
v = r[0]
if v != 0:
r = r/v
if suba is not None and subb is not None:
r = simplify(r.subs(suba, subb))
e.append(Eigenvector(name=name, ev=ev, r=r))
# now sort the system from smallest (u-c) to largest (u+c)
e.sort()
# now let's construct the R with this sorting
for n in range(A.shape[0]):
R[:,n] = e[n].r
# the left eigenvector matrix, L, is just the inverse
L = R**-1
for n in range(A.shape[0]):
e[n].l = L.row(n)
return e
The Euler equations in primitive variable form, $q = (\rho, u, p, (\rho e))^\intercal$ appear as: \begin{align} \frac{\partial \rho}{\partial t} &= -u \frac{\partial \rho}{\partial x} - \rho \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \\ \frac{\partial (\rho e)}{\partial t} &= -u \frac{\partial (\rho e)}{\partial x} - \rho h \frac{\partial u}{\partial x} \end{align}
In vector form, we have: $$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{cccc} u & \rho & 0 & 0\\ 0 & u & 1/\rho & 0\\ 0 & \rho c^2 & u & 0\\ 0 & \rho h & 0 & u\end{array} \right ) $$The sound speed is related to the adiabatic index, $\Gamma_1$, as $c^2 = \Gamma_1 p /\rho$.
We can represent this matrix symbolically in SymPy and explore its eigensystem.
In [4]:
q = Matrix([rho, u, p, re]).transpose()
A = Matrix([[u, rho, 0, 0], [0, u, rho**-1, 0], [0, c**2 * rho, u, 0], [0, rho*h, 0, u]])
A
Out[4]:
The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.
In [5]:
A.eigenvals()
Out[5]:
We see that there are 2 eigenvalues $u$ -- the addition of $(\rho e)$ to the system adds this degeneracy.
The right eigenvectors are defined for a given eigenvalue, $\lambda$, as: $$A r = \lambda r$$ and the left eigenvectors satisfy: $$l A = \lambda l$$
Note that the left and right eigenvectors are orthogonal to those corresponding to a different eigenvalue, and usually normalized so: $$l^i \cdot r^j = \delta_{ij}$$
In [6]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)
In [7]:
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
print(e.name)
display(e.r, e.l)
The final interface state is writen by projecting the jump in primitive variables, $\Delta q$, into characteristic variables (as $l \cdot \Delta q$), and then adding up all the jumps that reach the interface.
The convention is to write $\beta^\nu = l^\nu \cdot \Delta q$, where the superscript identifies which eigenvalue (and corresponding eigenvectors) we are considering. Note, that often a reference state is used, and the jump, $\Delta q$, will be the difference with respect to this reference state. For PPM, the $\Delta q$ will take the form of the integral under the parabola over the range that each wave can reach.
The final interface state is then: $$ q_\mathrm{int} = q - \sum_\nu \beta^\nu r^\nu $$
The tracing projects the primitive variables into characteristic variables by defining $$ \Delta q = \left ( \begin{array}{c} \Delta \rho \\ \Delta u \\ \Delta p \\ \Delta (\rho e) \end{array} \right ) $$ and then $\beta^\nu = l^\nu \cdot \Delta q$
In [8]:
from sympy.abc import delta
dr = symbols(r"\Delta\rho")
du = symbols(r"\Delta{}u")
dp = symbols(r"\Delta{}p")
dre = symbols(r"\Delta(\rho{}e)")
rhoi = symbols(r"\rho_\mathrm{int}")
ui = symbols(r"u_\mathrm{int}")
pri = symbols(r"p_\mathrm{int}")
rei = symbols(r"(\rho{}e)_\mathrm{int}")
# this is the jump
dq = Matrix([[dr, du, dp, dre]]).transpose()
# this is the interface state vector
qint = Matrix([[rhoi, ui, pri, rei]]).transpose()
Now compute the $\beta$s
In [10]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
print(e.name)
beta = e.l.dot(dq)
display(Eq(betas[n], simplify(beta)))
and now the final interface states
In [11]:
for n in range(len(eigen)):
rhs = q[n]
for m in range(len(eigen)):
rhs -= betas[m]*eigen[m].r[n]
display(Eq(qint[n],rhs))
We can define $\gamma_e = p/(\rho e) + 1$ and differentiate it to get: $$ \frac{\partial \gamma_e}{\partial t} = -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x} $$ The original CG paper used $\tau = 1/\rho$ in place of density. With this, the continuity equation becomes: $$ \frac{\partial \tau}{\partial t} = -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} $$
The Euler equations with this set of primitive variables, $q = (\tau, u, p, \gamma_e)^\intercal$ appear as: \begin{align} \frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c^2 \frac{\partial u}{\partial x} \\ \frac{\partial \gamma_e}{\partial t} &= -u \frac{\partial \gamma_e}{\partial x} + (\gamma_e - 1)(\gamma_e - \Gamma_1) \frac{\partial u}{\partial x} \end{align}
For convenience, we define $$ \alpha = (\gamma_e - 1)(\gamma_e - \Gamma_1) $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{cccc} u & -\tau & 0 & 0\\ 0 & u & \tau & 0\\ 0 & c^2/\tau & u & 0\\ 0 & -\alpha & 0 & u\end{array} \right ) $$We can represent this matrix symbolically in SymPy and explore its eigensystem.
In [12]:
q = Matrix([[tau, u, p, ge]]).transpose()
A = Matrix([[u, -tau, 0, 0], [0, u, tau, 0], [0, c**2/tau, u, 0], [0, -alpha, 0, u]])
A
Out[12]:
The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.
In [13]:
A.eigenvals()
Out[13]:
We see that there are 2 eigenvalues $u$ -- the addition of $\gamma_e$ to the system adds this degeneracy.
In [14]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A)
In [15]:
# printing them out for inspection
from sympy.printing.mathml import mathml
for e in eigen:
print(e.name)
display(e.r, e.l)
The tracing projects the primitive variables into characteristic variables by defining $$ \Delta q = \left ( \begin{array}{c} \Delta \tau \\ \Delta u \\ \Delta p \\ \Delta \gamma_e \end{array} \right ) $$ and then $\beta^\nu = l^\nu \cdot \Delta q$
In [16]:
dtau = symbols(r"\Delta\tau")
dge = symbols(r"\Delta\gamma_e")
dq = Matrix([[dtau, du, dp, dge]]).transpose()
taui = symbols(r"\tau_\mathrm{int}")
gei = symbols(r"(\gamma_e)_\mathrm{int}")
qint = Matrix([[taui, ui, pri, gei]]).transpose()
In [18]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
print(e.name)
beta = e.l.dot(dq)
display(Eq(betas[n], simplify(beta)))
and now the final interface states
In [19]:
for n in range(len(eigen)):
rhs = q[n]
for m in range(len(eigen)):
rhs -= betas[m]*eigen[m].r[n]
display(Eq(qint[n],rhs))
Following Zhang et al. (2011), the equations of gray radiation hydrodynamics with primitive variables $q = (\rho, u, p, (\rho e)_g, E_r)^\intercal$ are: \begin{align} \frac{\partial \rho}{\partial t} &= -u\frac{\partial \rho}{\partial x} - \rho \frac{\partial u}{\partial x} \ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \frac{1}{\rho} \frac{\partial p}{\partial x}
- \frac{\lambda_f}{\rho} \frac{\partial E_r}{\partial x}\\
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \rho c_g^2 \frac{\partial u}{\partial x} \ \frac{\partial (\rho e)_g}{\partial t} &= -u \frac{\partial (\rho e)_g}{\partial x} - \rho h \frac{\partial u}{\partial x} \ \frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x} \end{align} where $(\rho e)_g$ is the gas internal energy density, $h_g$ is the gas specific enthalpy, $c_g$ is the gas sound speed (obeying $c_g^2 = \Gamma_1 p /\rho$), $E_r$ is the radiation energy density, $f$ is the Eddington factor, and $\lambda_f$ is the flux limiter.
Following Zhang et al., we make the approximation that $$ \frac{3-f}{2} = \lambda_f + 1 $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{ccccc} u & \rho & 0 & 0 & 0 \\ 0 & u & 1/\rho & 0 & \lambda_f/\rho\\ 0 & \rho c_g^2 & u & 0 & 0 \\ 0 & \rho h_g & 0 & u & 0 \\ 0 & (3-f)E_r/2 & 0 & 0 & u\end{array} \right ) $$We can represent this matrix symbolically in SymPy and explore its eigensystem.
In [20]:
cg = symbols(r"c_g", real=True, positive=True)
hg = symbols(r"h_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True)
f = symbols("f", real=True, positive=True)
q = Matrix([[rho, u, p, re, Er]]).transpose()
A = Matrix([[u, rho, 0, 0, 0],
[0, u, rho**-1, 0, lf/rho],
[0, rho*cg**2, u, 0, 0],
[0, rho*hg, 0, u, 0],
[0, (lf+1)*Er, 0, 0, u]])
A
Out[20]:
The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.
In [21]:
A.eigenvals()
Out[21]:
We see that there are 3 eigenvalues $u$. We identify the total sound speed (radiation + gas) as: $$c^2 = c_g^2 + (\lambda_f + 1)\frac{\lambda_f E_r}{\rho}$$ We can simplify these by substituting in that relationship
In [22]:
cc = c**2 - (lf +1)*lf*Er/rho
evs = A.eigenvals()
for e in evs.keys():
display(powsimp(simplify(e.subs(cg**2, cc))))
We see that these are the same form of the eigenvalues for the pure hydrodynamics system
In [23]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))
In [24]:
# printing them out for inspection
for e in eigen:
print(e.name)
display(e.r)
display(simplify(e.l))
The tracing projects the primitive variables into characteristic variables by defining $\beta^\nu = l^\nu \cdot \Delta q$
In [25]:
dEr = symbols(r"\Delta{}E_r")
Eri = symbols(r"{E_r}_\mathrm{int}")
dq = Matrix([[dr, du, dp, dre, dEr]]).transpose()
qint = Matrix([[rhoi, ui, pri, rei, Eri]])
In [27]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
print(e.name)
beta = e.l.dot(dq)
display(Eq(betas[n], simplify(beta)))
and now the final interface states
In [28]:
for n in range(len(eigen)):
rhs = q[n]
for m in range(len(eigen)):
rhs -= betas[m]*eigen[m].r[n]
display(Eq(qint[n],rhs))
In [ ]:
We now look at the same system with a different auxillary thermodynamic variable (as we did with pure hydro), using $q = (\tau, u, p, {\gamma_e}_g, E_r)^\intercal$: \begin{align} \frac{\partial \tau}{\partial t} &= -u\frac{\partial \tau}{\partial x} + \tau \frac{\partial u}{\partial x} \ \frac{\partial u}{\partial t} &= -u \frac{\partial u}{\partial x} - \tau \frac{\partial p}{\partial x}
- \tau \lambda_f \frac{\partial E_r}{\partial x}\\
\frac{\partial p}{\partial t} &= -u \frac{\partial p}{\partial x} - \frac{c_g^2}{\tau} \frac{\partial u}{\partial x} \ \frac{\partial {\gamma_e}_g}{\partial t} &= -u \frac{\partial {\gamma_e}_g}{\partial x} + \alpha \frac{\partial u}{\partial x} \ \frac{\partial E_r}{\partial t} &= -\frac{3-f}{2} E_r \frac{\partial u}{\partial x} - \left ( \frac{3-f}{2} - \lambda_f\right ) u \frac{\partial E_r}{\partial x} \end{align} here, ${\gamma_e}_g$ is defined solely in terms of the gas pressure and energy, the remaining variables have the same meaning as above. We again make the approximation that $$ \frac{3-f}{2} = \lambda_f + 1 $$ and then in vector form, we have: $$q_t + A(q) q_x = 0$$
with the matrix $A(q)$:
$$A(q) = \left ( \begin{array}{ccccc} u & -\tau & 0 & 0 & 0 \\ 0 & u & \tau & 0 & \tau \lambda_f\\ 0 & c_g^2/\tau & u & 0 & 0 \\ 0 & -\alpha & 0 & u & 0 \\ 0 & (3-f)E_r/2 & 0 & 0 & u\end{array} \right ) $$We can represent this matrix symbolically in SymPy and explore its eigensystem.
In [29]:
geg = symbols(r"{\gamma_e}_g", real=True, positive=True)
Er = symbols(r"E_r", real=True, positive=True)
lf = symbols(r"\lambda_f", real=True, positive=True)
cg = symbols(r"c_g", real=True, positive=True)
q = Matrix([[tau, u, p, geg, Er]]).transpose()
A = Matrix([[u, -tau, 0, 0, 0],
[0, u, tau, 0, tau*lf],
[0, cg**2/tau, u, 0, 0],
[0, -alpha, 0, u, 0],
[0, (lf+1)*Er, 0, 0, u]])
A
Out[29]:
The eigenvalues are the speeds at which information propagates with. SymPy returns them as a dictionary, giving the multiplicity for each eigenvalue.
In [30]:
A.eigenvals()
Out[30]:
We see that there are 3 eigenvalues $u$. We identify the total sound speed (radiation + gas) as: $$c^2 = c_g^2 + (\lambda_f + 1)\lambda_f E_r\tau$$ We can simplify these by substituting in that relationship
In [31]:
cc = c**2 - (lf +1)*lf*Er*tau
evs = A.eigenvals()
for e in evs.keys():
display(powsimp(simplify(e.subs(cg**2, cc))))
We see that these are the same form of the eigenvalues for the pure hydrodynamics system
In [32]:
# we use the helper rountines above to find the orthogonal left and right eigenvectors
eigen = eigensystem(A, suba=cg, subb=sqrt(cc))
In [33]:
# printing them out for inspection
for e in eigen:
print(e.name)
display(e.r)
display(simplify(e.l))
The tracing projects the primitive variables into characteristic variables by defining $\beta^\nu = l^\nu \cdot \Delta q$
In [34]:
dgeg = symbols(r"\Delta{\gamma_e}_g")
dEr = symbols(r"\Delta{}E_r")
dq = Matrix([[dtau, du, dp, dgeg, dEr]]).transpose()
Eri = symbols(r"{E_r}_\mathrm{int}")
qint = Matrix([[taui, ui, pri, gei, Eri]])
In [36]:
betas = [symbols(r"\beta^-")]
for n in range(len([e for e in eigen if e.d == 1])):
betas += [symbols(r"\beta^0_{}".format(n))]
betas += [symbols(r"\beta^+")]
for n, e in enumerate(eigen):
print(e.name)
beta = e.l.dot(dq)
display(Eq(betas[n], simplify(beta)))
and now the final interface states
In [37]:
for n in range(len(eigen)):
rhs = q[n]
for m in range(len(eigen)):
rhs -= betas[m]*eigen[m].r[n]
display(Eq(qint[n],rhs))
In [ ]: