Uniform rain on an idealised catchment¶

This 2D rainfall runoff test is based on the paper of (Cea, 2008) (https://iwra.org/congress/2008/resource/authors/abs478_article.pdf), which compare a numerical model against a laboratory experiment (Iwagaki,1955). This configuration consists in an ideal catchment with three land surface of constant slope, and two walls increasing the travel time of run-off. A constant rain is applied over the domain, associated with constant friction. The outlet hydrograph is used for comparing models results against experiments mesurements.

Model set-up¶

Python set-up¶

Installation¶

The following libraries will be used in this tutorial. They can be installed using the following requirement file (pip install -r requierements.txt), or any of your favorite method.

matplotlib-3.4.2 scipy-1.6.3 netCDF4-1.5.6 numpy-1.20.3 xarray-0.18.2 python-3.9.5

Loading libraries¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Rectangle
from scipy import interpolate
from scipy.interpolate import LinearNDInterpolator

import xarray as xr
import os
import shutil
from netCDF4 import Dataset
import netCDF4 as nc4
import utm

## Read text file given its name (and path), number of columns, number of header lines
def readCSV_Ncol(name,n_col,n_pass):
    X=[]
    for i in range(n_col):
        X.append([])

    f=open(name,'r')
    n_line=0
    for line in f.readlines():
        if n_line<n_pass:
            n_line = n_line+1
        else:
            sline=line.split(',')
            for i in range(n_col):
                X[i].append(float(sline[i]))
    f.close()
    return(X)

import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap from matplotlib.patches import Rectangle from scipy import interpolate from scipy.interpolate import LinearNDInterpolator import xarray as xr import os import shutil from netCDF4 import Dataset import netCDF4 as nc4 import utm ## Read text file given its name (and path), number of columns, number of header lines def readCSV_Ncol(name,n_col,n_pass): X=[] for i in range(n_col): X.append([]) f=open(name,'r') n_line=0 for line in f.readlines(): if n_line

Setting up the directories¶

In [ ]:

# Path to the folder containing the code executable
codeDir = r".\CodeDir"
# Creation (if needed) and access to the working directory
if not os.path.isdir ('auto_CEA2008'):
    os.mkdir('auto_CEA2008')
os.chdir('auto_CEA2008')
MyDir = os.getcwd()

# Path to the folder containing the code executable codeDir = r".\CodeDir" # Creation (if needed) and access to the working directory if not os.path.isdir ('auto_CEA2008'): os.mkdir('auto_CEA2008') os.chdir('auto_CEA2008') MyDir = os.getcwd()

Preparation of the bathymetry¶

The bathymetry is given in an ascii file (Download topo.asc), that can be directly read by BG_Flood. The initial model has been surronded by walls and a spillway has been added at the outlet to ensure a minimum water elevation in the inital domain (outlet is located at the top of the model).

In [4]:

grid=np.loadtxt('topo.asc',skiprows=6)
x=np.linspace(0, len(grid[1,:])/100, len(grid[1,:])+1)
y=np.linspace(0, len(grid[:,1])/100, len(grid[:,1])+1)
plt.pcolor(x,y,grid, vmin=-0.1, vmax=0.1)
plt.colorbar()

grid=np.loadtxt('topo.asc',skiprows=6) x=np.linspace(0, len(grid[1,:])/100, len(grid[1,:])+1) y=np.linspace(0, len(grid[:,1])/100, len(grid[:,1])+1) plt.pcolor(x,y,grid, vmin=-0.1, vmax=0.1) plt.colorbar()

Out[4]:

<matplotlib.colorbar.Colorbar at 0x294817b7d30>

Creation of the forcing¶

In this configuration, the code will be forced by an homogeneous rain during 45s then, switch off.

Writing the rain forcing file¶

The boundary condition file will be defined by a time serie with a column for time, and one for millimeters of rain: rainfall intensity $I = 317.0 mm/h$).

In [ ]:

Tr=[0.0,45.0,45.1,150.0]
Hr=[317.0,317.0,0.0,0.0]
filename = "cstRain_forcing.txt"
f=open(filename,'w')
for i in range(len(Tr)):
    f.write('%f \t' % Tr[i])
    f.write('%f \n' % Hr[i])
f.close()

plt.plot(Tr,Hr,'.-')
plt.title('Rain forcing')
plt.ylabel('I (mm/h)')
plt.xlabel('t(s)')
plt.savefig('rain-forcing.png')
plt.show()

Tr=[0.0,45.0,45.1,150.0] Hr=[317.0,317.0,0.0,0.0] filename = "cstRain_forcing.txt" f=open(filename,'w') for i in range(len(Tr)): f.write('%f \t' % Tr[i]) f.write('%f \n' % Hr[i]) f.close() plt.plot(Tr,Hr,'.-') plt.title('Rain forcing') plt.ylabel('I (mm/h)') plt.xlabel('t(s)') plt.savefig('rain-forcing.png') plt.show()

Writing the BG_Params file¶

This text file is the main interface to BG_Flood solver. It contains a list of key words and associated values, to modify the default topographic grid and time keeping parameters, boundary and external forcing, as well as output options and parameter values. Different keys can be accepted for the same parameter. (See List of inputs parameters for a detailled list of the input parameters). The minimum input is the bathymetry file. In this example, we will conserve a simple version of the BG_param.txt file where only the modified values will be listed (default values will be kept for the parameters not mentionned):

In [6]:

key=[]
value=[]

key=[] value=[]

• name of the bathymetry file (variable name associated to altitude is also needed the format is netCDF ("filename?variable"))

In [7]:

key.append('topofile')
value.append('topo.asc')

