From 11bb3d51f36f6f7726e08c7d1725df23e22855a6 Mon Sep 17 00:00:00 2001 From: Owen Evans Date: Wed, 16 Dec 2015 11:23:20 -0500 Subject: [PATCH 1/3] Adding Roe solver entropy fix and an exact Riemann solver --- src/shallow_1D_py.py | 144 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 131 insertions(+), 13 deletions(-) diff --git a/src/shallow_1D_py.py b/src/shallow_1D_py.py index 6a6b9fd9..a426393c 100644 --- a/src/shallow_1D_py.py +++ b/src/shallow_1D_py.py @@ -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) return wave, s, amdq, apdq @@ -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 From 2f1e3e1a8cd105721593588b0fe9beae02551ee8 Mon Sep 17 00:00:00 2001 From: Owen Evans Date: Fri, 18 Dec 2015 11:20:59 -0500 Subject: [PATCH 2/3] Adding line continuations --- src/shallow_1D_py.py | 72 ++++++++++++++++++++++++++++++++------------ 1 file changed, 52 insertions(+), 20 deletions(-) diff --git a/src/shallow_1D_py.py b/src/shallow_1D_py.py index a426393c..f42efede 100644 --- a/src/shallow_1D_py.py +++ b/src/shallow_1D_py.py @@ -117,14 +117,32 @@ def lamb(i,q): # 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])) + 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) + 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) return wave, s, amdq, apdq @@ -176,11 +194,11 @@ def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data): # Compute middle state q_hat = np.empty((2,num_rp)) - q_hat[0,:] = ((q_r[1,:] - q_l[1,:] - s[1,:] * q_r[0,:] - + s[0,:] * q_l[0,:]) / (s[0,:] - s[1,:])) - q_hat[1,:] = ((q_r[1,:]**2/q_r[0,:] + 0.5 * problem_data['grav'] * q_r[0,:]**2 - - (q_l[1,:]**2/q_l[0,:] + 0.5 * problem_data['grav'] * q_l[0,:]**2) - - s[1,:] * q_r[1,:] + s[0,:] * q_l[1,:]) / (s[0,:] - s[1,:])) + q_hat[0,:] = ((q_r[1,:] - q_l[1,:] - s[1,:] * q_r[0,:] \ + + s[0,:] * q_l[0,:]) / (s[0,:] - s[1,:])) + q_hat[1,:] = ((q_r[1,:]**2/q_r[0,:] + 0.5 * problem_data['grav'] * q_r[0,:]**2 \ + - (q_l[1,:]**2/q_l[0,:] + 0.5 * problem_data['grav'] * q_l[0,:]**2) \ + - s[1,:] * q_r[1,:] + s[0,:] * q_l[1,:]) / (s[0,:] - s[1,:])) # Compute each family of waves wave[:,0,:] = q_hat - q_l @@ -273,7 +291,6 @@ def shallow_exact_1D(q_l,q_r,aux_l,aux_r,problem_data): 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,:] @@ -296,13 +313,16 @@ def psi(x, h_l, h_r, u_l, u_r): # 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])) + h_m[i] = scipy.optimize.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]) + psi_min[i] = psi(h_min, h_l[i], h_r[i], u_l[i], u_r[i]) + psi_max[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,:])) ) + 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,:]) @@ -326,7 +346,7 @@ def psi(x, h_l, h_r, u_l, u_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 and qt2 are 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,:]))) @@ -345,10 +365,22 @@ def psi(x, h_l, h_r, u_l, u_r): # 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) + 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) ) From 11168f29b6174103358bc4a11c5a2c127a41371d Mon Sep 17 00:00:00 2001 From: Owen Evans Date: Fri, 18 Dec 2015 15:54:51 -0500 Subject: [PATCH 3/3] importing scipy newton solver --- src/shallow_1D_py.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/shallow_1D_py.py b/src/shallow_1D_py.py index f42efede..742c6a41 100644 --- a/src/shallow_1D_py.py +++ b/src/shallow_1D_py.py @@ -37,6 +37,7 @@ # ============================================================================ import numpy as np +from scipy.optimize import newton num_eqn = 2 num_waves = 2 @@ -313,7 +314,7 @@ def psi(x, h_l, h_r, u_l, u_r): # Newton solve to find intermediate state q_m for i in xrange(num_rp): - h_m[i] = scipy.optimize.newton(psi, 1.e-3, \ + 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])