Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding Roe solver entropy fix and an exact Riemann solver #109

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Adding Roe solver entropy fix and an exact Riemann solver #109

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
11bb3d5
Adding Roe solver entropy fix and an exact Riemann solver
OwenEvans Dec 16, 2015
2f1e3e1
Adding line continuations
orevans Dec 18, 2015
11168f2
importing scipy newton solver
orevans Dec 18, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 131 additions & 13 deletions src/shallow_1D_py.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,16 +92,39 @@ def shallow_roe_1D(q_l,q_r,aux_l,aux_r,problem_data):
wave[0,1,:] = a2
wave[1,1,:] = a2 * (ubar + cbar)
s[1,:] = ubar + cbar


s_index = np.zeros((2,num_rp))
for m in xrange(num_eqn):
for mw in xrange(num_waves):
s_index[0,:] = s[mw,:]
amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]

if problem_data['efix']:
raise NotImplementedError("Entropy fix has not been implemented.")
else:
s_index = np.zeros((2,num_rp))
for m in xrange(num_eqn):
for mw in xrange(num_waves):
s_index[0,:] = s[mw,:]
amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:]
apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:]
# Compute eigenvalues of f'(q)
def lamb(i,q):
if (i == 1):
return q[1,:]/q[0,:] - np.sqrt(problem_data['grav']*q[0,:])
else:
return q[1,:]/q[0,:] + np.sqrt(problem_data['grav']*q[0,:])

# Compute intermediate state
q_m = q_l + wave[:,0,:]

# Check for transonic rarefactions
transonic_1 = (lamb(1,q_l) < 0.0) * (lamb(1,q_m) > 0.0)
transonic_2 = (lamb(2,q_m) < 0.0) * (lamb(2,q_r) > 0.0)

# Harten-Hyman entropy fix parameter
beta = np.zeros(num_rp)
beta[transonic_1] = (lamb(1,q_m[:,transonic_1]) - s[0,transonic_1]) / (lamb(1,q_m[:,transonic_1]) - lamb(1,q_l[:,transonic_1]))
beta[transonic_2] = (lamb(2,q_r[:,transonic_2]) - s[1,transonic_2]) / (lamb(2,q_r[:,transonic_2]) - lamb(2,q_m[:,transonic_2]))

