path: root/makeBoundary.py

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
import bisect
import ast
import math

#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

# read boundary definition file.
definition_dict = read_definition('boundaryDefinition.txt')
#for dd in definition_dict:
#    print(definition_dict[dd])
slope       = abs(definition_dict["slope"])  # slope at top boundary
target_flow = definition_dict["target_flow"] # imposed discharge
location    = definition_dict["location"]    # boundary location
n_co_chan   = definition_dict["n_co_chan"]   # coefficient for inland water
n_co_west   = definition_dict["n_co_west"]   # coefficient for general surface
n_co_east   = definition_dict["n_co_east"]   # coefficient for general surface
# TODO: use weighted mean 'n' value.  See http://help.floodmodeller.com/isis/ISIS/River_Section.htm (Eq. 4)
# Note: weighted mean calculation requires roughness map.
height_data = definition_dict["height_data"] # topography
markers     = definition_dict["markers"]     # distances from corner point
channel     = definition_dict["channel"]     # identifier of channel panel
ztol        = definition_dict["ztol"]        # tolerance in overtopping height

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

with open(height_data, "r") as topo:
    xtp, ytp, ztp = np.loadtxt(topo, delimiter=' ', unpack=True)

# domain extent in x-direction:
xmax = (xtp[0]+xtp[-1])
# domain extent in y-direction:
ymax = (ytp[0]+ytp[-1])
# number of cells in x-direction:
ncols = int(math.sqrt(len(xtp)*xmax/ymax))
# number of cells in y-direction:
nrows = int(len(xtp)/ncols)
# cell size
dX = xmax/ncols
print('dX =', dX)

#print(ncols, nrows)

# extract slices from height data array.  Note: xyz format uses ncols
# blocks, with nrows lines per block.
if location == 'top':
    xin = xtp[nrows-1:len(xtp):nrows]
    yin = ytp[nrows-1:len(xtp):nrows]
    zin = ztp[nrows-1:len(xtp):nrows]
elif location == 'bottom':
    xin = xtp[0:len(xtp):nrows]
    yin = ytp[0:len(xtp):nrows]
    zin = ztp[0:len(xtp):nrows]
elif location == 'left':
    xin = xtp[:nrows]
    yin = ytp[:nrows]
    zin = ztp[:nrows]
elif location == 'right':
    xin = xtp[nrows*(ncols-1):]
    yin = ytp[nrows*(ncols-1):]
    zin = ztp[nrows*(ncols-1):]
    

print(xin)

num_panels = len(markers) + 1 # number of panels across boundary

# convert panel co-ordinates to array indices:
panel_ind = [0]
for i in range(len(markers)):
    panel_ind.append(int(markers[i]/dX - 1/2))
if location == 'left' or location == 'right':
    panel_ind.append(nrows-1)
elif location == 'top' or location == 'bottom':
    panel_ind.append(ncols-1)


# print(panel_ind)

xregion = []
zregion = []
for p in range(num_panels):
    xregion.append(xin[panel_ind[p]:panel_ind[p+1]])
    zregion.append(zin[panel_ind[p]:panel_ind[p+1]])

# 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])

# minimum height in channel:
zmin = zregion[channel].min()
# overtopping height (minimum of left bank and right bank height):
zmax = min(zregion[channel][0], zregion[channel][-1]) - ztol

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

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)

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