rename makeBoundaryFile.py to makeBoundary.py; add boundary definition file.

1 files changed, 0 insertions, 208 deletions

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()