debug Rubar3 (from scenario branch)

src/Solver/RubarBE.py

@@ -1,5 +1,5 @@
 # RubarBE.py -- Pamhyr
-# Copyright (C) 2024  INRAE
+# Copyright (C) 2024-2025  INRAE
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -20,7 +20,7 @@ import os
 import logging
 import numpy as np
 
-from tools import timer, trace, old_pamhyr_date_to_timestamp
+from tools import timer, trace, old_pamhyr_date_to_timestamp, logger_exception
 
 from Solver.CommandLine import CommandLineSolver
 
@@ -39,7 +39,7 @@ class Rubar3(CommandLineSolver):
         self._type = "rubar3"
 
         self._cmd_input = ""
-        self._cmd_solver = "@path @input -o @output"
+        self._cmd_solver = "@path @args @input"
         self._cmd_output = ""
 
     @classmethod
@@ -73,11 +73,12 @@ class Rubar3(CommandLineSolver):
             ("rubarbe_tf_4", "y"),
             ("rubarbe_tf_5", "y"),
             ("rubarbe_tf_6", "n"),
-            ("rubarbe_trased", "y"),
+            ("rubarbe_trased", "n"),
-            # ("rubarbe_optfpc", "0"),
+            ("rubarbe_optfpc", "0"),
-            # ("rubarbe_ros", "2650.0"),
+            # trased parameters
-            # ("rubarbe_dm", "0.1"),
+            ("rubarbe_ros", "2650.0"),
-            # ("rubarbe_segma", "1.0"),
+            ("rubarbe_dm", "0.1"),
+            ("rubarbe_segma", "1.0"),
             # Sediment parameters
             ("rubarbe_sediment_ros", "2650.0"),
             ("rubarbe_sediment_por", "0.4"),
@@ -129,11 +130,11 @@ class Rubar3(CommandLineSolver):
 
     def output_param(self):
         name = self._study.name
-        return f"{name}"
+        return f"profil.{name}"
 
     def log_file(self):
         name = self._study.name
-        return f"{name}"
+        return f"geomac-i.{name}"
 
     def export(self, study, repertory, qlog=None):
         self._study = study
@@ -143,10 +144,12 @@ class Rubar3(CommandLineSolver):
         self._export_ts(study, repertory, qlog, name=name)
         self._export_geomac_i(study, repertory, qlog, name=name)
         self._export_mail(study, repertory, qlog, name=name)
+        self._export_devers(study, repertory, qlog, name=name)
         self._export_condin(study, repertory, qlog, name=name)
         self._export_stricklers(study, repertory, qlog, name=name)
         self._export_hydro(study, repertory, qlog, name=name)
         self._export_condav(study, repertory, qlog, name=name)
+        # self._export_abshyd(study, repertory, qlog, name=name)
 
         return True
 
@@ -166,9 +169,8 @@ class Rubar3(CommandLineSolver):
 
             it = iter(params)
 
-            line = 0
+            param = next(it, None)
-            while line < 25:
+            while param is not None:
-                param = next(it)
                 name = param.name
                 value = param.value
 
@@ -196,12 +198,12 @@ class Rubar3(CommandLineSolver):
                         v2 = next(it).value
                         v3 = next(it).value
 
-                        f.write(f"{v2}{v3}")
+                        # f.write(f"{v2}{v3}")
 
                     # New line
                     f.write(f"\n")
 
-                line += 1
+                param = next(it, None)
 
     def _export_ts(self, study, repertory, qlog, name="0"):
         if qlog is not None:
@@ -330,6 +332,50 @@ class Rubar3(CommandLineSolver):
                     if ind % 3 != 0:
                         f.write("\n")
 
+    def _export_abshyd(self, study, repertory, qlog, name="0"):
+        if qlog is not None:
+            qlog.put("Export ABSHYD file")
+
+        with open(
+                os.path.join(
+                    repertory, f"abshyd.{name}"
+                ), "w+"
+        ) as f:
+            reach_ind = 1
+            for edge in study.river.enable_edges():
+                reach = edge.reach
+                lm = len(reach) + 1
+                f.write(f"{lm:>13}\n")
+
+                ind = 1
+                for mail in reach.inter_profiles_rk():
+                    f.write(f"{ind:>4} {mail:15.3f} {reach_ind:>4}\n")
+                    ind += 1
+
+                reach_ind += 1
+
+    def _export_devers(self, study, repertory, qlog, name="0"):
+        if qlog is not None:
+            qlog.put("Export DEVERS file")
+
+        with open(
+                os.path.join(
+                    repertory, f"devers.{name}"
+                ), "w+"
+        ) as f:
+            reach_ind = 1
+            for edge in study.river.enable_edges():
+                reach = edge.reach
+                lm = len(reach) + 1
+                f.write(f"{lm:>13}\n")
+
+                ind = 1
+                for mail in reach.inter_profiles_rk():
+                    f.write(f"{ind:>4} {0.0:15.3f} {0.0:15.3f} {1.0:>15.3f}\n")
+                    ind += 1
+
+                reach_ind += 1
+
     def _export_stricklers(self, study, repertory, qlog, name="0"):
         self._export_frot(study, repertory, qlog, name=name, version="")
         self._export_frot(study, repertory, qlog, name=name, version="2")
