diff options
| -rwxr-xr-x | makeBoundary.py | 175 |
1 files changed, 89 insertions, 86 deletions
diff --git a/makeBoundary.py b/makeBoundary.py index 8dc7ad0..639160f 100755 --- a/makeBoundary.py +++ b/makeBoundary.py @@ -16,6 +16,95 @@ def read_definition(filename): ddict[items[0]] = ast.literal_eval(items[1]) return ddict +def conveyance(numH, n_co, xregion, zregion, zmin, zmax): + p_i = [] # wetted perimeter + A_i = [] # area + r_h = [] # hydraulic radius + h_i = [] # list of heights + K_i = [] # conveyance + Q_i = [] # discharge + x_sub = [[] for i in range(numH)] # list of x values in subregion + z_sub = [[] for i in range(numH)] # list of z values in subregion + for i in range(numH): + h_i.append(zmin + (i+1)*(zmax-zmin)/numH) + #print(zregion[zregion < h_i[i]]) + booleanArray = zregion < h_i[i] + #print(booleanArray[i]) + x_sub[i] += list(xregion[booleanArray]) + z_sub[i] += list(zregion[booleanArray]) + for interval in range(len(xregion)-1): + if booleanArray[interval+1] != booleanArray[interval]: + x_extra = xregion[interval] + (h_i[i]-zregion[interval])*(xregion[interval+1]-xregion[interval])/(zregion[interval+1]-zregion[interval]) + bisect.insort(x_sub[i], x_extra) # add intercept value + ind_x = x_sub[i].index(x_extra) + z_sub[i].insert(ind_x, h_i[i]) # add height value + #print(z_sub[i]) + + dp = 0 + dA = 0 + eps = 1e-06 + for j in range(len(x_sub[i])-1): + if abs(z_sub[i][j+1] - h_i[i]) > eps or abs(z_sub[i][j] - h_i[i]) > eps: + dp += np.hypot(x_sub[i][j+1] - x_sub[i][j], + abs(z_sub[i][j+1] - z_sub[i][j])) + #print(dp) + # calculate area using trapezium rule + dA += (h_i[i] - (z_sub[i][j+1] + z_sub[i][j])/2)*(x_sub[i][j+1] - x_sub[i][j]) + #print('Area =', dA) + + p_i.append(dp) + A_i.append(dA) + + r_h.append(A_i[i]/p_i[i]) # ratio of area and wetted perimeter + #print('hydraulic radius =', r_h[i]) + K_i.append(A_i[i]*(1/n_co)*r_h[i]**(2/3)) # conveyance + Q_i.append(K_i[i]*slope**0.5) # discharge + + return p_i, A_i, r_h, h_i, K_i, Q_i + +def plot_region(xdata, labelx, + ydata1, labely1, + ydata2, labely2, + ydata3, labely3, titlep): + plt.xlabel(labelx) + plt.ylabel(labely1) + plt.title(titlep) + plt.plot(xdata, ydata1) + plt.show() + + plt.xlabel(labelx) + plt.ylabel(labely2) + plt.title(titlep) + plt.plot(xdata, ydata2) + plt.show() + + plt.xlabel(labelx) + plt.ylabel(labely3) + plt.title(titlep) + plt.plot(xdata, ydata3) + plt.show() + +def save_bc(): + with open(outputFilename, 'w') as f: + f.write('{:6} {:>2} {:>2} {:>10}\n'.format('#x', 'c', 'q', 'h')) + for ind_z, (xitem, zitem) in enumerate(zip(xin, zin)): + panel_x = bisect.bisect(markers, xitem) + if csa[ind_z] == 0: + # wall boundary condition + f.write('{:6.2f} {:>2}\n'.format(xitem, 2)) + elif panel_x == ind_p: + # imposed discharge within part-filled panel + f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( + xitem, 5, + -csa[ind_z]*panel_target_flow/csa_p[panel_x], + h_extra-zitem)) + else: + # imposed discharge within filled panels + f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( + xitem, 5, + -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x], + zmax-zitem)) + # read boundary definition file. definition_dict = read_definition('boundaryDefinition.txt') #for dd in definition_dict: @@ -126,75 +215,9 @@ print(zmin[1]) # channel overtopping height (minimum of left bank and right bank heights): zmax = min(zregion[panel[0]][0], zregion[panel[0]][-1]) - ztol -def conveyance(numH, n_co, xregion, zregion, zmin, zmax): - p_i = [] # wetted perimeter - A_i = [] # area - r_h = [] # hydraulic radius - h_i = [] # list of heights - K_i = [] # conveyance - Q_i = [] # discharge - x_sub = [[] for i in range(numH)] # list of x values in subregion - z_sub = [[] for i in range(numH)] # list of z values in subregion - for i in range(numH): - h_i.append(zmin + (i+1)*(zmax-zmin)/numH) - #print(zregion[zregion < h_i[i]]) - booleanArray = zregion < h_i[i] - #print(booleanArray[i]) - x_sub[i] += list(xregion[booleanArray]) - z_sub[i] += list(zregion[booleanArray]) - for interval in range(len(xregion)-1): - if booleanArray[interval+1] != booleanArray[interval]: - x_extra = xregion[interval] + (h_i[i]-zregion[interval])*(xregion[interval+1]-xregion[interval])/(zregion[interval+1]-zregion[interval]) - bisect.insort(x_sub[i], x_extra) # add intercept value - ind_x = x_sub[i].index(x_extra) - z_sub[i].insert(ind_x, h_i[i]) # add height value - #print(z_sub[i]) - - dp = 0 - dA = 0 - eps = 1e-06 - for j in range(len(x_sub[i])-1): - if abs(z_sub[i][j+1] - h_i[i]) > eps or abs(z_sub[i][j] - h_i[i]) > eps: - dp += np.hypot(x_sub[i][j+1] - x_sub[i][j], - abs(z_sub[i][j+1] - z_sub[i][j])) - #print(dp) - # calculate area using trapezium rule - dA += (h_i[i] - (z_sub[i][j+1] + z_sub[i][j])/2)*(x_sub[i][j+1] - x_sub[i][j]) - #print('Area =', dA) - - p_i.append(dp) - A_i.append(dA) - - r_h.append(A_i[i]/p_i[i]) # ratio of area and wetted perimeter - #print('hydraulic radius =', r_h[i]) - K_i.append(A_i[i]*(1/n_co)*r_h[i]**(2/3)) # conveyance - Q_i.append(K_i[i]*slope**0.5) # discharge - - return p_i, A_i, r_h, h_i, K_i, Q_i #print(h_i) -def plot_region(xdata, labelx, - ydata1, labely1, - ydata2, labely2, - ydata3, labely3, titlep): - plt.xlabel(labelx) - plt.ylabel(labely1) - plt.title(titlep) - plt.plot(xdata, ydata1) - plt.show() - - plt.xlabel(labelx) - plt.ylabel(labely2) - plt.title(titlep) - plt.plot(xdata, ydata2) - plt.show() - - plt.xlabel(labelx) - plt.ylabel(labely3) - plt.title(titlep) - plt.plot(xdata, ydata3) - plt.show() p_i = [[] for _ in range(num_panels)] A_i = [[] for _ in range(num_panels)] @@ -281,25 +304,5 @@ for i, p in enumerate(panel): #print('csa_p[0] = {} csa_p[1] = {}'.format(csa_p[0], csa_p[1])) #print('A_i_west = {} A_i = {} A_i_east = {}'.format(A_extra, A_i[-1], A_i_east[-1])) -def save_bc(): - with open(outputFilename, 'w') as f: - f.write('{:6} {:>2} {:>2} {:>10}\n'.format('#x', 'c', 'q', 'h')) - for ind_z, (xitem, zitem) in enumerate(zip(xin, zin)): - panel_x = bisect.bisect(markers, xitem) - if csa[ind_z] == 0: - # wall boundary condition - f.write('{:6.2f} {:>2}\n'.format(xitem, 2)) - elif panel_x == ind_p: - # imposed discharge within part-filled panel - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, 5, - -csa[ind_z]*panel_target_flow/csa_p[panel_x], - h_extra-zitem)) - else: - # imposed discharge within filled panels - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, 5, - -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x], - zmax-zitem)) save_bc() |