diff options
Diffstat (limited to 'makeBoundary')
| -rwxr-xr-x | makeBoundary | 364 |
1 files changed, 0 insertions, 364 deletions
diff --git a/makeBoundary b/makeBoundary deleted file mode 100755 index a0d3773..0000000 --- a/makeBoundary +++ /dev/null @@ -1,364 +0,0 @@ -#!/usr/bin/env python3 - -import matplotlib.pyplot as plt -import numpy as np -import bisect -import ast -import math -import argparse -import sys - -#TODO: use YAML/ruamel.yaml for configuration file. -def read_definition(filename): - ddict = {} - with open(filename, "r") as f: - for line in f: - items = line.split(': ', 1) - if len(items) == 2: - 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(outputfile): - with open(outputfile, '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 == panel[ind_p]: - # imposed discharge within part-filled panel - f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format( - xitem, btype, - -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, btype, - -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x], - zmax-zitem)) - -def interp(extra2, max1, min1, max2, min2): - # use similar triangles to perform linear interpolation - extra1 = min1 + (max1 - min1)*(extra2 - min2)/(max2 - min2) - - return extra1 - -# read command line argument: -parser = argparse.ArgumentParser( - description="generate FullSWOF boundary files") -parser.add_argument("location", help="boundary location") -args = parser.parse_args() - -if args.location == 'top': - inputFilename = "boundaryTop.txt" - outputFilename = "BCTop.txt" -elif args.location == 'bottom': - inputFilename = "boundaryBottom.txt" - outputFilename = "BCBottom.txt" -elif args.location == 'left': - inputFilename = "boundaryLeft.txt" - outputFilename = "BCLeft.txt" -elif args.location == 'right': - inputFilename = "boundaryRight.txt" - outputFilename = "BCRight.txt" - -# read boundary definition file: -definition_dict = read_definition(inputFilename) -#for dd in definition_dict: -# print(definition_dict[dd]) -btype = definition_dict["type"] # boundary type (1--5) -slope = abs(definition_dict["slope"]) # slope at top boundary -target_flow = definition_dict["target_flow"] # imposed discharge -plotting = definition_dict["plotting"] # enable or disable plotting -n_co = definition_dict["n_co"] # Manning's 'n' coefficients -# TODO: use weighted mean 'n' values. See -# http://help.floodmodeller.com/isis/ISIS/River_Section.htm (Eq. 4) -# Note: weighted mean calculation requires roughness map. -markers = definition_dict["markers"] # distances from corner point -panel = definition_dict["panel"] # panel fill order -ztol = definition_dict["ztol"] # tolerance in overtopping height -numH = definition_dict["numH"] # number of height intervals - - -# print(len(markers)) - -# 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) - -# read topography: -with open("./topography.txt", "r") as topo: - xtp, ytp, ztp = np.loadtxt(topo, delimiter=' ', unpack=True) - - -xmax = (xtp[0]+xtp[-1]) # domain extent in x-direction -ymax = (ytp[0]+ytp[-1]) # domain extent in y-direction -ncols = int(math.sqrt(len(xtp)*xmax/ymax)) # number of cells in x-direction -nrows = int(len(xtp)/ncols) # number of cells in y-direction -dX = xmax/ncols # cell size -print('dX =', dX) - -#print(ncols, nrows) - -# extract slices from height data array. Note: xyz format uses ncols -# blocks, with nrows lines per block. -if args.location == 'top': - xin = xtp[nrows-1:len(xtp):nrows] - yin = ytp[nrows-1:len(xtp):nrows] - zin = ztp[nrows-1:len(xtp):nrows] -elif args.location == 'bottom': - xin = xtp[0:len(xtp):nrows] - yin = ytp[0:len(xtp):nrows] - zin = ztp[0:len(xtp):nrows] -elif args.location == 'left': - xin = xtp[:nrows] - yin = ytp[:nrows] - zin = ztp[:nrows] -elif args.location == 'right': - xin = xtp[nrows*(ncols-1):] - yin = ytp[nrows*(ncols-1):] - zin = ztp[nrows*(ncols-1):] - - -# print(xin) - -num_panels = len(panel) # number of panels across boundary - -# convert marker co-ordinates to array indices: -marker_ind = [0] -for i in range(len(markers)): - marker_ind.append(int(markers[i]/dX)) -if args.location == 'left' or args.location == 'right': - marker_ind.append(nrows) -elif args.location == 'top' or args.location == 'bottom': - marker_ind.append(ncols) - - -# print(marker_ind) - -xregion = [] -zregion = [] -zmin = [] -for p in range(num_panels): - # identify regions: - xregion.append(xin[marker_ind[p]:marker_ind[p+1]]) - zregion.append(zin[marker_ind[p]:marker_ind[p+1]]) - # identify minimum heights within each panel: - zmin.append(zregion[p].min()) - -# xregion_west = xin[100:281] -# zregion_west = zin[100:281] - -# xregion_east = xin[300:408] -# zregion_east = zin[300:408] - -# print(zregion) -# print(xin[12:20]) - -print('zmin =', zmin) - -# channel overtopping height (minimum of left bank and right bank heights): -zmax = min(zregion[panel[0]][0], zregion[panel[0]][-1]) - ztol - -print('zmax =', zmax) - -#print(h_i) - - - -p_i = [[] for _ in range(num_panels)] -A_i = [[] for _ in range(num_panels)] -r_h = [[] for _ in range(num_panels)] -h_i = [[] for _ in range(num_panels)] -K_i = [[] for _ in range(num_panels)] -Q_i = [[] for _ in range(num_panels)] -for p in range(num_panels): - if p == panel[0]-1 and zregion[p][-1] < zmax: - # ensure end node in region to the left of channel is dry: - xregion[p] = np.append(xregion[p], xin[marker_ind[p]]) - zregion[p] = np.append(zregion[p], zin[marker_ind[p]]) - if p == panel[0]+1 and zregion[p][0] < zmax: - # ensure start node in region to the right of channel is dry: - xregion[p] = np.insert(xregion[p], 0, xin[marker_ind[p]-1]) - zregion[p] = np.insert(zregion[p], 0, zin[marker_ind[p]-1]) - if zmax > zmin[p]: - p_i[p], A_i[p], r_h[p], h_i[p], K_i[p], Q_i[p] = conveyance( - numH, - n_co[p], - xregion[p], - zregion[p], - zmin[p], - zmax) - if plotting: - plot_region( - h_i[p]-zmin[p], 'maximum depth / m', - r_h[p], 'hydraulic radius / m', - K_i[p], r'conveyance / $m^3/s$', - Q_i[p], r'discharge / $m^3/s$', - 'Panel {}'.format(p)) - else: - p_i[p], A_i[p], r_h[p], h_i[p], K_i[p], Q_i[p] = [ - [0] * numH for _ in range(6)] - -# sort list of discharge lists according to panel fill order: -sortedQ = [Q_i[i] for i in panel] -# create cumulative discharge list: -total_flow = np.cumsum([item[-1] for item in sortedQ]) -print('total_flow = ', total_flow) -# target_flow_west = target_flow - Q_i[-1] - Q_i_east[-1] -# calculate velocity: note dependence on hydraulic radius -velocity_channel = Q_i[panel[0]][-1]/A_i[panel[0]][-1] -# velocity_east = Q_i_east[-1]/A_i_east[-1] - -# print(target_flow_west) -# find part-filled panel: -if total_flow[-1] > target_flow: - ind_p = bisect.bisect(total_flow, target_flow) -else: - print('Error: imposed discharge is higher than total capacity of panels.') - sys.exit() - -print('index of part-filled panel:', ind_p) - -# calculate target flow in part-filled panel: -if ind_p == 0: - panel_target_flow = target_flow -else: - panel_target_flow = target_flow - total_flow[ind_p-1] - -# find insertion point for target flow value: -ind_q = bisect.bisect(Q_i[panel[ind_p]], panel_target_flow) - -print('insertion point =', ind_q) - -# find height at target flow by linear interpolation -h_extra = interp( - panel_target_flow, - h_i[panel[ind_p]][ind_q], - h_i[panel[ind_p]][ind_q-1], - Q_i[panel[ind_p]][ind_q], - Q_i[panel[ind_p]][ind_q-1]) - -print('heights:', h_i[panel[ind_p]][ind_q-1], h_extra, h_i[panel[ind_p]][ind_q]) - -# find area at target flow by linear interpolation -A_extra = interp( - h_extra, - A_i[panel[ind_p]][ind_q], - A_i[panel[ind_p]][ind_q-1], - h_i[panel[ind_p]][ind_q], - h_i[panel[ind_p]][ind_q-1]) - -print('hydraulic radii:', r_h[panel[ind_p]][ind_q-1], r_h[panel[ind_p]][ind_q]) - -velocity_panel = panel_target_flow/A_extra -print('velocities:', velocity_channel, velocity_panel) - -csa = np.zeros(len(xin)) # cross-sectional area of element -csa_p = np.zeros(num_panels) # cross-sectional area of panel -for i, p in enumerate(panel): - if i < ind_p: - # panels are filled - area_sum = 0 - for m in range(marker_ind[p], marker_ind[p+1]): - csa[m] = max(0, (zmax - zin[m])*dX) - area_sum += csa[m] - csa_p[p] = area_sum - elif i == ind_p: - # panel is part-filled - area_sum = 0 - for m in range(marker_ind[p], marker_ind[p+1]): - csa[m] = max(0, (h_extra - zin[m])*dX) - area_sum += csa[m] - csa_p[p] = area_sum - else: - # panel is empty - for m in range(marker_ind[p], marker_ind[p+1]): - csa[m] = 0 - csa_p[p] = 0 - - -#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])) - - -save_bc(outputFilename) |