@@ -348,24 +394,57 @@ class Rubar3(CommandLineSolver):
                 lm = len(reach) + 1
                 f.write(f"{lm:>6}\n")
 
-                def get_stricklers_from_rk(rk):
+                lst = list(filter(
-                    return next(
+                    lambda f: f.is_full_defined(),
-                        map(
+                    edge.frictions.frictions
-                            lambda s: (
+                ))
-                                s.begin_strickler.medium if version == "2"
+
-                                else s.begin_strickler.minor
+                rk_min = 9999999.9
-                            ),
+                rk_max = -9999999.9
-                            filter(
+                coeff_min = -1.0
-                                lambda f: rk in f,
+                coeff_max = -1.0
-                                edge.frictions.lst
+                for s in lst:  # TODO optimise ?
-                            )
+                    if s.begin_rk > rk_max:
-                        )
+                        rk_max = s.begin_rk
+                        coeff_max = s.begin_strickler
+                    if s.begin_rk < rk_min:
+                        rk_min = s.begin_rk
+                        coeff_min = s.begin_strickler
+                    if s.end_rk > rk_max:
+                        rk_max = s.end_rk
+                        coeff_max = s.end_strickler
+                    if s.end_rk < rk_min:
+                        rk_min = s.end_rk
+                        coeff_min = s.end_strickler
+
+                def get_stricklers_from_rk(rk, lst):
+
+                    coeff = None
+                    if rk > rk_max:
+                        coeff = coeff_max
+                    elif rk < rk_min:
+                        coeff = coeff_min
+                    else:
+                        for s in lst:
+                            if (rk >= s.begin_rk and rk <= s.end_rk or
+                                    rk <= s.begin_rk and rk >= s.end_rk):
+                                coeff = s.begin_strickler  # TODO: inerpolate
+                                break
+
+                    # TODO interpolation if rk is not in frictons
+
+                    if coeff is None:
+                        logger.error(
+                            "Study frictions are not fully defined"
                         )
+                        return None
+
+                    return coeff.medium if version == "2" else coeff.minor
 
                 ind = 1
                 for mail in edge.reach.inter_profiles_rk():
-                    coef = get_stricklers_from_rk(mail)
+                    coef = get_stricklers_from_rk(mail, lst)
-
+                    if coef is not None:
                         f.write(f"{ind:>6} {coef:>12.5f}")
 
                     ind += 1
@@ -383,7 +462,7 @@ class Rubar3(CommandLineSolver):
             for edge in study.river.enable_edges():
                 reach = edge.reach
 
-                f.write(f"0.0\n")
+                f.write(f"{0.0}\n")
 
                 ics = study.river.initial_conditions.get(edge)
                 data = self._export_condin_init_data(ics)
@@ -392,7 +471,9 @@ class Rubar3(CommandLineSolver):
                 first = profiles[0]
                 last = profiles[-1]
 
+                # if first not in data or last not in data:
                 if first.rk not in data or last.rk not in data:
+                    print(data)
                     logger.error(
                         "Study initial condition is not fully defined"
                     )
@@ -493,8 +574,103 @@ class Rubar3(CommandLineSolver):
 
             for bc in bcs:
                 f.write(f"{len(bc)}\n")
+                if bc.bctype == "ZD":
                     for d0, d1 in bc.data:
                         f.write(f"{d1} {d0}\n")
+                else:
+                    for d0, d1 in bc.data:
+                        f.write(f"{d0} {d1}\n")
+
+    def read_profil(self, study, fname, results, qlog=None, name="0"):
+        logger.info(f"read: profil.{name}")
+        with open(fname, "r") as f:
+            reachs = []
+            for i, edge in enumerate(study.river.enable_edges()):
+                reach = edge.reach
+                # Add results reach to reach list
+                res_reach = results.river.add(i)
+                reachs.append((res_reach, len(reach) + 1))
+
+            def read_data_line(f):
+                line = f.readline().split()
+                if len(line) > 3:
+                    return tuple(map(float, line))
+                else:
+                    return (0.0, 0.0, 0.0, 0.0)
+
+            def set_and_compute_limites(reach, ind, z, q, s):
+                reach.set(ind, timestamp, "Z", z)
+                reach.set(ind, timestamp, "Q", q)
+                reach.set(ind, timestamp, "V", s)
+
+                profile = reach.profile(ind)
+                limits = profile.geometry.get_water_limits(z)
+                reach.set(
+                    ind, timestamp, "water_limits", limits
+                )
+
+            ts = set()
+            end = False
+            while True:
+                line = f.readline()
+                if line == "":
+                    results.set("timestamps", ts)
+                    return
+
+                timestamp = float(line)
+                ts.add(timestamp)
+
+                logger.info(f"read: profil.{name}: timestamp = {timestamp}")
+
+                for reach, lm in reachs:
+                    # First profile
+                    h, s, q, z = read_data_line(f)
+                    set_and_compute_limites(reach, 0, z+h, q, s)
+
+                    # For each profile
+                    ind = 1
+                    ph, ps, pq, pz = read_data_line(f)
+                    while ind < lm - 2:
+                        h, s, q, z = read_data_line(f)
+                        set_and_compute_limites(
+                            reach, ind,
+                            (pz + z + ph + h) / 2,
+                            (pq + q) / 2,
+                            (ps + s) / 2
+                        )
+
+                        ph = h
+                        ps = s
+                        pq = q
+                        pz = z
+                        ind += 1
+
+                    # Last profile
+                    h, s, q, z = read_data_line(f)
+                    set_and_compute_limites(reach, ind, z+h, q, s)
+
+    @timer
+    def results(self, study, repertory, qlog=None, name="0"):
+        results = Results(
+            study=study,
+            solver=self,
+            repertory=repertory,
+            name=name,
+        )
+        results_file = f"profil.{name}"
+
+        fname = os.path.join(repertory, results_file)
+        if not os.path.isfile(fname):
+            logger.warning(f"Result file {results_file} does not exist")
+            return None
+        try:
+            self.read_profil(study, fname, results, qlog, name=name)
+        except Exception as e:
+            logger.error(f"Failed to read results")
+            logger_exception(e)
+            return None
+
+        return results
 
 
 class RubarBE(Rubar3):
@@ -506,5 +682,5 @@ class RubarBE(Rubar3):
         self._type = "rubarbe"
 
         self._cmd_input = ""
-        self._cmd_solver = "@path @input -o @output"
+        self._cmd_solver = "@path @args @input"
         self._cmd_output = ""

src/View/MainWindow.py

@@ -30,6 +30,7 @@ from platformdirs import user_cache_dir
 
 from Solver.AdisTS import AdisTS
 from Solver.Mage import Mage8
+from Solver.RubarBE import Rubar3
 from tools import logger_exception
 
 from PyQt5 import QtGui
@@ -1414,7 +1415,10 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
         if type(results) is str:
             logger.info(f"Open results from {os.path.dirname(results)}")
 
+            if ".BIN" in results:
                 name = os.path.basename(results).replace(".BIN", "")
+            elif "profil." in results:
+                name = os.path.basename(results).replace("profil.", "")
 
             def reading_fn():
                 self._tmp_results = solver.results(
@@ -1545,7 +1549,7 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
         dialog.setFileMode(QFileDialog.FileMode.ExistingFile)
         dialog.setDefaultSuffix(".BIN")
         # dialog.setFilter(dialog.filter() | QtCore.QDir.Hidden)
-        dialog.setNameFilters(['Mage (*.BIN)'])
+        dialog.setNameFilters(['Mage (*.BIN)', 'Rubar3 (profil.*)'])
 
         if self._study.filename is not None:
             if self._last_solver is None:
@@ -1560,8 +1564,15 @@ class ApplicationWindow(QMainWindow, ListedSubWindow, WindowToolKit):
         if dialog.exec_():
             file_name = dialog.selectedFiles()
             logger.info(f"Select results: {file_name}")
+
+            solver = None
+            if ".BIN" in file_name:
+                solver = Mage8("Mage8")
+            else:
+                solver = Rubar3("Rubar3")
+
             self.open_solver_results(
-                Mage8("Mage"),
+                solver,
                 results=file_name[0]
             )
 

src/View/Results/PlotAC.py

@@ -116,7 +116,7 @@ class PlotAC(PamhyrPlot):
 
         self.collection = self.canvas.axes.fill_between(
             x, z, water_z,
-            where=z <= water_z,
+            where=list(map(lambda x: x <= water_z, z)),
             color=self.color_plot_river_water_zone,
             alpha=0.7, interpolate=True
         )
@@ -263,7 +263,7 @@ class PlotAC(PamhyrPlot):
         self.collection.remove()
         self.collection = self.canvas.axes.fill_between(
             x, z, water_z,
-            where=z <= water_z,
+            where=list(map(lambda x: x <= water_z, z)),
             color=self.color_plot_river_water_zone,
             alpha=0.7, interpolate=True
         )

src/View/RunSolver/Window.py

@@ -486,7 +486,8 @@ class SolverLogWindow(PamhyrWindow):
         if self._results is None:
             def reading_fn():
                 self._results = self._solver.results(
-                    self._study, self._workdir, qlog=self._output
+                    self._study, self._workdir,
+                    qlog=self._output, name=self._study.name
                 )
 
             dlg = ReadingResultsDialog(reading_fn=reading_fn, parent=self)