Merge branch 'master' of https://github.com/rustychris/stompy

Author: Rusty Holleman

examples/openfoam.ipynb

@@ -15,6 +15,7 @@
 import copy
 import xarray as xr
 import numpy as np
+from .. import utils
 from . import unstructured_grid
 import logging
 log=logging.getLogger('multi_ugrid')
@@ -417,7 +418,7 @@ def __init__(self,paths,cleanup_dfm='auto',xr_kwargs={},grid=None,
 
     def load(self,**xr_kwargs):
         if isinstance(self.paths[0],str):
-            return [xr.open_dataset(p,**xr_kwargs) for p in self.paths]
+            return [xr.open_dataset(p,**xr_kwargs) for p in utils.progress(self.paths)]
         elif isinstance(self.paths[0],xr.Dataset):
             return self.paths
         else:

stompy/grid/multi_ugrid.py

@@ -14,7 +14,7 @@
 
 def dd_to_synonyms(code):
     # print "Looking for synonyms of ",code
-    m = re.match('\d\d_(.+)',code) 
+    m = re.match(r'\d\d_(.+)',code) 
     if m:
         short_code = m.group(1)
         # print "Got a short code ",short_code

stompy/io/rdb_datadescriptors.py

@@ -87,12 +87,10 @@ def build_process_db(self):
             shutil.rmtree(self.proc_table_dst_dir)
         shutil.copytree(self.proc_table_src_dir,self.proc_table_dst_dir)
 
-        print("First call to waqpb_export")
+        print("First call to waqpb_export - creating procesm.asc from decomposed tables")
         sys.stdout.flush()
         output = utils.call_with_path(self.waqpbexport,self.proc_table_dst_dir).decode('latin1') 
         print("Suppressed output")
-        # That says Normal end, but then goes on to
-        # make proces.asc
         sys.stdout.flush()
         print("First call to waqpb_export DONE")
         sys.stdout.flush()
@@ -101,6 +99,7 @@ def build_process_db(self):
         # First line of proces_asc includes a process count. Rewrite that line, shove in
         # our new process, and copy the rest of the file.
 
+        print(f"Transcribing procesm.asc -> proces.asc, adding {len(self.custom_procs)} custom processes")
         procesm_asc=os.path.join(self.proc_table_dst_dir,'procesm.asc')
         proces_asc=os.path.join(self.proc_table_dst_dir,'proces.asc')
 
@@ -116,18 +115,22 @@ def build_process_db(self):
                     fp.write(line)
         sys.stdout.flush()
 
-        print("Calling waqpb_import")
+        print("Calling waqpb_import to convert proces.asc back to tables")
         sys.stdout.flush()
-        # Pretty sure this is the one that fails.
+        # This is the one likely to fail.
         output = utils.call_with_path(self.waqpbimport,self.proc_table_dst_dir).decode('latin1') 
         print("Suppressed output")
+        #print(output) # uncomment this line to diagnose Fatal Error 
         sys.stdout.flush()
         print("Return from waqpb_import")
         sys.stdout.flush()
-        print("Second call to waqpb_export")
+
+        print("Second call to waqpb_export to take the updated tables and write compiled form of tables")
         sys.stdout.flush()
         output = utils.call_with_path(self.waqpbexport,self.proc_table_dst_dir).decode('latin1') 
         print( "Suppressed output")
+        #print("Output from waqpb_export")
+        #print(output)
         sys.stdout.flush()
         print("Return from second call to waqpb_export")
         sys.stdout.flush()
@@ -285,3 +288,62 @@ def fall_velocity(self, **kw):
 """
         return process
 
+    def custom_DynDen(self,**kw):
+        idx=self.copy_count['DynDen']
+        self.copy_count['DynDen']+=1
+        if self.copy_count['DynDyn']>1:
+            raise Exception("custom_DynDen can only be used once")
+        
+        p=dict(name="DynDen")
+        p.update(kw)
+
+        # Presumably want to enable it, though could be optional
+        if isinstance(self,dflow_model.DFlowModel):
+            self.dwaq.add_process(p['name'])
+        else: # WaqModel or related
+            self.add_process(p['name'])
+        
+        # With SW_Uitz=0.0:
+        # SD = PAConstant / ExtVl
+        process=f"""
+{p['name']:10}                    Reuse Secchi depth for denit rate
+SECCHI    ; module name
+123       ; TRswitch
+        22; # input items for segments
+InvDenRate     -999.000     x POC                                                    (1/m)               
+IM1             0.00000       inorganic matter (IM1)                                 (gDM/m3)            
+IM2             0.00000       inorganic matter (IM2)                                 (gDM/m3)            
+IM3             0.00000       inorganic matter (IM3)                                 (gDM/m3)            
+POC1            0.00000       POC1 (fast decomposing fraction)                       (gC/m3)             
+POC2            0.00000       POC2 (medium decomposing fraction)                     (gC/m3)             
+POC3            0.00000       POC3 (slow decomposing fraction)                       (gC/m3)             
+POC4            0.00000       POC4 (particulate refractory fraction)                 (gC/m3)             
+ExtVlODS        0.00000       VL extinction by DOC                                   (1/m)               
+Chlfa           0.00000       Chlorophyll-a concentration                            (mg/m3)             
+SW_Uitz         0.00000       Extinction by Uitzicht On (1) or Off (0)               (-)                 
+UitZDEPT1       1.20000       Z1 (depth)                                             (m)                 
+UitZDEPT2       1.00000       Z2 (depth)                                             (m)                 
+UitZCORCH       2.50000       CORa correction factor                                 (-)                 
+UitZC_DET      0.260000E-01   C3 coeff. absorption ash weight & detritus             (-)                 
+UitZC_GL1      0.730000       C1 coeff. absorption ash weight & detritus             (-)                 
+UitZC_GL2       1.00000       C2 coeff. absorption ash weight & detritus             (-)                 
+UitZHELHM      0.140000E-01   Hel_h constant                                         (1/nm)              
+UitZTAU         7.80000       Tau constant calculation transparency                  (-)                 
+UitZangle       30.0000       Angle of incidence solar radiation                     (degrees)           
+DMCFDetC        2.50000       DM:C ratio DetC                                        (gDM/gC)            
+DetCS1          1.70000     x Poole-Atkins constant                                  (-)                 
+         0; # input items for exchanges
+         1; # output items for segments
+RcDenSed                    x 1st order denit m/d                                    (m)                 
+         0; # output items for exchanges
+         1; # fluxes
+dDumDynDen                  x dummy flux to access Secchi                            (-)                 
+         2; # stoichiometry lines
+IM1         dDumDynDen     0.00000
+POC1        dDumDynDen     0.00000
+         0; # stoichiometry lines dispersion arrays
+         0; # stoichiometry lines velocity arrays
+END
+"""
+        self.custom_procs.append(process)
+

stompy/model/delft/custom_process.py

@@ -104,7 +104,7 @@ def decstrip(s):
 
     time0,time_unit = parse_time0(time_meta)
     times=time0 + time_unit*frames['tsec']
-    ds['time']=( ('time',), times)
+    ds['time']=( ('time',), times.astype('M8[ns]'))
     ds['tsec']=( ('time',), frames['tsec'])
 
     region_names=[decstrip(s) for s in regions['name']]

stompy/model/delft/io.py

@@ -1040,7 +1040,7 @@ def check_volume_conservation_incr(self,seg_select=slice(None),
                     ('type','S20')]
     @property
     def n_boundaries(self):
-        return -self.pointers[:,:2].min()
+        return max(0,-self.pointers[:,:2].min())
 
     # for a while, was using 'grouped', but that relies on boundary segments,
     # which didn't transfer well to aggregated domains which usually have many
@@ -1737,7 +1737,7 @@ def trav(c,mark):
                 ngroups+=1
         return groups
 
-    def extract_transect_flow(self,transect,func=False,time_range=None):
+    def extract_transect_flow(self, transect, func=False, time_range=None, quantity='Q'):
         """
         Extract time series of discharge through a transect as a function of time.
         transect is expected to have the same structure as it is configured on
@@ -1749,15 +1749,24 @@ def extract_transect_flow(self,transect,func=False,time_range=None):
         time_range: [np.datetime64,np.datetime64]
 
         if func is True, instead return a function that takes a datetime64 and returns flow.
+
+        quantity: Q or area, determines the variable to extract
         """
         exchs=transect[1]
 
         def fn(t,exchs=exchs,hydro=self):
             t_secs = (t - np.datetime64(hydro.time0))/np.timedelta64(1,'s')
-            flo = hydro.flows(t_secs)
-            signs=np.sign(exchs)
             exch0=np.abs(exchs)-1
-            return np.sum(flo[exch0]*signs)
+            signs=np.sign(exchs)
+
+            if quantity=='Q':
+                flo = hydro.flows(t_secs)
+                return np.sum(flo[exch0]*signs)
+            elif quantity=='area':
+                area = hydro.areas(t_secs)
+                return np.sum(area[exch0])
+            else:
+                raise Exception("Only Q and area are supported outputs")
         if func:
             return fn
 
@@ -1772,15 +1781,20 @@ def fn(t,exchs=exchs,hydro=self):
         t_secs=self.t_secs[tidx_start:tidx_stop]
         times=np.datetime64(self.time0)+t_secs*np.timedelta64(1,'s')
 
-        Q=np.zeros( len(t_secs), np.float64)
-        for i,t in enumerate(times):
-            Q[i] = fn(t)
+        data=np.zeros( len(t_secs), np.float64)
+        for i,t in utils.progress(enumerate(times)):
+            data[i] = fn(t)
 
-        da= xr.DataArray(data=Q, dims=['time'], 
+        da= xr.DataArray(data=data, dims=['time'], 
                          coords=dict( time=times,
                                       time_seconds=("time",t_secs) ),
-                         attrs=dict(transect=transect[0], units='m3 s-1'))
-        da.name="discharge"
+                         attrs=dict(transect=transect[0]))
+        if quantity=='Q':
+            da.name="discharge"
+            da.attrs['units']='m3 s-1'
+        elif quantity=='area':
+            da.name="area"
+            da.attrs['units']='m2'
         return da
 
     # Data formats on disk
@@ -2447,7 +2461,12 @@ def write_geom(self):
 
         dest=self.flowgeom_filename
         if os.path.exists(dest) or os.path.lexists(dest):
-            assert not os.path.samefile(dest,orig)
+            # RH: originally checked for not samefile
+            # then it was commented out, but that seems dangerous
+            # If this check is somehow problematic update comments here
+            # to explain
+            assert dest != orig
+#            assert not os.path.samefile(dest,orig)
             if self.overwrite:
                 self.log.warning("Removing old geom file %s"%dest)
                 os.unlink(dest)

stompy/model/delft/waq_scenario.py

@@ -0,0 +1,54 @@
+from . import depth_average
+from ... import utils
+import os
+
+def main(argv=None):
+    import argparse
+    
+    parser = argparse.ArgumentParser(description='Postprocess OpenFOAM results, esp. depth-averaging')
+
+    parser.add_argument("-c", "--case", help="Path to case, default is working directory",
+                        default=".")
+    parser.add_argument("-U", "--velocity", help="Enable depth-averaged velocity output",
+                        action="store_true")
+    parser.add_argument("-s","--stage", help="Enable stage (water surface elevation) output",
+                        action="store_true")
+    parser.add_argument("-t","--time", help="Select output times as comma-separated time names, or 'all'",
+                        default="all")
+    parser.add_argument("-r","--res", help="Output raster cell size (positive), or count in x (negative)",
+                        default=20.0, type=float)
+    parser.add_argument("-f","--force", help="Force overwrite existing files",action='store_true')
+
+    args = parser.parse_args(argv)
+
+    pf = depth_average.PostFoam(sim_dir=args.case)
+    
+    if args.time=='all':
+        times = pf.available_times(omit_zero=True)
+    else:
+        times = args.time.split(",")
+    print(f"{len(times)} time step(s) to process")
+
+    raster_info = pf.set_raster_parameters(dx=args.res)
+
+    print(f"Output rasters will be {raster_info['ny']}x{raster_info['nx']}")
+    output_path=os.path.join(args.case,"postProcessing",'python')
+    if not os.path.exists(output_path):
+        os.makedirs(output_path)
+    print(f"Writing output to {output_path}")
+        
+    for t in utils.progress(times,msg="Time steps: %s"):
+        if args.velocity:
+            fn=os.path.join(output_path,f"U_{t}.tif")
+            if args.force or not os.path.exists(fn):
+                fld_U=pf.to_raster('U',t)
+                # should put the velocity components into bands.
+                fld_U.write_gdal(fn,overwrite=True)
+        if args.stage:
+            fn=os.path.join(output_path,f"stage_{t}.tif")
+            if args.force or not os.path.exists(fn):
+                fld_stage=pf.to_raster('wse',t)
+                fld_stage.write_gdal(fn,overwrite=True)
+        
+if __name__ == '__main__':
+    main()