key.append('topofile') value.append('topo.asc')

• mask to remove parts of the domain with elevation over 0.1m

In [8]:

key.append('mask')
value.append('0.1')

key.append('mask') value.append('0.1')

• Timekeeping parameters
  • "inittime" : Start time for the simulation
  • "endtime" : End of simulation (in s)

In [ ]:

key.append('inittime')
value.append(Tr[0])
key.append('endtime')
value.append(Tr[-1])

key.append('inittime') value.append(Tr[0]) key.append('endtime') value.append(Tr[-1])

• Resolution of the simulation (by default, the resolution of the bathymetry or topography file is used)

In [10]:

key.append('dx')
value.append('0.005')

key.append('dx') value.append('0.005')

• Choice of the GPU/CPU mode (-1 if no GPU else 0 for GPU, by default if present on the computer) and precison of the model

In [11]:

#key.append('gpudevice')
#value.append(0)
key.append('doubleprecision')
value.append('1')

#key.append('gpudevice') #value.append(0) key.append('doubleprecision') value.append('1')

• Flow and models parameters
  • "frictionmodel": bottom friction model (-1: Manning model, 0:, 1: Smart model (default))
  • "cf": friction coefficient
  • "eps": minimum water elevation
  • "theta": minmod limiter that can be used to tune the momentum dissipation (theta=1 gives minmod, the most dissipative limiter and theta = 2 gives superbee, the least dissipative) (by default 1.3).

In [12]:

key.append('frictionmdel')
value.append('1')
key.append('cf')
value.append('0.000004')
key.append('eps')
value.append('0.00000100')
key.append('theta')
value.append('1.14')

key.append('frictionmdel') value.append('1') key.append('cf') value.append('0.000004') key.append('eps') value.append('0.00000100') key.append('theta') value.append('1.14')

• Boundaries condition

Boundaries are refered by their position, using "top/bottom/right/left" keywords. They are associated to a boundary type ( 0:wall; 1: Neumann (Default); 2:Dirichlet (zs); 3: abs1d) and possibly a file containing a time serie. In this case, the file name is placed before the type, coma-separated. Here, the by default Neumann type boundary conditions are conserved.

• Exterior forcings

Different forcings are available (river, rain, atmospherique pressure...). In this example, we will use the rain forcing. This can be a time serie in a text file if the forcing is uniform on all the domain (like in this case) or a netcdf file if the forcing is space dependent (and time dependent).

In [13]:

key.append('rainfile')
value.append('cstRain_forcing.txt')

key.append('rainfile') value.append('cstRain_forcing.txt')

• Initial conditions
  • "zinit" : initial water level

In [14]:

key.append('zsinit')
value.append('-0.5')

key.append('zsinit') value.append('-0.5')

• Outputs:

  The code can output 2D map of 2D time-varying fields in netCDF format and time series at a point location (not used here).

  • "outputtimestep" : Time step for 2D fields outputs
  • "smallnc": to scaled and saved 2D fields as a short integer if 1 (default value)
  • "outvars": to choose the variable to output (by default: $hh$, $uu$, $vv$, $zb$, $zs$)

In [15]:

key.append('outputtimestep')
value.append('1.0')
key.append('smallnc')
value.append('0')
key.append('outvars')
value.append('zs,h,u,v,zb,hmax,zsmax')

key.append('outputtimestep') value.append('1.0') key.append('smallnc') value.append('0') key.append('outvars') value.append('zs,h,u,v,zb,hmax,zsmax')

In [16]:

key

key

Out[16]:

['topofile',
 'mask',
 'totaltime',
 'endtime',
 'dx',
 'doubleprecision',
 'frictionmdel',
 'cf',
 'eps',
 'theta',
 'rainfile',
 'zsinit',
 'outputtimestep',
 'smallnc',
 'outvars']

In [17]:

value

value

Out[17]:

['topo.asc',
 '0.1',
 0.0,
 60.0,
 '0.005',
 '1',
 '1',
 '0.000004',
 '0.00000100',
 '1.14',
 'cstRain_forcing.txt',
 '-0.5',
 '1.0',
 '0',
 'zs,h,u,v,zb,hmax,zsmax']

These data are then saved in a BG_Param.txt file.

In [18]:

data = np.column_stack([key, value])
datafile_path = "BG_param.txt"
np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=')

data = np.column_stack([key, value]) datafile_path = "BG_param.txt" np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=')

Launching the code BG_flood¶

All the needed links to the libraries (.dll files) and the BG_flood executable are copied in the working folder.

In [43]:

myFiles=[f for f in os.listdir(codeDir) if f.endswith(".dll")]
myFiles.append('BG_flood.exe')
for f in myFiles:
    origin=os.path.join(codeDir,f)
    destination=os.path.join(MyDir,f)
    shutil.copy(origin, destination) 

myFiles=[f for f in os.listdir(codeDir) if f.endswith(".dll")] myFiles.append('BG_flood.exe') for f in myFiles: origin=os.path.join(codeDir,f) destination=os.path.join(MyDir,f) shutil.copy(origin, destination)

Launch of the code¶

The code can be launched by a double clic on the executable or from python using the following line:

In [ ]:

os.system('BG_flood.exe >> BG_out.txt')

os.system('BG_flood.exe >> BG_out.txt')

Comparaison of the data output with the experimental results¶

After identifying the position of the section at the end of the slope ($y=0.1$), we will calculate the flux through this section summing, through all the cells ($C_s$) at this location, the water elevation ($h$) times the velocity normal to the section ($v$):

$$ F_{out}(t) = \sum \limits _{C_s} h \, v \, dx$$

This will be calculated at each steps where the 2D fields have been outputted in the netCDF file.

