Major update V3

.gitignore

@@ -105,4 +105,7 @@ cython_debug/
 .vscode/
 
 # Developer Notes
-dev-notes/
+dev-notes/
+*.nc
+*.mat
+*.log

aeolis/avalanching.py

@@ -184,6 +184,10 @@ def avalanche(s, p):
 
                 s['zb'] += E * (q_in - q_out)
 
+        # Ensure water level is up-to-date with bed level
+        s['zs'] = s['SWL'].copy()
+        ix = (s['zb'] > s['zs'])
+        s['zs'][ix] = s['zb'][ix]
     return s	    
 
 

aeolis/bed.py

@@ -83,7 +83,12 @@ def initialize(s, p):
     if isinstance(p['grain_dist'], str):
             logger.log_and_raise('Grain size file not recognized as array, check file path and whether all values have been filled in.', exc=ValueError) 
 
-    if p['bedcomp_file'] is None and p['grain_dist'].ndim == 1 and p['grain_dist'].dtype == 'float64' or p['grain_dist'].dtype == 'int': 
+    if p['bedcomp_file'] is not None and p['supply_file'] is not None :
+            logger.log_and_raise('Conflict in input definition, cannot define supply_file and bedcomp_file simultaneously', exc=ValueError) 
+
+    if p['supply_file'] is not None:
+        s['mass'][:,:,:,:] = 0 #p['supply_file'].reshape(s['mass'].shape)                
+    elif p['bedcomp_file'] is None and p['grain_dist'].ndim == 1 and p['grain_dist'].dtype == 'float64' or p['grain_dist'].dtype == 'int': 
         # Both float and int are included as options for the grain dist to make sure there is no error when grain_dist is filled in as 1 instead of 1.0. 
         for i in range(nl):
             gs = makeiterable(p['grain_dist'])
@@ -213,7 +218,7 @@ def wet_bed_reset(s, p):
 
         Tbedreset = p['dt_opt'] / p['Tbedreset']
 
-        ix = s['zs'] > (s['zb'] + 0.01)
+        ix = s['TWL'] > (s['zb'])
         s['zb'][ix] += (s['zb0'][ix] - s['zb'][ix]) * Tbedreset
 
     return s
