Merge branch 'main' into aeolis-dynarevdev

Author: christavanijzendoorn

.devcontainer/Dockerfile

@@ -0,0 +1,21 @@
+FROM continuumio/miniconda3
+
+# Set the working directory
+
+# WORKDIR /workspace
+
+# Install Python
+RUN conda install python -y
+
+# Copy the project files
+# COPY . .
+
+# Install the necessary packages
+
+# RUN pip install .
+
+# Clean up
+RUN conda clean -afy
+
+# Set the default command
+CMD ["bash"]

.devcontainer/devcontainer.json

@@ -0,0 +1,7 @@
+{
+  "name": "AEOLIS Development Container",
+  "build": {
+    "dockerfile": "Dockerfile",
+    "context": ".."
+  }
+}

.gitignore

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

CITATION.cff

@@ -69,5 +69,5 @@ keywords:
   - sediment
 license: GPL-3.0
 commit: 0adaedfc8015f18b3b551bed1dda38d630cd8c95
-version: 3.0.0.rc1
-date-released: '2023-10-18'
+version: 3.0.1
+date-released: '2025-07-11'

README.md

@@ -91,7 +91,7 @@ Read our [Contribution Guidelines](CONTRIBUTING.md) to know how you can help to
 
 Please, cite this software as follows:
 
-*de Vries, S., Hallin, C., van IJzendoorn, C., van Westen, B., Cohn, N., Strypsteen, G., Skaden, J., Agrawal, N., & Garcia Alvarez, M. (2023). AeoLiS (Version 3.0.0.rc1) [Computer software]. https://github.com/openearth/aeolis-python*
+*de Vries, S., Hallin, C., van IJzendoorn, C., van Westen, B., Cohn, N., Strypsteen, G., Skaden, J., Agrawal, N., & Garcia Alvarez, M. (2023). AeoLiS (Version 3.0.1) [Computer software]. https://github.com/openearth/aeolis-python*
 
 ## Acknowlegdements
 

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'])
@@ -100,8 +105,14 @@ def initialize(s, p):
             for j in range(nf):
                 s['mass'][:,:,i,j] = p['rhog'] * (1. - p['porosity']) \
                                      * s['thlyr'][:,:,i] * gs[j]
-    else:
-        s['mass'][:,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape)                
+    else:      
+        # if a quasi 2D domain is used, the bedcomp_file is reshaped to fit the domain
+        if s['mass'].shape[0] == 3:
+            s['mass'][0,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape[1:])
+            s['mass'][1,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape[1:])
+            s['mass'][2,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape[1:])                                                      
+        else:
+            s['mass'][:,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape)                
 
     # initialize masks
     for k, v in p.items():
@@ -200,7 +211,7 @@ def mixtoplayer(s, p):
 
 
 def wet_bed_reset(s, p):
-    ''' Text
+    ''' Reset wet bed to initial bed level if the total water level is above the bed level.
 
 
 
@@ -222,7 +233,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
@@ -257,6 +268,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']
@@ -274,7 +295,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
@@ -298,8 +320,11 @@ def update(s, p):
     if p['grain_dist'].ndim == 2: 
         m[ix_ero,-1,:] -= dm[ix_ero,:] * normalize(p['grain_dist'][-1,:])[np.newaxis,:].repeat(np.sum(ix_ero), axis=0)
     elif type(p['bedcomp_file']) == np.ndarray:
-        gs = p['bedcomp_file'].reshape((-1,nl,nf))
-        m[ix_ero,-1,:] -= dm[ix_ero,:] * normalize(gs[ix_ero,-1, :], axis=1)
+        gs = np.zeros(s['mass'].shape)
+        gs[0,:,:,:] = p['bedcomp_file'].reshape((-1,nl,nf))
+        gs[1,:,:,:] = p['bedcomp_file'].reshape((-1,nl,nf))
+        gs[2,:,:,:] = p['bedcomp_file'].reshape((-1,nl,nf))
+        m[ix_ero,-1,:] -= dm[ix_ero,:] * normalize(gs.reshape((-1,nl,nf))[ix_ero,-1, :], axis=1)
     else:
         m[ix_ero,-1,:] -= dm[ix_ero,:] * normalize(p['grain_dist'])[np.newaxis,:].repeat(np.sum(ix_ero), axis=0)
     # remove tiny negatives
@@ -316,6 +341,9 @@ def update(s, p):
     # reshape mass matrix
     s['mass'] = m.reshape((ny+1,nx+1,nl,nf))
 
+    # Store toplayer of 'mass' variable (ilayer = 0)
+    s['masstop'][:,:,:] = s['mass'][:,:,0,:].copy()
+
     # update bathy
     if p['process_bedupdate']:
 
@@ -475,8 +503,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