Reading the benchmark data¶

This data (Download reference data) has been extracted from the initial paper (Iwagaki, 1955).

In [21]:

[xC1,yC1,xC3,yC3]=readCSV_Ncol('CEA_C1-C3_datasets.csv',4,2)

[xC1,yC1,xC3,yC3]=readCSV_Ncol('CEA_C1-C3_datasets.csv',4,2)

Reading the BG_flood outputs¶

Here, the model is uniform (same grid size everywhere). The handling of BG_Flood results is then relatively easy. When an adaptive grid with different level of refinements is used, different variables, ex: hh_P0, hh_P1, hh_N2, will be saved in the output file. Each corresponds to a layer of refinement and have different dimensions, XX_P0 corresponding to the reference level, XX_P1 corresponding to one level up (+1), XX_N2 corresponding to 2 levels down (-2). One level up correspond to dividing the resolution by two.

In [31]:

fn = 'Output_50.nc'
ds = nc4.Dataset(fn)
ds

fn = 'Output_50.nc' ds = nc4.Dataset(fn) ds

Out[31]:

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    maxlevel: 0
    minlevel: 0
    xmin: -0.005
    xmax: 2.075
    ymin: -0.005
    ymax: 2.555
    dimensions(sizes): time(61), blockid(832), xx_P0(416), yy_P0(512)
    variables(dimensions): float32 time(time), int32 blockid(blockid), float32 blockxo(blockid), float32 blockyo(blockid), float32 blockwidth(blockid), int32 blocklevel(blockid), int32 blockstatus(blockid), float64 xx_P0(xx_P0), float64 yy_P0(yy_P0), float32 zb_P0(time, yy_P0, xx_P0), float32 zs_P0(time, yy_P0, xx_P0), float32 u_P0(time, yy_P0, xx_P0), float32 v_P0(time, yy_P0, xx_P0), float32 h_P0(time, yy_P0, xx_P0)
    groups: 

Calculating the flux at the output¶

We calculate the flux through a $1m$-long section ($x \in [0.5,1.5]$)

In [33]:

xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zs=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]

F=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dy
    F.append(-ft)
                    
plt.plot(tt,F,'.-')
plt.title('Flux at the output of the domain')
plt.ylabel('Flux ($m^3$/s)')
plt.xlabel('Time (s)')

xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zs=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] F=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

Out[33]:

Text(0.5, 0, 'Time (s)')

In [34]:

plt.pcolormesh(xx,yy,hh[50,:,:]-hh[0,:,:], shading='auto')
#plt.plot(xx[200],yy[50],'o')
#plt.plot(xx[200],yy[250],'o')
#plt.plot(xx[200],yy[400],'o')
plt.title('Water elevation in the idealized catchment')
plt.colorbar()

plt.pcolormesh(xx,yy,hh[50,:,:]-hh[0,:,:], shading='auto') #plt.plot(xx[200],yy[50],'o') #plt.plot(xx[200],yy[250],'o') #plt.plot(xx[200],yy[400],'o') plt.title('Water elevation in the idealized catchment') plt.colorbar()

Out[34]:

<matplotlib.colorbar.Colorbar at 0x2610a346a30>

In [9]:

plt.pcolormesh(xx,yy,hh[50,:,:], shading='auto', vmax=0.0)
plt.plot(xx[200],yy[50],'o')
plt.plot(xx[200],yy[250],'o')
plt.plot(xx[200],yy[400],'o')
plt.colorbar()

plt.pcolormesh(xx,yy,hh[50,:,:], shading='auto', vmax=0.0) plt.plot(xx[200],yy[50],'o') plt.plot(xx[200],yy[250],'o') plt.plot(xx[200],yy[400],'o') plt.colorbar()

Out[9]:

<matplotlib.colorbar.Colorbar at 0x294804e0d90>

In [10]:

plt.plot(tt, hh[:,50,200],'o')
plt.plot(tt, hh[:,250,200],'d')
plt.plot(tt, hh[:,400,200])

plt.plot(tt, hh[:,50,200],'o') plt.plot(tt, hh[:,250,200],'d') plt.plot(tt, hh[:,400,200])

Out[10]:

[<matplotlib.lines.Line2D at 0x29483215eb0>]

In [11]:

ax=plt.subplot(1,1,1)
P1,=ax.plot(tt, zb[:,250,200],'o', label='zb')
P2,=ax.plot(tt, zs[:,250,200],'d', label='zs')
plt.legend(handles=[P1, P2])
ax.set_ylabel('z(m)')
plt.xlabel('time(s)')
plt.title('Time evolution of the surface elevation compare to the bottom elevation');

ax=plt.subplot(1,1,1) P1,=ax.plot(tt, zb[:,250,200],'o', label='zb') P2,=ax.plot(tt, zs[:,250,200],'d', label='zs') plt.legend(handles=[P1, P2]) ax.set_ylabel('z(m)') plt.xlabel('time(s)') plt.title('Time evolution of the surface elevation compare to the bottom elevation');

Test of the new Buttinger scheme¶

The simulatuion has been restarted (from the previous simulation) from the end of the rain period (45s) for the old Kurgonov and new Buttinger scheme for surface reconstruction.

In [12]:

fn = 'Output_55.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsk=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]

fn = 'Output_55.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsk=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0]

In [13]:

fn = 'Output_56.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsb=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]

ax=plt.subplot(1,1,1)
P1,=ax.plot(tt, zb[:,250,200],'o', label='zb')
P2,=ax.plot(tt, zsk[:,250,200],'s', label='zs using Kurgonov')
P3,=ax.plot(tt, zsb[:,250,200],'d', label='zs using Buttinger')
plt.legend(handles=[P1, P2, P3])
ax.set_ylabel('z(m)')
plt.xlabel('time(s)')
plt.title('Time evolution of the surface elevation compare to the bottom elevation');