# Update fluctuations
amdq[:,transonic_1] += beta[transonic_1] * lamb(1,q_l[:,transonic_1]) * wave[:,0,transonic_1] * (s[0,transonic_1] >= 0.0) + (beta[transonic_1] * lamb(1,q_l[:,transonic_1]) - s[0,transonic_1]) * wave[:,0,transonic_1] * (s[0,transonic_1] < 0.0)
amdq[:,transonic_2] += beta[transonic_2] * lamb(2,q_m[:,transonic_2]) * wave[:,1,transonic_2] * (s[1,transonic_2] >= 0.0) + (beta[transonic_2] * lamb(2,q_m[:,transonic_2]) - s[1,transonic_2]) * wave[:,1,transonic_2] * (s[1,transonic_2] < 0.0)
apdq[:,transonic_1] += (1. - beta[transonic_1]) * lamb(1,q_m[:,transonic_1]) * wave[:,0,transonic_1] * (s[0,transonic_1] < 0.0) + ((1. - beta[transonic_1]) * lamb(1,q_m[:,transonic_1]) - s[0,transonic_1]) * wave[:,0,transonic_1] * (s[0,transonic_1] >= 0.0)
apdq[:,transonic_2] += (1. - beta[transonic_2]) * lamb(2,q_r[:,transonic_2]) * wave[:,1,transonic_2] * (s[1,transonic_2] < 0.0) + ((1. - beta[transonic_2]) * lamb(2,q_r[:,transonic_2]) - s[1,transonic_2]) * wave[:,1,transonic_2] * (s[1,transonic_2] >= 0.0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add some continuations to these lines so they are not quite as long?


return wave, s, amdq, apdq

Expand Down Expand Up @@ -233,9 +256,104 @@ def shallow_fwave_1d(q_l, q_r, aux_l, aux_r, problem_data):
def shallow_exact_1D(q_l,q_r,aux_l,aux_r,problem_data):
r"""
Exact shallow water Riemann solver

"""
num_eq = 2
num_waves = 2

.. warning::
This solver has not been implemented.
# Parameters
g = problem_data['grav']

# Array shapes
num_rp = q_l.shape[1]

# Output arrays
wave = np.zeros( (num_eqn, num_waves, num_rp) )
s = np.zeros( (num_waves, num_rp) )
sm = np.zeros( (num_waves, num_rp) )
amdq = np.zeros( (num_eqn, num_rp) )
apdq = np.zeros( (num_eqn, num_rp) )
h_m, u_m = np.zeros(num_rp), np.zeros(num_rp)

# Set heights and velocities
h_l, h_r = q_l[0,:], q_r[0,:]
u_l, u_r = q_l[1,:] / q_l[0,:], q_r[1,:] / q_r[0,:]

# Set intermediate states
h_m, u_m = np.zeros(num_rp), np.zeros(num_rp)

# Functions defined in George 2008 (Appendix B)
def phi(x, h_p):
if (x <= h_p):
return 2.*(np.sqrt(g*x) - np.sqrt(g*h_p))
else:
return (x - h_p)*np.sqrt(0.5*g*(1./x + 1./h_p))

def psi(x, h_l, h_r, u_l, u_r):
return phi(x, h_r) + phi(x, h_l) + u_r - u_l

psi_min, psi_max = np.zeros(num_rp), np.zeros(num_rp)

# Newton solve to find intermediate state q_m
for i in xrange(num_rp):
h_m[i] = newton(psi, 1.e-3, args=(h_l[i],h_r[i],u_l[i],u_r[i]))
u_m[i] = (u_l[i] - phi(h_m[i], h_l[i]))
h_min, h_max = min(h_l[i], h_r[i]), max(h_l[i], h_r[i])
psi_min[i], psi_max[i] = psi(h_min, h_l[i], h_r[i], u_l[i], u_r[i]), psi(h_max, h_l[i], h_r[i], u_l[i], u_r[i])

# Compute Roe and right and left speeds
ubar = ( (q_l[1,:]/np.sqrt(q_l[0,:]) + q_r[1,:]/np.sqrt(q_r[0,:])) / (np.sqrt(q_l[0,:]) + np.sqrt(q_r[0,:])) )
cbar = np.sqrt(0.5*g*(q_l[0,:] + q_r[0,:]))
u_r = q_r[1,:]/q_r[0,:]
c_r = np.sqrt(g*q_r[0,:])
u_l = q_l[1,:]/q_l[0,:]
c_l = np.sqrt(g*q_l[0,:])

# Compute Einfeldt speeds
s_index = np.empty((4,num_rp))
s_index[0,:] = ubar+cbar
s_index[1,:] = ubar-cbar
s_index[2,:] = u_l + c_l
s_index[3,:] = u_l - c_l
s[0,:] = np.min(s_index,axis=0)
s_index[2,:] = u_r + c_r
s_index[3,:] = u_r - c_r
s[1,:] = np.max(s_index,axis=0)

"""
raise NotImplementedError("The exact swe solver has not been implemented.")
# Determine characteristic structure for each Riemann problem
all_shock = (psi_min <= psi_max)*(psi_max <= 0.0)
one_rar = (psi_min < 0.0)*(psi_max >= 0.0)*(h_l > h_r)
two_rar = (psi_min < 0.0)*(psi_max > 0.0)*(h_l < h_r)
all_rar = (0.0 <= psi_min)*(psi_min < psi_max)

# qt1 and qt2 are solutions that correspond to transonic rarefactions in the 1- and 2-wave, respectively.
qt1, qt2 = np.zeros( (num_eqn, num_rp) ), np.zeros( (num_eqn, num_rp) )
qt1[0,:] =(1./(9.*g))*(u_l + 2.*np.sqrt(g*h_l))**2
qt1[1,:] = qt1[0,:]*(u_l + 2.*(np.sqrt(g*h_l) - np.sqrt(g*qt1[0,:])))
qt2[0,:] =(1./(9.*g))*(u_r - 2.*np.sqrt(g*h_r))**2
qt2[1,:] = qt2[0,:]*(u_r + 2.*(np.sqrt(g*qt2[0,:]) - np.sqrt(g*h_r)))

# Compute q_m and associated eigenvalues
q_m = np.zeros( (num_eqn, num_rp ) )
q_m[0,:], q_m[1,:] = h_m, h_m*u_m
sm[0,:] = q_m[1,:]/q_m[0,:] - np.sqrt(g*q_m[0,:])
sm[1,:] = q_m[1,:]/q_m[0,:] + np.sqrt(g*q_m[0,:])

# Compute waves
wave[:,0,:] = q_m - q_l
wave[:,1,:] = q_r - q_m

# Evaluate q at the interface
q = 0.5*(q_l + q_r)
q[:,all_shock] = q_r[:, all_shock]*(s[1,all_shock] <= 0.0) + q_l[:,all_shock]*(s[0,all_shock] >= 0.0) + q_m[:,all_shock]*(s[0,all_shock] < 0.0)*(0.0 < s[1,all_shock])
q[:,one_rar] = (q_m[:,one_rar]*(sm[0,one_rar] <= 0.0) + qt1[:,one_rar]*(sm[0,one_rar] >= 0.0))*(s[0,one_rar] <= 0.0)*(0.0 <= s[1,one_rar]) + q_r[:,one_rar]*(s[1,one_rar] < 0.0) + q_l[:,one_rar]*(s[0,one_rar] > 0.0)
q[:,two_rar] = (q_m[:,two_rar]*(sm[1,two_rar] >= 0.0) + qt2[:,two_rar]*(sm[1,two_rar] < 0.0))*(s[0,two_rar] <= 0.0)*(0.0 <= s[1,two_rar]) + q_r[:,two_rar]*(s[1,two_rar] < 0.0) + q_l[:,two_rar]*(s[0,two_rar] > 0.0)
q[:,all_rar] = q_m[:,all_rar]*(sm[0,all_rar] <= 0.0)*(0.0 <= sm[1,all_rar]) + qt1[:,all_rar]*(sm[0,all_rar] > 0.0)*(s[0,all_rar] <= 0.0) + qt2[:,all_rar]*(sm[1,all_rar] < 0.0)*(s[1,all_rar] >= 0.0) + q_r[:,all_rar]*(s[1,all_rar] < 0.0) + q_l[:,all_rar]*(s[0,all_rar] > 0.0)

# Compute fluctuations amdq = f(q) and apdq = -f(q)
f = np.zeros( (num_eqn, num_rp) )
f[0,:] = q[1,:]
f[1,:] = ((q[1,:])**2)/q[0,:] + 0.5*g*(q[0,:])**2
amdq, apdq = f, -f

return wave, s, amdq, apdq