aboutsummaryrefslogtreecommitdiff
path: root/makeBoundary
diff options
context:
space:
mode:
Diffstat (limited to 'makeBoundary')
-rwxr-xr-xmakeBoundary364
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)