fn = 'Output_56.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsb=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] ax=plt.subplot(1,1,1) P1,=ax.plot(tt, zb[:,250,200],'o', label='zb') P2,=ax.plot(tt, zsk[:,250,200],'s', label='zs using Kurgonov') P3,=ax.plot(tt, zsb[:,250,200],'d', label='zs using Buttinger') plt.legend(handles=[P1, P2, P3]) ax.set_ylabel('z(m)') plt.xlabel('time(s)') plt.title('Time evolution of the surface elevation compare to the bottom elevation');

Adding fix in interpolation of the forcing¶

A small error of typo on the identification of the running time, causing an error during the interpolation, was creating a negative rain (causing the observed error of negative water elevation).

In [14]:

fn = 'Output_59.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]

ax=plt.subplot(1,1,1)
P1,=ax.plot(tt, zb[:,250,200],'o', label='zb')
P2,=ax.plot(tt, zsk[:,250,200],'s', label='zs using Kurgonov')
P3,=ax.plot(tt, zsb[:,250,200],'d', label='zs using Buttinger')
P4,=ax.plot(tt, zsf[:,250,200],'.', label='zs using Buttinger after fix')
plt.legend(handles=[P1, P2, P3, P4])
ax.set_ylabel('z(m)')
plt.xlabel('time(s)')
plt.title('Time evolution of the surface elevation compare to the bottom elevation');

fn = 'Output_59.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] ax=plt.subplot(1,1,1) P1,=ax.plot(tt, zb[:,250,200],'o', label='zb') P2,=ax.plot(tt, zsk[:,250,200],'s', label='zs using Kurgonov') P3,=ax.plot(tt, zsb[:,250,200],'d', label='zs using Buttinger') P4,=ax.plot(tt, zsf[:,250,200],'.', label='zs using Buttinger after fix') plt.legend(handles=[P1, P2, P3, P4]) ax.set_ylabel('z(m)') plt.xlabel('time(s)') plt.title('Time evolution of the surface elevation compare to the bottom elevation');

Comparison between the Kurgonov and the Buttinger model¶

The initial implementation of the water surface reconstruction was based on Kurganov&Petrova2007. It is compared to a new surface reconstruction method proposed by Buttinger-Kreuzhuber-et-al2019.

In [15]:

fn = 'Output_ButtingerC1.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_b=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_b.append(-ft)

fn = 'Output_ButtingerC1.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_b=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [20]:

fn = 'Output_KurgonovC1.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_k=[]
for k in range(len(tt)):
    ft=0;
    LL=0
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    LL=LL+dx
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_k.append(-ft)

fn = 'Output_KurgonovC1.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_k=[] for k in range(len(tt)): ft=0; LL=0 for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [17]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P2,=ax1.plot(tt, F_k,'.-', label='zs using Kurganov')
P3,=ax1.plot(tt, F_b,'.-', label='zs using Buttinger')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
plt.legend(handles=handles)
plt.xlabel('Time(s)')
plt.title('Time evolution of the flux at the output of the domain, C1');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P2,=ax1.plot(tt, F_k,'.-', label='zs using Kurganov') P3,=ax1.plot(tt, F_b,'.-', label='zs using Buttinger') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) plt.legend(handles=handles) plt.xlabel('Time(s)') plt.title('Time evolution of the flux at the output of the domain, C1');

Test of the Case 3 of CEA¶

For this configuration, 328mm/h rain is applied two times 25s, with a 7s cut-off in between.

Writing the rain forcing file¶

The boundary condition file will be defined by a time serie with a column for time, and on for millimeters of rain.

In [23]:

Tr3=[0.0,25.0,25.0001,31.99999,32.0,57.0,57.1,150.0]
Hr3=[328.0,328.0,0.0,0.0,328.0,328.0,0.0,0.0]
filename = "cstRain_forcingC3.txt"
f=open(filename,'w')
for i in range(len(Tr3)):
    f.write('%f \t' % Tr3[i])
    f.write('%f \n' % Hr3[i])
f.close()

plt.plot(Tr3,Hr3,'.-')
plt.title('Rain forcing for test case 3')
plt.ylabel('h(mm)')
plt.xlabel('t(s)')
plt.savefig('rain-forcing_C3.png')
plt.show()

Tr3=[0.0,25.0,25.0001,31.99999,32.0,57.0,57.1,150.0] Hr3=[328.0,328.0,0.0,0.0,328.0,328.0,0.0,0.0] filename = "cstRain_forcingC3.txt" f=open(filename,'w') for i in range(len(Tr3)): f.write('%f \t' % Tr3[i]) f.write('%f \n' % Hr3[i]) f.close() plt.plot(Tr3,Hr3,'.-') plt.title('Rain forcing for test case 3') plt.ylabel('h(mm)') plt.xlabel('t(s)') plt.savefig('rain-forcing_C3.png') plt.show()

Modification of the Parameter files and running the code¶

In [ ]:

#For Kurgonov
    #Modify the BG_param file
    value[key.index('outfile')] = 'Output_KurgonovC3.nc'
    value[key.index('rain')] = filename

    
    #Save it
    data = np.column_stack([key, value])
    datafile_path = 'BG_param_K3.txt'
    np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=')
    

    #run BG_Flood for the appropriate Parameters file
    shutil.copy(datafile_path,'BG_param.txt')
    process=Popen('BG_flood_Kur.exe',stdout=PIPE, stderr=PIPE)
    stdout, stderr = process.communicate()
    
