From 391d514d04c39806d659a6c2534470f222d041f4 Mon Sep 17 00:00:00 2001 From: Paul Garlick Date: Mon, 30 Sep 2019 16:30:16 +0100 Subject: rename makeBoundaryFile.py to makeBoundary.py; add boundary definition file. --- makeBoundaryFile.py | 208 ---------------------------------------------------- 1 file changed, 208 deletions(-) delete mode 100755 makeBoundaryFile.py (limited to 'makeBoundaryFile.py') diff --git a/makeBoundaryFile.py b/makeBoundaryFile.py deleted file mode 100755 index e3d30cb..0000000 --- a/makeBoundaryFile.py +++ /dev/null @@ -1,208 +0,0 @@ -#!/usr/bin/env python3 - -import matplotlib.pyplot as plt -import numpy as np -import bisect - -with open('./1D_top.txt', "r") as data: - xch, ych, zch = np.loadtxt(data, delimiter=' ', unpack=True) - -# Fit with polyfit -m, c = np.polyfit(ych, zch, 1) -print('gradient =', m, 'intercept =', c) - -slope = abs(m) # slope at top boundary -target_flow = 2.0 # imposed discharge - -with open('../topography/top_boundary.xyz', "r") as topo: - xin, yin, zin = np.loadtxt(topo, delimiter=' ', unpack=True) - -# array index and co-ordinates are related by: -# index = (co-ord - 0.25)*2 -xregion = xin[280:301] -zregion = zin[280:301] - -xregion_west = xin[100:281] -zregion_west = zin[100:281] - -xregion_east = xin[300:408] -zregion_east = zin[300:408] - -#print(zregion) - -zmin = zregion.min() # minimum height -zmax = zregion[-1]-0.01 # overtopping height -zmax_west = zmax -zmax_east = zmax - -zmin_west = zregion_west.min() -zmin_east = zregion_east.min() - -print(zmin_east) - -numH = 50 # number of height intervals -n_co_chan = 0.035 # Manning's coefficient for inland water -n_co_west = 0.040 # Manning's coefficient for general surface -n_co_east = 0.040 # Manning's coefficient for general surface -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, A_i, r_h, h_i, K_i, Q_i = conveyance( - numH, n_co_chan, xregion, zregion,zmin, zmax) - -plot_region(h_i-zmin, 'maximum depth / m', - r_h, 'hydraulic radius / m', - K_i, r'conveyance / $m^3/s$', - Q_i, r'discharge / $m^3/s$', 'main channel flow') - -p_i_west, A_i_west, r_h_west, h_i_west, K_i_west, Q_i_west = conveyance( - numH, n_co_west, xregion_west, zregion_west, zmin_west, zmax_west) - -plot_region(h_i_west-zmin_west, 'maximum depth / m', - r_h_west, 'hydraulic radius / m', - K_i_west, r'conveyance / $m^3/s$', - Q_i_west, r'discharge / $m^3/s$', 'west overland flow') - -p_i_east, A_i_east, r_h_east, h_i_east, K_i_east, Q_i_east = conveyance( - numH, n_co_east, xregion_east, zregion_east, zmin_east, zmax_east) - -plot_region(h_i_east-zmin_east, 'maximum depth / m', - r_h_east, 'hydraulic radius / m', - K_i_east, r'conveyance / $m^3/s$', - Q_i_east, r'discharge / $m^3/s$', 'east overland flow') - -target_flow_west = target_flow - Q_i[-1] - Q_i_east[-1] -# calculate velocity: note dependence on hydraulic radius -velocity_channel = Q_i[-1]/A_i[-1] -velocity_east = Q_i_east[-1]/A_i_east[-1] - -print(target_flow_west) - -# find insertion point for target flow value -ind_q = bisect.bisect(Q_i_west, target_flow_west) -print(ind_q) - -# find height at target flow by linear interpolation -h_extra = h_i_west[ind_q-1] + (h_i_west[ind_q]-h_i_west[ind_q-1])*(target_flow_west-Q_i_west[ind_q-1])/(Q_i_west[ind_q]-Q_i_west[ind_q-1]) -print(h_i_west[ind_q-1], h_extra, h_i_west[ind_q]) -# find area at target flow by linear interpolation -A_extra = A_i_west[ind_q-1] + (h_extra-h_i_west[ind_q-1])*(A_i_west[ind_q]-A_i_west[ind_q-1])/(h_i_west[ind_q]-h_i_west[ind_q-1]) -print(r_h_west[ind_q-1], r_h_west[ind_q]) - -velocity_west = target_flow_west/A_extra -print(velocity_channel, velocity_east, velocity_west) - -dX = (xin[0]+xin[-1])/len(xin) # cell size -print('dX =', dX) -csa = np.zeros(len(xin)) # cross-sectional area -csa_west = 0 -csa_chan = 0 -csa_east = 0 -for index, xitem in enumerate(xin): - if xitem < xin[280]: - csa[index] = max(0, (h_extra - zin[index])*dX) - csa_west += csa[index] - elif xitem < xin[301]: - csa[index] = max(0, (zmax - zin[index])*dX) - csa_chan += csa[index] - else: - csa[index] = max(0, (zmax_east - zin[index])*dX) - csa_east += csa[index] - -#print('csa_west = {} csa_chan = {} csa_east = {}'.format(csa_west, csa_chan, csa_east)) -#print('A_i_west = {} A_i = {} A_i_east = {}'.format(A_extra, A_i[-1], A_i_east[-1])) - -def save_bc(): - with open('BCTop.txt', '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)): - if csa[ind_z] == 0: - # wall boundary condition - f.write('{:6.2f} {:>2}\n'.format(xitem, 2)) - elif xitem < xin[280]: - # imposed discharge on western side - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, 5, - -csa[ind_z]*target_flow_west/csa_west, - h_extra-zitem)) - elif xitem < xin[301]: - # imposed discharge within channel - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, 5, - -csa[ind_z]*Q_i[-1]/csa_chan, - zmax-zitem)) - else: - # imposed discharge on eastern side - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, 5, - -csa[ind_z]*Q_i_east[-1]/csa_west, - zmax_east-zitem)) - -save_bc() -- cgit