@@ -248,6 +253,16 @@ def update(s, p):
         Spatial grids
 
     '''
+    # this is where a supply file is used, this in only for simple cases.
+    if type(p['supply_file']) == np.ndarray:
+        # in descrete supply limited conditions the bed bed layer operations are not valid. 
+        s['mass'][:,:,0,0] -= s['pickup'][:,:,0]
+        s['mass'][:,:,0,0] += p['supply_file']*p['dt_opt']
+        # reset supply under water if process tide is active
+        if p['process_tide']:
+            s['mass'][(s['zb']< s['zs']),0,0]=0
+        return s
+
 
     nx = p['nx']
     ny = p['ny']
@@ -265,7 +280,8 @@ def update(s, p):
     ix_dep = dm[:,0] > 0.
 
     # reshape mass matrix
-    m = s['mass'].reshape((-1,nl,nf))
+    m = s['mass'].reshape((-1,nl,nf)).copy()
+
 
     # negative mass may occur in case of deposition due to numerics,
     # which should be prevented
@@ -462,8 +478,11 @@ def average_change(l, s, p):
     s['dzbavg'] = n*s['dzbyear']+(1-n)*l['dzbavg']
 
     # Calculate average bed level change as input for vegetation growth [m/year]
-    # s['dzbveg'] = s['dzbavg'].copy()
-    s['dzbveg'] = s['dzbyear'].copy()
+    s['dzbveg'] = s['dzbavg'].copy()
+    # s['dzbveg'] = s['dzbyear'].copy()
+
+    if p['_time'] < p['avg_time']:
+        s['dzbveg'] *= 0.
 
 
     return s

aeolis/constants.py

@@ -106,8 +106,9 @@
         'hveg',                             # [m] height of vegetation
         'dhveg',                            # [m] Difference in vegetation height per time step
         'dzbveg',                           # [m] Bed level change used for calculation of vegetation growth
-        'germinate',                        # vegetation germination
-        'lateral',                          # vegetation lateral expansion
+        'germinate',                        # [bool] Newly vegetated due to germination (or establishment) 
+        'lateral',                          # [bool] Newly vegetated due to lateral propagation 
+        'vegetated',                        # [bool] Vegetated, determines if vegetation growth or burial is allowed
         'vegfac',                           # Vegetation factor to modify shear stress by according to Raupach 1993
         'fence_height',                     # Fence height
         'R',                                # [m] wave runup
@@ -174,15 +175,15 @@
     'process_meteo'                 : False,              # Enable the process of meteo
     'process_salt'                  : False,              # Enable the process of salt
     'process_humidity'              : False,              # Enable the process of humidity
-    'process_groundwater'           : False,        #NEWCH      # Enable the process of groundwater
-    'process_scanning'              : False,         #NEWCH      # Enable the process of scanning curves
+    'process_groundwater'           : False,              #NEWCH      # Enable the process of groundwater
+    'process_scanning'              : False,              #NEWCH      # Enable the process of scanning curves
     'process_inertia'               : False,              # NEW
     'process_separation'            : False,               # Enable the including of separation bubble
     'process_vegetation'            : False,              # Enable the process of vegetation
     'process_fences'                : False,              # Enable the process of sand fencing
     'process_dune_erosion'          : False,              # Enable the process of wave-driven dune erosion
     'process_seepage_face'          : False,              # Enable the process of groundwater seepage (NB. only applicable to positive beach slopes)
-    'visualization'                 : False,               # Boolean for visualization of model interpretation before and just after initialization
+    'visualization'                 : False,              # Boolean for visualization of model interpretation before and just after initialization
     'xgrid_file'                    : None,               # Filename of ASCII file with x-coordinates of grid cells
     'ygrid_file'                    : None,               # Filename of ASCII file with y-coordinates of grid cells
     'bed_file'                      : None,               # Filename of ASCII file with bed level heights of grid cells
@@ -195,11 +196,13 @@
     'fence_file'                    : None,               # Filename of ASCII file with sand fence location/height (above the bed)
     'ne_file'                       : None,               # Filename of ASCII file with non-erodible layer
     'veg_file'                      : None,               # Filename of ASCII file with initial vegetation density
+    'supply_file'                   : None,               # Filename of ASCII file with a manual definition of sediment supply (mainly used in academic cases)
     'wave_mask'                     : None,               # Filename of ASCII file with mask for wave height
     'tide_mask'                     : None,               # Filename of ASCII file with mask for tidal elevation
     'runup_mask'                    : None,               # Filename of ASCII file with mask for run-up
     'threshold_mask'                : None,               # Filename of ASCII file with mask for the shear velocity threshold
     'gw_mask'                       : None,         #NEWCH      # Filename of ASCII file with mask for the groundwater level
+    'vver_mask'                     : None,         #NEWBvW      # Filename of ASCII file with mask for the vertical vegetation growth    
     'nx'                            : 0,                  # [-] Number of grid cells in x-dimension
     'ny'                            : 0,                  # [-] Number of grid cells in y-dimension
     'dt'                            : 60.,                # [s] Time step size

aeolis/inout.py

@@ -426,19 +426,18 @@ def visualize_timeseries(p, t):
     axs[1].set_title('Wind direction, udir (deg)')
 
     # Read the user input (waves)
-    if p['wave_file'] is not None:
+
+    if np.shape(p['wave_file'])[1]== 3:
         w_t = p['wave_file'][:,0]
         w_Hs = p['wave_file'][:,1]
+        w_Tp = p['wave_file'][:,2]
         axs[2].plot(w_t, w_Hs, 'k')
+        axs[3].plot(w_t, w_Tp, 'k')
         axs[2].set_title('Wave height, Hs (m)')
-        if np.shape(p['wave_file'])[1] == 3:
-            w_Tp = p['wave_file'][:,2]
-            axs[3].plot(w_t, w_Tp, 'k')
-            axs[3].set_title('Wave period, Tp (sec)')
-
+        axs[3].set_title('Wave period, Tp (sec)')
 
     # Read the user input (tide)
-    if p['tide_file'] is not None:
+    if np.shape(p['tide_file'])[1]==2:
         T_t = p['tide_file'][:,0]
         T_zs = p['tide_file'][:,1]
         axs[4].plot(T_t, T_zs, 'k')
@@ -501,7 +500,8 @@ def visualize_spatial(s, p):
         pcs[0][2] = axs[0,2].pcolormesh(x, y, s['rhoveg'], cmap='Greens', clim= [0, 1])
         pcs[1][0] = axs[1,0].pcolormesh(x, y, s['uw'], cmap='plasma')
         pcs[1][1] = axs[1,1].pcolormesh(x, y, s['ustar'], cmap='plasma')
-        pcs[1][2] = axs[1,2].pcolormesh(x, y, s['tau'], cmap='plasma')
+        # pcs[1][2] = axs[1,2].pcolormesh(x, y, s['tau'], cmap='plasma')
+        pcs[1][2] = axs[1,2].pcolormesh(x, y, s['u'][:, :, 0], cmap='plasma')
         pcs[2][0] = axs[2,0].pcolormesh(x, y, s['moist'], cmap='Blues', clim= [0, 0.4])
         pcs[2][1] = axs[2,1].pcolormesh(x, y, s['gw'], cmap='viridis')
         pcs[2][2] = axs[2,2].pcolormesh(x, y, s['uth'][:,:,0], cmap='plasma')
@@ -532,15 +532,17 @@ def visualize_spatial(s, p):
     skip = 10
     axs[1,0].quiver(x[::skip, ::skip], y[::skip, ::skip], s['uws'][::skip, ::skip], s['uwn'][::skip, ::skip])
     axs[1,1].quiver(x[::skip, ::skip], y[::skip, ::skip], s['ustars'][::skip, ::skip], s['ustarn'][::skip, ::skip])
-    axs[1,2].quiver(x[::skip, ::skip], y[::skip, ::skip], s['taus'][::skip, ::skip], s['taun'][::skip, ::skip])
+    # axs[1,2].quiver(x[::skip, ::skip], y[::skip, ::skip], s['taus'][::skip, ::skip], s['taun'][::skip, ::skip])
+    axs[1,2].quiver(x[::skip, ::skip], y[::skip, ::skip], s['us'][::skip, ::skip, 0], s['un'][::skip, ::skip, 0])
 
     # Adding titles to the plots
     axs[0,0].set_title('Bed level, zb (m)')
     axs[0,1].set_title('Non-erodible layer, zne (m)')
     axs[0,2].set_title('Vegetation density, rhoveg (-)')
     axs[1,0].set_title('Wind velocity, uw (m/s)')
     axs[1,1].set_title('Shear velocity, ustar (m/s)')
-    axs[1,2].set_title('Shear stress, tau (N/m2)')
+    # axs[1,2].set_title('Shear stress, tau (N/m2)')
+    axs[1,2].set_title('Grain velocity, u (m/s)')
     axs[2,0].set_title('Soil moisture content, (-)')
     axs[2,1].set_title('Ground water level, gw (m)')
     axs[2,2].set_title('Velocity threshold (0th fraction), uth (m/s)')