#For Buttinger
    #Modify the BG_param file
    value[key.index('outfile')] = 'Output_ButtingerC3.nc'

    
    #Save it
    data = np.column_stack([key, value])
    datafile_path = 'BG_param_B3.txt'
    np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=')
    

    #run BG_Flood for the appropriate Parameters file
    shutil.copy(datafile_path,'BG_param.txt')
    process=Popen('BG_flood_But.exe',stdout=PIPE, stderr=PIPE)
    stdout, stderr = process.communicate()
      

#For Kurgonov #Modify the BG_param file value[key.index('outfile')] = 'Output_KurgonovC3.nc' value[key.index('rain')] = filename #Save it data = np.column_stack([key, value]) datafile_path = 'BG_param_K3.txt' np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=') #run BG_Flood for the appropriate Parameters file shutil.copy(datafile_path,'BG_param.txt') process=Popen('BG_flood_Kur.exe',stdout=PIPE, stderr=PIPE) stdout, stderr = process.communicate() #For Buttinger #Modify the BG_param file value[key.index('outfile')] = 'Output_ButtingerC3.nc' #Save it data = np.column_stack([key, value]) datafile_path = 'BG_param_B3.txt' np.savetxt(datafile_path , data, fmt=['%s ',' %s ;'], delimiter='=') #run BG_Flood for the appropriate Parameters file shutil.copy(datafile_path,'BG_param.txt') process=Popen('BG_flood_But.exe',stdout=PIPE, stderr=PIPE) stdout, stderr = process.communicate()

  File "<ipython-input-19-84f3ecb6676b>", line 3
    value[key.index('outfile')] = 'Output_KurgonovC3.nc'
    ^
IndentationError: unexpected indent

Results¶

In [20]:

fn = 'Output_ButtingerC3.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_b3=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_b3.append(-ft)

fn = 'Output_ButtingerC3.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_b3=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [22]:

fn = 'Output_KurgonovC3.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]

F_k3=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_k3.append(-ft)

fn = 'Output_KurgonovC3.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] F_k3=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [23]:

ax=plt.subplot(1,1,1)
P1,=ax.plot(xC3,yC3,'o', label='Experiment')
P2,=ax.plot(tt, F_k3,'.-', label='zs using Kurganov')
P3,=ax.plot(tt, F_b3,'.-', label='zs using Buttinger')
ax2=ax.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax.set_ylabel('Flux ($m^3$/s)')
ax.set_xlim([0.0, 150.0])
ax.set_ylim([0.0, 0.0005])
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
plt.xlabel('Time(s)')
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');

ax=plt.subplot(1,1,1) P1,=ax.plot(xC3,yC3,'o', label='Experiment') P2,=ax.plot(tt, F_k3,'.-', label='zs using Kurganov') P3,=ax.plot(tt, F_b3,'.-', label='zs using Buttinger') ax2=ax.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax.set_ylabel('Flux ($m^3$/s)') ax.set_xlim([0.0, 150.0]) ax.set_ylim([0.0, 0.0005]) ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) plt.xlabel('Time(s)') handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3');

In [29]:

ax=plt.subplot(1,1,1)
P1,=ax.plot(xC3,yC3,'o', label='Experiment')
#P2,=ax.plot(tt, F_k3,'.-', label='zs using Kurganov')
P3,=ax.plot(tt, F_b3,'.-', label='BG_Flood')
ax2=ax.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax.set_ylabel('Flux ($m^3$/s)')
ax.set_xlim([0.0, 150.0])
ax.set_ylim([0.0, 0.0005])
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
plt.xlabel('Time(s)')
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
#handles.append(P2)
handles.append(P3)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');

ax=plt.subplot(1,1,1) P1,=ax.plot(xC3,yC3,'o', label='Experiment') #P2,=ax.plot(tt, F_k3,'.-', label='zs using Kurganov') P3,=ax.plot(tt, F_b3,'.-', label='BG_Flood') ax2=ax.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax.set_ylabel('Flux ($m^3$/s)') ax.set_xlim([0.0, 150.0]) ax.set_ylim([0.0, 0.0005]) ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) plt.xlabel('Time(s)') handles, labels = ax2.get_legend_handles_labels() handles.append(P1) #handles.append(P2) handles.append(P3) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3');

Comparing with the new friction implementation from Xia2018¶

(from [Xia2018] (https://www.sciencedirect.com/science/article/pii/S0309170818302124))

In [24]:

fn = 'Output_But-XiaC1.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_bx1=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_bx1.append(-ft)

fn = 'Output_But-XiaC1.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_bx1=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [25]:

fn = 'Output_But-XiaC3.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_bx3=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_bx3.append(-ft)

fn = 'Output_But-XiaC3.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_bx3=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [26]:

fn = 'Output_Kur-XiaC1.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kx1=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kx1.append(-ft)

fn = 'Output_Kur-XiaC1.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kx1=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [27]:

fn = 'Output_Kur-XiaC3.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kx3=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kx3.append(-ft)

fn = 'Output_Kur-XiaC3.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kx3=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [28]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P2,=ax1.plot(tt, F_kx1,'.-', label='zs using Kurganov and Xia')
P3,=ax1.plot(tt, F_bx1,'.-', label='zs using Buttinger and Xia')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
plt.legend(handles=handles)
plt.xlabel('Time(s)')
plt.title('Time evolution of the flux at the output of the domain, C1');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P2,=ax1.plot(tt, F_kx1,'.-', label='zs using Kurganov and Xia') P3,=ax1.plot(tt, F_bx1,'.-', label='zs using Buttinger and Xia') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) plt.legend(handles=handles) plt.xlabel('Time(s)') plt.title('Time evolution of the flux at the output of the domain, C1');

In [29]:

ax=plt.subplot(1,1,1)
P1,=ax.plot(xC3,yC3,'o', label='Experiment')
P2,=ax.plot(tt, F_kx3,'.-', label='zs using Kurganov and Xia')
P3,=ax.plot(tt, F_bx3,'.-', label='zs using Buttinger and Xia')
ax2=ax.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax.set_ylabel('Flux ($m^3$/s)')
ax.set_xlim([0.0, 150.0])
ax.set_ylim([0.0, 0.0005])
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
plt.xlabel('Time(s)')
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');

ax=plt.subplot(1,1,1) P1,=ax.plot(xC3,yC3,'o', label='Experiment') P2,=ax.plot(tt, F_kx3,'.-', label='zs using Kurganov and Xia') P3,=ax.plot(tt, F_bx3,'.-', label='zs using Buttinger and Xia') ax2=ax.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax.set_ylabel('Flux ($m^3$/s)') ax.set_xlim([0.0, 150.0]) ax.set_ylim([0.0, 0.0005]) ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) plt.xlabel('Time(s)') handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3');

