Comput loop fix

smoderp2d/core/surface.py

@@ -3,6 +3,7 @@
 
 import numpy as np
 import numpy.ma as ma
+import numpy.ma as ma
 
 from smoderp2d.core.general import Globals, GridGlobals
 if Globals.isStream:
@@ -167,19 +168,22 @@ def return_str_vals(self, i, j, sep, dt, extra_out):
                 )
                 bil_ = ''
             else:
+                h_sheet = min(arr.h_total_new[i,j],arr.h_crit[i,j])
+                h_rill = max(arr.h_total_new[i,j]-arr.h_crit[i,j],0)
                 velocity = ma.where(
                     arr.h_sheet == 0,
                     0,
-                    arr.vol_runoff / dt / (arr.h_sheet*GridGlobals.dx)
+                    arr.vol_runoff / dt / (h_sheet*GridGlobals.dx)
                 )
+
                 # if profile1d provider - the data in extra output are the unit
                 #                          width data
                 #                     if you need runoff from non-unit slope and
                 #                     with extra output calculate it yourself
                 line = '{0:.4e}{sep}{1:.4e}{sep}{2:.4e}{sep}{3:.4e}{sep}' \
                        '{4:.4e}{sep}{5:.4e}{sep}{6:.4e}{sep}{7:.4e}{sep}' \
                        '{8:.4e}{sep}{9:.4e}'.format(
-                    arr.h_sheet[i, j],
+                    h_sheet,
                     vol_runoff / dt[i, j],
                     vol_runoff,
                     velocity[i, j],
@@ -196,7 +200,7 @@ def return_str_vals(self, i, j, sep, dt, extra_out):
                     line += '{sep}{0:.4e}{sep}{1:.4e}{sep}{2:.4e}{sep}{3:.4e}' \
                             '{sep}{4:.4e}{sep}{5:.4e}{sep}{6:.4e}{sep}' \
                             '{7:.4e}'.format(
-                        arr.h_rill[i, j],
+                        h_rill,
                         arr.rillWidth[i, j],
                         vol_runoff_rill / dt[i, j],
                         vol_runoff_rill,
@@ -231,7 +235,7 @@ def __runoff(sur, dt, effect_vrst, ratio):
 
     :return: TODO
     """
-    h_total_pre = sur.h_total_pre
+    h_total_pre = ma.copy(sur.h_total_new)
     h_crit = sur.h_crit
     state = sur.state  # da se tady podivat v jakym jsem casovym kroku a jak
     # se a

smoderp2d/courant.py

@@ -22,7 +22,7 @@ def __init__(self):
         """TODO."""
         self.cour_speed = 0
         # citical courant value
-        self.cour_crit = 0.95
+        self.cour_crit = 0.3
         self.cour_most = self.cour_crit + 1
         self.cour_most_rill = self.cour_crit + 1.0
         self.cour_coef = 0.5601
@@ -97,7 +97,6 @@ def CFL(self, v, delta_t, effect_cont, co, rill_courant):
         :param rill_courant: TODO
         """
         cour = v / self.cour_coef * delta_t / effect_cont
-        cour = ma.maximum(cour, rill_courant)
         if ma.any(cour > self.cour_most):
             self.i = np.unravel_index(ma.argmax(cour), cour.shape)[0]
             self.j = np.unravel_index(ma.argmax(cour), cour.shape)[1]

smoderp2d/io_functions/hydrographs.py

@@ -112,7 +112,7 @@ def __init__(self):
                               '{linesep}'.format(sep=SEP, linesep=os.linesep)
                 else:
                     header += 'time[s]{sep}deltaTime[s]{sep}rainfall[m]{sep}'\
-                              'waterLevel[m]{sep}sheetFlow[m3/s]{sep}' \
+                              'sheetWaterLevel[m]{sep}sheetFlow[m3/s]{sep}' \
                               'sheetVRunoff[m3]{sep}sheetFlowVelocity[m/s]' \
                               '{sep}sheetVRest[m3]{sep}infiltration[m]{sep}' \
                               'surfaceRetetion[m]{sep}state{sep}vInflow[m3]' \

smoderp2d/processes/surface.py

@@ -2,4 +2,5 @@
 
 
 def shallowSurfaceKinematic(a, b, h_sheet):
+    a = 160
     return ma.power(h_sheet, b) * a

smoderp2d/runoff.py

@@ -11,12 +11,14 @@
 
 """
 
+import sys
 import time
 import numpy as np
 import numpy.ma as ma
 
 from smoderp2d.core.general import Globals, GridGlobals
 from smoderp2d.core.vegetation import Vegetation
+from smoderp2d.core.surface import runoff
 from smoderp2d.core.surface import get_surface
 from smoderp2d.core.subsurface import Subsurface
 from smoderp2d.core.cumulative_max import Cumulative
@@ -264,14 +266,22 @@ def run(self):
             self.delta_t, mask=GridGlobals.masks
         )
 
+        time_ = time.time()
         while ma.any(self.flow_control.compare_time(end_time)):
 
+
+
             self.flow_control.save_vars()
             self.flow_control.refresh_iter()
 
             # iteration loop
+            print ('=++++++++++++++++++++++++++++++')
+            timepreiter = time.time()
             while self.flow_control.max_iter_reached():
 
+                viter1 = time.time()
+                #print ('                        ')
+
                 self.flow_control.update_iter()
                 self.flow_control.restore_vars()
 
@@ -288,6 +298,47 @@ def run(self):
                     self.courant
                 )
 
+
+
+                viter2 = time.time()
+                #print ('dt            {}'.format(self.delta_t[6,2]))
+                #print ('time          {}'.format(self.flow_control.total_time[6,2]))
+                #print ('iter          {}'.format(self.flow_control.iter_))
+                #print ('hcir          {}'.format(self.surface.arr.h_crit[6,2]))
+                #print ('courant       {}'.format(self.courant.cour_most))
+                #print ('cour_speed    {}'.format(self.courant.cour_speed))
+
+                # Calculate actual rainfall and adds up interception todo:
+                # AP - actual is not storred in hydrographs
+                actRain = self.time_step.do_next_h(
+                    self.surface,
+                    self.subsurface,
+                    self.rain_arr,
+                    self.cumulative,
+                    self.hydrographs,
+                    self.flow_control,
+                    self.courant,
+                    potRain,
+                    self.delta_t
+                )
+
+                viter3 = time.time()
+                runoff_return = runoff(
+                    self.surface.arr, self.delta_t,
+                    Globals.get_mat_effect_cont(), self.flow_control.ratio
+                )
+
+
+                viter4 = time.time()
+                v = ma.maximum(runoff_return[0], runoff_return[1])
+                self.courant.CFL(
+                    v, self.delta_t,
+                    Globals.get_mat_effect_cont(),
+                    '---',
+                    None
+                )
+                # print ('courant       {}'.format(self.courant.cour_most))
+
                 # stores current time step
                 delta_t_tmp = self.delta_t
 
@@ -311,21 +362,26 @@ def run(self):
                     ma.logical_and(delta_t_tmp == self.delta_t,
                                    self.flow_control.compare_ratio())
                 ):
+                    potRain = self.time_step.do_flow(
+                        self.surface,
+                        self.subsurface,
+                        self.delta_t,
+                        self.flow_control,
+                        self.courant
+                    )
+                    viter5 = time.time()
+                    print (viter2-viter1, viter3-viter2, viter4-viter3,
+                            viter5-viter4)
                     break
+                viter5 = time.time()
+                print (viter2-viter1, viter3-viter2, viter4-viter3,
+                        viter5-viter4)
 
-            # Calculate actual rainfall and adds up interception todo:
-            # AP - actual is not storred in hydrographs
-            actRain = self.time_step.do_next_h(
-                self.surface,
-                self.subsurface,
-                self.rain_arr,
-                self.cumulative,
-                self.hydrographs,
-                self.flow_control,
-                self.courant,
-                potRain,
-                self.delta_t
-            )
+            timepostiter = time.time()
+
+            print ('iter          {}'.format(self.flow_control.iter_))
+            print ('dt            {}'.format(self.delta_t[6,2]))
+            print ('vsechny iterace   {}'.format(timepostiter - timepreiter))
 
             # if the iteration exceed the maximal amount of iteration
             # last results are stored in hydrographs
@@ -425,6 +481,23 @@ def run(self):
                 self.surface.arr.state
             )
 
+            self.cumulative.update_cumulative(
+                self.surface.arr,
+                self.subsurface.arr,
+                self.delta_t)
+            self.hydrographs.write_hydrographs_record(
+                None,
+                None,
+                self.flow_control,
+                self.courant,
+                self.delta_t,
+                self.surface,
+                self.cumulative,
+                actRain)
+
+            #print (self.surface.arr.h_total_pre[6,2])
+            #print (self.surface.arr.h_total_new[6,2])
+
             self.surface.arr.h_total_pre = ma.copy(self.surface.arr.h_total_new)
 
             timeperc = 100 * (self.flow_control.total_time + self.delta_t) / end_time
@@ -441,6 +514,10 @@ def run(self):
             # proceed to next time
             self.flow_control.update_total_time(self.delta_t)
 
+            print ('jeden krok v case {}'.format(time.time() - time_))
+            time_ = time.time()
+            if self.flow_control.total_time[6,2] > 300 : sys.exit()
+
     def save_output(self):
         """TODO."""
         Logger.info('Saving output data...')

smoderp2d/time_step.py

@@ -51,6 +51,7 @@ def do_flow(surface, subsurface, delta_t, flow_control, courant):
         runoff_return = runoff(
             surface.arr, delta_t, mat_effect_cont, fc.ratio
         )
+
 
         cond_state_flow = surface_state > Globals.streams_flow_inc
         v_sheet = ma.where(cond_state_flow, 0, runoff_return[0])
@@ -93,11 +94,13 @@ def do_flow(surface, subsurface, delta_t, flow_control, courant):
             cond_state_flow, surface.arr.vel_rill, runoff_return[11]
         )
 
+        #print ('vel_sheet {}'.format(v_sheet[6,2]))
+        #print ('vel_rill {}'.format(v_rill[6,2]))
+
         v = ma.maximum(v_sheet, v_rill)
-        co = 'sheet'
+        co = '---'
         courant.CFL(
-            v,
-            delta_t,
+            v, delta_t,
             mat_effect_cont,
             co,
             rill_courant
@@ -248,18 +251,10 @@ def do_next_h(self, surface, subsurface, rain_arr, cumulative, hydrographs,
         subsurface.fill_slope()
         """
 
-        cumulative.update_cumulative(
-            surface.arr,
-            subsurface.arr,
-            delta_t)
-        hydrographs.write_hydrographs_record(
-            None,
-            None,
-            flow_control,
-            courant,
-            delta_t,
-            surface,
-            cumulative,
-            actRain)
+        #print ('h_sheet       {}'.format(surface.arr.h_sheet[6,2]))
+        #print ('h_rill        {}'.format(surface.arr.h_rill[6,2]))
+        #print ('h_total_pre   {}'.format(surface.arr.h_total_pre[6,2]))
+        #print ('h_total_new   {}'.format(surface.arr.h_total_new[6,2]))
+        # input('press')
 
         return actRain

utils/hydragramy.R

@@ -10,8 +10,7 @@ library('manipulate')
 # install package is missing with: install.packages("manipulate")
 #
 # root dir
-root  <-  "d:/0_Smoderp/00_QGtest_ds_plocha/out2"
-root  <-  "d:/0_Smoderp/02_AGPro_provider/out19_AG"
+root  <-  c("tests/data/output/", 'tests/data/output-PR321/')
 
 root  <-  "tests/data/output/"
 #root  <-  "d:/2_granty_projekty/2_Bezici/2022_RAGO/01_reseni_projektu/00_test_Smoderp/out2"
@@ -21,8 +20,8 @@ outdir <- 'control_point'
 # point000.dat -> id = 1
 # point001.dat -> id = 2
 # atd...
-id1_ = 1
-id2_ = 3
+id1_ = 1 + 6
+id2_ = id1_ 
 #2+6;1+4
 # End setting  
 #
@@ -69,14 +68,16 @@ pp = function(t1,t2,sel,add_,sel2,od,do,stejny,titles)
   plot(t1[,1],t1[[sel]],
        ylab = '',type = 'o',lwd=2,xlim = c(od,do),ylim=r1,cex=0.5)
   grid()
-  mtext(paste(basename(titles[1]),":",sel),side = 3,line = 0.8,adj = 0,cex = 1.5)
+  # mtext(paste(basename(titles[1]),":",sel),side = 3,line = 0.8,adj = 0,cex = 1.5)
+  mtext(paste(titles[1],":",sel),side = 3,line = 0.8,adj = 0,cex = 1.5)
   mtext(names1_[sel],side = 2,line = 3)
   if (add_) {
     par(new=TRUE)
     plot(t2[,1],t2[[sel2]],
          axes = FALSE, ylab = '',type = 'o',col=2,lwd=2,xlim = c(od,do),ylim=r2,cex=0.5)
     axis(4,col.ticks = 2, col = 2,col.axis=2)
-    mtext(paste(basename(titles[2]),":",sel2),side = 3,line = 2,adj = 1,cex = 1.5, col=2)
+    # mtext(paste(basename(titles[2]),":",sel2),side = 3,line = 2,adj = 1,cex = 1.5, col=2)
+    mtext(paste(titles[2],":",sel2),side = 3,line = 2,adj = 1,cex = 1.5, col=2)
     mtext(names2_[sel2],side = 4,line = 3,col = 2)
   }
 }
@@ -97,9 +98,9 @@ plot_ = function(id1,id2,title='')
              # sel = slider(initial = 5,1,n1,label = 'spoupec v levem grafu'),
              sel = picker(as.list(names1_), initial = 'wLevelTotal.m.'),
              add_= checkbox(TRUE,'pridat druhy graf'),
-             stejny= checkbox(FALSE,'stejny meritka'),
+             stejny= checkbox(TRUE,'stejny meritka'),
              # sel2 = slider(initial = n2, 1,n2,label = 'spoupec v pravem grafu'),
-             sel2 = picker(as.list(names2_),initial = 'wLevelTotal.m.'),
+             sel2 = picker(as.list(names2_), initial = 'wLevelTotal.m.'),
              od = slider(initial = 0     ,0,maxCas,label = 'cas od'),
              do = slider(initial = maxCas,0,maxCas,label = 'cas do')
              )