1d2d patch (#207)

fixes issue #190 

* Update model.py

* initiate quasi 2D implementation

* quasi 2D with sorting

* handy caching for numba speeds up compilation

* write quasi 2D output in 1D nc-file

* New rotation scheme

Capabilities added to port 1D parameterizations into 2D using a central rotation grid

---------

Co-authored-by: Nick <[email protected]>

Author: Sierd

aeolis/bed.py

@@ -105,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():
@@ -305,8 +311,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

aeolis/constants.py

@@ -322,6 +322,7 @@
     'method_transport'              : 'bagnold',          # Name of method to compute equilibrium sediment transport rate
     'method_roughness'              : 'constant',         # Name of method to compute the roughness height z0, note that here the z0 = k, which does not follow the definition of Nikuradse where z0 = k/30.
     'method_grainspeed'             : 'windspeed',        # Name of method to assume/compute grainspeed (windspeed, duran, constant)
+    'method_shear'                  : 'fft',              # Name of method to compute topographic effects on wind shear stress (fft, quasi2d, duna2d (experimental))
     'max_error'                     : 1e-6,               # [-] Maximum error at which to quit iterative solution in implicit numerical schemes
     'max_iter'                      : 1000,               # [-] Maximum number of iterations at which to quit iterative solution in implicit numerical schemes
     'max_iter_ava'                  : 1000,               # [-] Maximum number of iterations at which to quit iterative solution in avalanching calculation

aeolis/examples/grainsizevariations/aeolis_horizontalgradient1.txt

@@ -87,7 +87,7 @@ Cdk                     = 5.                 					% [-] Constant in DK formulati
 
 %% -------------------- [Solver] ----------------------------- %%
 T 			= 1.							% [s] Adaptation time scale in advection equation
-solver 			= trunk							% [-] Numerical solver of advection scheme
+solver 			= trunk   						% [-] Numerical solver of advection scheme
 CFL                     = 1.                 					% [-] CFL number to determine time step in explicit scheme
 accfac 			= 1.							% [-] Numerical acceleration factor
 scheme 			= euler_backward 					% [-] Name of numerical scheme (euler_forward, euler_backward or crank_nicolson)

aeolis/examples/grainsizevariations/aeolis_horizontalgradient2.txt

@@ -87,7 +87,7 @@ Cdk                     = 5.                 					% [-] Constant in DK formulati
 
 %% -------------------- [Solver] ----------------------------- %%
 T 			= 1.							% [s] Adaptation time scale in advection equation
-solver 			= trunk							% [-] Numerical solver of advection scheme
+solver 			= trunk                                                 % [-] Numerical solver of advection scheme
 CFL                     = 1.                 					% [-] CFL number to determine time step in explicit scheme
 accfac 			= 1.							% [-] Numerical acceleration factor
 scheme 			= euler_backward 					% [-] Name of numerical scheme (euler_forward, euler_backward or crank_nicolson)

aeolis/examples/grainsizevariations/aeolis_verticallayering1.txt

@@ -86,7 +86,7 @@ Cdk                     = 5.                 					% [-] Constant in DK formulati
 
 %% -------------------- [Solver] ----------------------------- %%
 T 			= 1.							% [s] Adaptation time scale in advection equation
-solver 			= trunk							% [-] Numerical solver of advection scheme
+solver 			= trunk 						% [-] Numerical solver of advection scheme
 CFL                     = 1.                 					% [-] CFL number to determine time step in explicit scheme
 accfac 			= 1.							% [-] Numerical acceleration factor
 scheme 			= euler_backward 					% [-] Name of numerical scheme (euler_forward, euler_backward or crank_nicolson)

aeolis/examples/grainsizevariations/aeolis_verticallayering2.txt

@@ -86,7 +86,7 @@ Cdk                     = 5.                 					% [-] Constant in DK formulati
 
 %% -------------------- [Solver] ----------------------------- %%
 T 			= 1.							% [s] Adaptation time scale in advection equation
-solver 			= trunk							% [-] Numerical solver of advection scheme
+solver 			= trunk  						% [-] Numerical solver of advection scheme
 CFL                     = 1.                 					% [-] CFL number to determine time step in explicit scheme
 accfac 			= 1.							% [-] Numerical acceleration factor
 scheme 			= euler_backward 					% [-] Name of numerical scheme (euler_forward, euler_backward or crank_nicolson)

aeolis/inout.py

@@ -117,7 +117,7 @@ def read_configfile(configfile, parse_files=True, load_defaults=True):
     # set default for nsavetimes, if not given
     if 'nsavetimes' in p and not p['nsavetimes']:
         p['nsavetimes'] = int(p['dzb_interval']/p['dt'])
-
+    
     return p
 
 

aeolis/model.py

@@ -189,7 +189,7 @@ def initialize(self)-> None:
         self.p = aeolis.inout.read_configfile(self.configfile)
         aeolis.inout.check_configuration(self.p)
 
-        # set nx, ny and nfractions
+        # set nx, ny and nfractions and make 1D into quasi 2D
         if self.p['xgrid_file'].ndim == 2:
             self.p['ny'], self.p['nx'] = self.p['xgrid_file'].shape
 
@@ -198,9 +198,29 @@ def initialize(self)-> None:
             self.p['ny'] -= 1
 
         else:
-            self.p['nx'] = len(self.p['xgrid_file'])
-            self.p['nx'] -= 1 
-            self.p['ny'] = 0
+            if 0:
+                #this is the old 1D stuff
+                self.p['nx'] = len(self.p['xgrid_file'])
+                self.p['nx'] -= 1 
+                self.p['ny'] = 0
+
+            if 1:
+                # this is the new quasi 2D stuff
+                # this is where we make the 1D grid into a 2D grid to ensure process compatibility
+                self.p['xgrid_file'] = np.transpose(np.stack((self.p['xgrid_file'], self.p['xgrid_file'], self.p['xgrid_file']),axis=1))
+                self.p['ny'], self.p['nx'] = self.p['xgrid_file'].shape
+
+                dy = self.p['xgrid_file'][1,2]-self.p['xgrid_file'][1,1]
+
+                # repeat the above for ygrid_file
+                self.p['ygrid_file'] = np.transpose(np.stack((self.p['ygrid_file'], self.p['ygrid_file']+dy, self.p['ygrid_file']+2*dy),axis=1))
+                
+                # repeat the above for bed_file
+                self.p['bed_file'] = np.transpose(np.stack((self.p['bed_file'], self.p['bed_file'], self.p['bed_file']),axis=1))
+
+                # change from number of points to number of cells
+                self.p['nx'] -= 1  
+                self.p['ny'] -= 1
 
         #self.p['nfractions'] = len(self.p['grain_dist'])
         self.p['nfractions'] = len(self.p['grain_size'])
@@ -266,10 +286,9 @@ def update(self, dt:float=-1) -> None:
 
         '''
 
-        self.p['_time'] = self.t    
+        self.p['_time'] = self.t 
+  
 
-        # here are going to make a change
-        #     
 
         # store previous state
         self.l = self.s.copy()
@@ -1631,7 +1650,12 @@ def solve_SS(self, alpha:float=0., beta:float=1.) -> dict:
                 w = w_init.copy()
             else:
                 # use initial guess for first time step
-                w = p['grain_dist'].reshape((1,1,-1))
+                # when p['grain_dist'] has 2 dimensions take the first row otherwise take the only row
+                if len(p['grain_dist'].shape) == 2:
+                    w = p['grain_dist'][0,:].reshape((1,1,-1))
+                else:
+                    w = p['grain_dist'].reshape((1,1,-1))
+                    
                 w = w.repeat(p['ny']+1, axis=0)
                 w = w.repeat(p['nx']+1, axis=1)
         else:
@@ -1678,7 +1702,7 @@ def solve_SS(self, alpha:float=0., beta:float=1.) -> dict:
                     Ct[-1,:,0] =  -2
 
                 # Ct, pickup = sweep(s['Cu'].copy(), s['mass'].copy(), self.dt, p['T'], s['ds'], s['dn'], s['us'], s['un'] )
-                Ct, pickup = sweep3(Ct, s['Cu'].copy(), s['mass'].copy(), self.dt, p['T'], s['ds'], s['dn'], s['us'], s['un'] )
+                Ct, pickup = sweep3(Ct, s['Cu'].copy(), s['mass'].copy(), self.dt, p['T'], s['ds'], s['dn'], s['us'], s['un'],w)
 
             if 0:
                 #define 4 quadrants based on wind directions

aeolis/netcdf.py

@@ -84,8 +84,14 @@ def initialize(outputfile, outputvars, s, p, dimensions):
     with netCDF4.Dataset(outputfile, 'w') as nc:
 
         # add dimensions
-        nc.createDimension('s', p['nx']+1)
-        nc.createDimension('n', p['ny']+1)
+        nc.createDimension('s', p['nx']+1)  
+
+        # this is where a quasi 2D model is stored in 1D    
+        if p['ny']+1 == 3:
+            nc.createDimension('n', 1)
+        else:
+            nc.createDimension('n', p['ny']+1)
+        
         nc.createDimension('time', 0)
         nc.createDimension('nv', 2)
         nc.createDimension('nv2', 4)
@@ -287,9 +293,17 @@ def initialize(outputfile, outputvars, s, p, dimensions):
 
         # store static data
         nc.variables['s'][:] = np.arange(p['nx']+1)
-        nc.variables['n'][:] = np.arange(p['ny']+1)
-        nc.variables['x'][:,:] = s['x']
-        nc.variables['y'][:,:] = s['y']
+
+        # this is where a quasi 2D model is stored in 1D
+        if p['ny']+1 == 3:
+            nc.variables['n'][:] = np.arange(1)
+            nc.variables['x'][:,:] = s['x'][0,:]
+            nc.variables['y'][:,:] = s['y'][0,:]
+        else:
+            nc.variables['n'][:] = np.arange(p['ny']+1)
+            nc.variables['x'][:,:] = s['x']
+            nc.variables['y'][:,:] = s['y']
+
         nc.variables['layers'][:] = np.arange(p['nlayers'])
         nc.variables['fractions'][:] = p['grain_size']
 
@@ -344,14 +358,21 @@ def append(outputfile, variables):
     # abort if netCDF4 is not available
     if not HAVE_NETCDF:
         return
-
+    
     with netCDF4.Dataset(outputfile, 'a') as nc:
         i = nc.variables['time'].shape[0]
         nc.variables['time'][i] = variables['time']
         for k, v in variables.items():
             if k == 'time':
                 continue
-            nc.variables[k][i,...] = v
+            # this is where a quasi 2D model is stored in 1D
+            if v.shape[0]==3: 
+                nc.variables[k][i,...] = v[1,:]
+            # this is where a 2D model is stored in 2D
+            else: 
+                nc.variables[k][i,...] = v
+            
+            
 
     set_bounds(outputfile)