In [30]:

ax=plt.subplot(1,1,1)
#P1,=ax.plot(xC3,yC3,'o', label='Experiment')
P1,=ax.plot(tt, (np.asanyarray(F_k)-np.asanyarray(F_kx1))/max(np.asanyarray(F_k))*100,'.-', label='Kurganov, C1, base - Xia (zs)')
P2,=ax.plot(tt, (np.asanyarray(F_b)-np.asanyarray(F_bx1))/max(np.asanyarray(F_k))*100,'.-', label='Buttinger, C1, base-Xia ')
P3,=ax.plot(tt, (np.asanyarray(F_k3)-np.asanyarray(F_kx3))/max(np.asanyarray(F_k))*100,'.-', label='Kurganov, C3, base - Xia (zs)')
P4,=ax.plot(tt, (np.asanyarray(F_b3)-np.asanyarray(F_bx3))/max(np.asanyarray(F_k))*100,'.-', label='Buttinger, C3, base-Xia ')
#ax2=ax.twinx()
#ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax.set_ylabel('% in Flux ($m^3$/s) difference between friction models')
#ax.set_xlim([0.0, 150.0])
#ax.set_ylim([0.0, 0.0005])
#ax2.set_ylabel('P (mm/h)')
#ax2.set_ylim([0.0, 500])
plt.xlabel('Time(s)')
#handles, labels = ax2.get_legend_handles_labels()
#handles.append(P1)
#handles.append(P2)
#handles.append(P3)
plt.legend(handles=[P1,P2,P3,P4])
plt.title('Comparison on the friction implementation : Base-Xia');

ax=plt.subplot(1,1,1) #P1,=ax.plot(xC3,yC3,'o', label='Experiment') P1,=ax.plot(tt, (np.asanyarray(F_k)-np.asanyarray(F_kx1))/max(np.asanyarray(F_k))*100,'.-', label='Kurganov, C1, base - Xia (zs)') P2,=ax.plot(tt, (np.asanyarray(F_b)-np.asanyarray(F_bx1))/max(np.asanyarray(F_k))*100,'.-', label='Buttinger, C1, base-Xia ') P3,=ax.plot(tt, (np.asanyarray(F_k3)-np.asanyarray(F_kx3))/max(np.asanyarray(F_k))*100,'.-', label='Kurganov, C3, base - Xia (zs)') P4,=ax.plot(tt, (np.asanyarray(F_b3)-np.asanyarray(F_bx3))/max(np.asanyarray(F_k))*100,'.-', label='Buttinger, C3, base-Xia ') #ax2=ax.twinx() #ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax.set_ylabel('% in Flux ($m^3$/s) difference between friction models') #ax.set_xlim([0.0, 150.0]) #ax.set_ylim([0.0, 0.0005]) #ax2.set_ylabel('P (mm/h)') #ax2.set_ylim([0.0, 500]) plt.xlabel('Time(s)') #handles, labels = ax2.get_legend_handles_labels() #handles.append(P1) #handles.append(P2) #handles.append(P3) plt.legend(handles=[P1,P2,P3,P4]) plt.title('Comparison on the friction implementation : Base-Xia');

In [31]:

ax=plt.subplot(1,1,1)
#P1,=ax.plot(xC3,yC3,'o', label='Experiment')
P1,=ax.plot(tt, (np.asanyarray(F_k)-np.asanyarray(F_b))/max(np.asanyarray(F_k))*100,'.-', label='Base friction, C1, Kur.-But. (zs)')
P2,=ax.plot(tt, (np.asanyarray(F_k3)-np.asanyarray(F_b3))/max(np.asanyarray(F_k3))*100,'.-', label='Base friction, C3, Kur.-But. ')
P3,=ax.plot(tt, (np.asanyarray(F_kx1)-np.asanyarray(F_bx1))/max(np.asanyarray(F_kx1))*100,'.-', label='Xia, C1, Kur.-But.')
P4,=ax.plot(tt, (np.asanyarray(F_kx3)-np.asanyarray(F_bx3))/max(np.asanyarray(F_kx3))*100,'.-', label='Xia, C3, Kur.-But. ')
#ax2=ax.twinx()
#ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax.set_ylabel('% in Flux ($m^3$/s) difference between surface reconstruction models')
#ax.set_xlim([0.0, 150.0])
#ax.set_ylim([0.0, 0.0005])
#ax2.set_ylabel('P (mm/h)')
#ax2.set_ylim([0.0, 500])
plt.xlabel('Time(s)')
#handles, labels = ax2.get_legend_handles_labels()
#handles.append(P1)
#handles.append(P2)
#handles.append(P3)
plt.legend(handles=[P1,P2,P3,P4])
plt.title('Comparison on the surface reconstruction: Kurganov-Buttinger');

ax=plt.subplot(1,1,1) #P1,=ax.plot(xC3,yC3,'o', label='Experiment') P1,=ax.plot(tt, (np.asanyarray(F_k)-np.asanyarray(F_b))/max(np.asanyarray(F_k))*100,'.-', label='Base friction, C1, Kur.-But. (zs)') P2,=ax.plot(tt, (np.asanyarray(F_k3)-np.asanyarray(F_b3))/max(np.asanyarray(F_k3))*100,'.-', label='Base friction, C3, Kur.-But. ') P3,=ax.plot(tt, (np.asanyarray(F_kx1)-np.asanyarray(F_bx1))/max(np.asanyarray(F_kx1))*100,'.-', label='Xia, C1, Kur.-But.') P4,=ax.plot(tt, (np.asanyarray(F_kx3)-np.asanyarray(F_bx3))/max(np.asanyarray(F_kx3))*100,'.-', label='Xia, C3, Kur.-But. ') #ax2=ax.twinx() #ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax.set_ylabel('% in Flux ($m^3$/s) difference between surface reconstruction models') #ax.set_xlim([0.0, 150.0]) #ax.set_ylim([0.0, 0.0005]) #ax2.set_ylabel('P (mm/h)') #ax2.set_ylim([0.0, 500]) plt.xlabel('Time(s)') #handles, labels = ax2.get_legend_handles_labels() #handles.append(P1) #handles.append(P2) #handles.append(P3) plt.legend(handles=[P1,P2,P3,P4]) plt.title('Comparison on the surface reconstruction: Kurganov-Buttinger');

Testing the Smart Friction model¶

The friction model has been changed from Manning (used precendently as in the CEA paper) to the Smart model, used here with a roughness height h=1E-6m (midrange estimation for stainless steel).

In [32]:

fn = 'Output_KurC1-Smart1-6.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS1=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS1.append(-ft)

fn = 'Output_KurC1-Smart1-6.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS1=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [36]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P2,=ax1.plot(tt, F_k,'.-', label='zs using Kurganov')
P3,=ax1.plot(tt, F_b,'.-', label='zs using Buttinger')
P4,=ax1.plot(tt, F_kx1,'.-', label='zs using Kurganov and Xia')
P5,=ax1.plot(tt, F_bx1,'.-', label='zs using Buttinger and Xia')
P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. basic frict. with Smart')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
handles.append(P4)
handles.append(P5)
handles.append(P6)
plt.legend(handles=handles)
plt.xlabel('Time(s)')
plt.title('Time evolution of the flux at the output of the domain, C1');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P2,=ax1.plot(tt, F_k,'.-', label='zs using Kurganov') P3,=ax1.plot(tt, F_b,'.-', label='zs using Buttinger') P4,=ax1.plot(tt, F_kx1,'.-', label='zs using Kurganov and Xia') P5,=ax1.plot(tt, F_bx1,'.-', label='zs using Buttinger and Xia') P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. basic frict. with Smart') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) handles.append(P4) handles.append(P5) handles.append(P6) plt.legend(handles=handles) plt.xlabel('Time(s)') plt.title('Time evolution of the flux at the output of the domain, C1');

In [43]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. basic frict. with Smart')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P6)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C1');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. basic frict. with Smart') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P6) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C1');

In [34]:

fn = 'Output_KurC3-Smart1-6.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS3=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS3.append(-ft)

fn = 'Output_KurC3-Smart1-6.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS3=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [38]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC3,yC3,'o', label='Experiment')
P2,=ax1.plot(tt, F_k3,'.-', label='zs using Kurganov')
P3,=ax1.plot(tt, F_b3,'.-', label='zs using Buttinger')
P4,=ax1.plot(tt, F_kx3,'.-', label='zs using Kurganov and Xia')
P5,=ax1.plot(tt, F_bx3,'.-', label='zs using Buttinger and Xia')
P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. basic fric. Smart')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P2)
handles.append(P3)
handles.append(P4)
handles.append(P5)
handles.append(P6)
plt.legend(handles=handles)
plt.xlabel('Time(s)')
plt.title('Time evolution of the flux at the output of the domain, C3');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC3,yC3,'o', label='Experiment') P2,=ax1.plot(tt, F_k3,'.-', label='zs using Kurganov') P3,=ax1.plot(tt, F_b3,'.-', label='zs using Buttinger') P4,=ax1.plot(tt, F_kx3,'.-', label='zs using Kurganov and Xia') P5,=ax1.plot(tt, F_bx3,'.-', label='zs using Buttinger and Xia') P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. basic fric. Smart') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P2) handles.append(P3) handles.append(P4) handles.append(P5) handles.append(P6) plt.legend(handles=handles) plt.xlabel('Time(s)') plt.title('Time evolution of the flux at the output of the domain, C3');

In [45]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC3,yC3,'o', label='Experiment')
P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. basic fric. Smart')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P6)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC3,yC3,'o', label='Experiment') P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. basic fric. Smart') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P6) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3');

Test by increasing the roughtness height¶

In [46]:

fn = 'Output_KurC1-Smart5-6.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS1b=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS1b.append(-ft)

fn = 'Output_KurC1-Smart5-6.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS1b=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [53]:

fn = 'Output_KurC1-Smart5-6_eps-5.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS1c=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS1c.append(-ft)

fn = 'Output_KurC1-Smart5-6_eps-5.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS1c=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [55]:

fn = 'Output_KurC1-Smart5-6_theta-14.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS1d=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS1d.append(-ft)

fn = 'Output_KurC1-Smart5-6_theta-14.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS1d=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [62]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. Smart, 1E-6')
P7,=ax1.plot(tt, F_kS1b,'.-', label='zs using Kur. Smart,5E-6')
P8,=ax1.plot(tt, F_kS1c,'.-', label='zs using Kur. Smart,5E-6, eps=1E-5')
P9,=ax1.plot(tt, F_kS1d,'.-', label='zs using Kur., theta=1.4')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P6)
handles.append(P7)
handles.append(P8)
handles.append(P9)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C1');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P6,=ax1.plot(tt, F_kS1,'.-', label='zs using Kur. Smart, 1E-6') P7,=ax1.plot(tt, F_kS1b,'.-', label='zs using Kur. Smart,5E-6') P8,=ax1.plot(tt, F_kS1c,'.-', label='zs using Kur. Smart,5E-6, eps=1E-5') P9,=ax1.plot(tt, F_kS1d,'.-', label='zs using Kur., theta=1.4') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P6) handles.append(P7) handles.append(P8) handles.append(P9) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C1');

In [48]:

fn = 'Output_KurC3-Smart5-6.nc'
ds = nc4.Dataset(fn)
ds
xx=ds['xx_P0'][:]
yy=ds['yy_P0'][:]
tt=ds['time'][:]
hh=ds['h_P0'][:,:,:]
zsf=ds['zs_P0'][:,:,:]
zb=ds['zb_P0'][:,:,:]
vv=ds['v_P0'][:]
dy=yy[1]-yy[0]
dx=xx[1]-xx[0]

F_kS3b=[]
for k in range(len(tt)):
    ft=0;
    for i in range(len(xx)):
        for j in range(len(yy)):
            if xx[i]>0.5 and xx[i]<1.5:
                if abs(yy[j]-0.1)<dy/2:
                    ft=ft+hh[k,j,i]*vv[k,j,i]*dx
    F_kS3b.append(-ft)

fn = 'Output_KurC3-Smart5-6.nc' ds = nc4.Dataset(fn) ds xx=ds['xx_P0'][:] yy=ds['yy_P0'][:] tt=ds['time'][:] hh=ds['h_P0'][:,:,:] zsf=ds['zs_P0'][:,:,:] zb=ds['zb_P0'][:,:,:] vv=ds['v_P0'][:] dy=yy[1]-yy[0] dx=xx[1]-xx[0] F_kS3b=[] for k in range(len(tt)): ft=0; for i in range(len(xx)): for j in range(len(yy)): if xx[i]>0.5 and xx[i]<1.5: if abs(yy[j]-0.1)

In [70]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC3,yC3,'o', label='Experiment')
#P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. Smart, 1E-6')
P7,=ax1.plot(tt, F_kS3b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
#handles.append(P6)
handles.append(P7)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC3,yC3,'o', label='Experiment') #P6,=ax1.plot(tt, F_kS3,'.-', label='zs using Kur. Smart, 1E-6') P7,=ax1.plot(tt, F_kS3b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) #handles.append(P6) handles.append(P7) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3');

Figures for the paper "Coast and Port"¶

In [77]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC1,yC1,'o', label='Experiment')
P7,=ax1.plot(tt, F_kS1b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P7)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C1');
plt.savefig('CEA_C1_Kur-Smart.png',bbox_inches='tight')

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC1,yC1,'o', label='Experiment') P7,=ax1.plot(tt, F_kS1b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr, Hr,alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P7) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C1'); plt.savefig('CEA_C1_Kur-Smart.png',bbox_inches='tight')

In [78]:

ax1=plt.subplot(1,1,1)
P1,=ax1.plot(xC3,yC3,'o', label='Experiment')
P7,=ax1.plot(tt, F_kS3b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m')
ax1.set_ylabel('Flux ($m^3$/s)')
ax1.set_xlabel('Time(s)')
ax1.set_xlim([0.0, 150.0])
ax1.set_ylim([0.0, 0.0005])
ax2=ax1.twinx()
ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing')
ax2.set_ylabel('P (mm/h)')
ax2.set_ylim([0.0, 500])
#ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500))
handles, labels = ax2.get_legend_handles_labels()
handles.append(P1)
handles.append(P7)
plt.legend(handles=handles)
plt.title('Time evolution of the flux at the output of the domain, C3');
plt.savefig('CEA_C3_Kur-Smart.png',bbox_inches='tight')

ax1=plt.subplot(1,1,1) P1,=ax1.plot(xC3,yC3,'o', label='Experiment') P7,=ax1.plot(tt, F_kS3b,'.-', label='BG_flood, $ks=5\, 10^{-6}$m') ax1.set_ylabel('Flux ($m^3$/s)') ax1.set_xlabel('Time(s)') ax1.set_xlim([0.0, 150.0]) ax1.set_ylim([0.0, 0.0005]) ax2=ax1.twinx() ax2.fill_between(Tr3, Hr3, [0]*len(Tr3), alpha=0.1, color='blue', label='Rain forcing') ax2.set_ylabel('P (mm/h)') ax2.set_ylim([0.0, 500]) #ax2 = plt.axes(xlim=(0.0, 150.0), ylim=(0.0, 500)) handles, labels = ax2.get_legend_handles_labels() handles.append(P1) handles.append(P7) plt.legend(handles=handles) plt.title('Time evolution of the flux at the output of the domain, C3'); plt.savefig('CEA_C3_Kur-Smart.png',bbox_inches='tight')

In [ ]: