Skip to content

Hydraulic structure - #524

Open
maximejay wants to merge 8 commits into
DassHydro:mainfrom
maximejay:barrage
Open

maximejay wants to merge 8 commits into
DassHydro:mainfrom
maximejay:barrage

Conversation

@maximejay

Copy link
Copy Markdown
Collaborator

Implémentation de structures hydrauliques dans smash:

Deux ouvrages sont actuellement codées :

  • les barrages
  • des injections d'eau

Dans le code fortran:

  • Nouveau type dérivé : hydraulic_structureDT
  • Nouveaux sous type dérivés : DamDT et InflowDT

Ces type dérivés permettent de stocker des data relative aux fonctionnements des ouvrages hydrauliques dans input_dataDT. D'autres structure peuvent être développé facilement sur ce principe (prélèvement, écrêtement)

Le calcul lié à ces ouvrages est réalisé dans la routine de routage, pour l'instant uniquement pour le module lr !! . Une nouvelle variable d'états permet de stocker la cote d'eau des barrages

Au niveau du maillage:

  • chaque ouvrage est défini comme un exutoire
  • une nouvelle liste dans generate_mesh permet de définir sir un exutoire est: un exutoire, une jauge, un barrage, une injection d'eau ...

A prévoir:

  • ne pas lire les débits si l'exutoire est typé autre que 'jauge'
  • standardiser les données relatives aux barrages et aux injection d'eau pour permettre leurs lectures automatiques dans smash (et aux future ouvrage) : necessite rélfexion commune.
  • Comment sortir le débit avant ouvrage si un exutoire est défini, et après l'ouvrage si une structure hydraulique est définie ? utilitée ?
  • Doublons dans le tableau response.q jauge/structurehydraulique si des jauge et des ouvrages sont identifiés au memes coordonnées ? est-ce problématique ?
  • Documentation
  • Tests: nécessite de regénérer la baseline car nouvelle variable d'état. Pb test de prcp_indices qui ne passe pas => Raison inconnu (on dirait que les valeurs de np.nan ne passe pas le test !)

Pour tester le code voici un script d'exemple:

import smash
import numpy as np
from smash.factory import load_dataset, generate_mesh
import matplotlib.pyplot as plt
import pandas as pd
import os
import smashbox as sb #seulement pour les graphiques et pour lire les paramètres

flwdir = load_dataset("flwdir")

mesh = generate_mesh(
    flwdir_path=flwdir,
    x=[778306.4, 778306.4],
    y=[6349743, 6349743],
    area=[115 * 1e6, 115 * 1e6],
    code=["V7155035", "Barrage_StCecile"],
    outlet_type=["gauge", "dam"],
)

mesh_no_dam = generate_mesh(
    flwdir_path=flwdir,
    x=[778306.4],
    y=[6349743],
    area=[115 * 1e6],
    code=["V7155035"],
    outlet_type=["gauge"],
)

setup = {}
setup["hydrological_module"] = "gr4"
setup["routing_module"] = "lr"
setup["qobs_directory"] = "/home/maxime/DATA/QOBS_60min"
setup["prcp_directory"] = "/home/maxime/DATA/PLUIE"
setup["pet_directory"] = "/home/maxime/DATA/ETP-SFR-FRA-INTERA_L93"
setup["prcp_conversion_factor"] = 0.1
setup["read_prcp"] = True
setup["read_pet"] = True
setup["read_qobs"] = True
setup["daily_interannual_pet"] = True

# période de chauffe
setup["start_time"] = "2014-01-01 00:00"
setup["end_time"] = "2015-10-27 00:00"
# setup["end_time"] = "2014-02-27 00:00"

model = smash.Model(setup, mesh)
model_no_dam = smash.Model(setup, mesh_no_dam)

smash.io.read_grid_parameters(
    model, path="/home/maxime/DEV/smashbox/smashbox/asset/params"
)
smash.io.read_grid_parameters(
    model_no_dam, path="/home/maxime/DEV/smashbox/smashbox/asset/params"
)

# lecture des données du barrage
rel_hv = pd.read_csv(
    os.path.join("/home/maxime/DATA/Barrage", "loi_H_V.txt"), sep=" "
).to_numpy()
rel_hq = pd.read_csv(
    os.path.join("/home/maxime/DATA/Barrage", "loi_Z_Q.txt"), sep=" "
).to_numpy()
arr_rel_hv = np.expand_dims(rel_hv.transpose(), axis=0)
arr_rel_hq = np.expand_dims(rel_hq.transpose(), axis=0)

# copie des donnée du barrage dans smash
model.set_dam_structure_data("dam_hv", arr_rel_hv)
model.set_dam_structure_data("dam_hq", arr_rel_hq)

model.forward_run()
model_no_dam.forward_run()

final_states = model.rr_final_states.copy()
final_states_no_dam = model_no_dam.rr_final_states.copy()

# simulation d'un evenement
setup["start_time"] = "2015-10-27 00:00"
setup["end_time"] = "2015-10-30 00:00"

model_event = smash.Model(setup, mesh)
model_event_no_dam = smash.Model(setup, mesh_no_dam)

smash.io.read_grid_parameters(
    model_event, path="/home/maxime/DEV/smashbox/smashbox/asset/params"
)
smash.io.read_grid_parameters(
    model_event_no_dam, path="/home/maxime/DEV/smashbox/smashbox/asset/params"
)

model_event.set_dam_structure_data("dam_hv", arr_rel_hv)
model_event.set_dam_structure_data("dam_hq", arr_rel_hq)

# copy des états du modèle
model_event.rr_initial_states = final_states.copy()
model_event_no_dam.rr_initial_states = final_states_no_dam.copy()

model_event.forward_run()
model_event_no_dam.forward_run()

# test en calage
model_event_calib = model_event.copy()

model_event_calib.optimize(
    mapping="distributed",
    optimizer="lbfgsb",
    optimize_options={
        "bounds": {
            "cp": (1, 4000),
            "ct": (1, 2000),
            "kexc": (-10, 5),
            "llr": (1, 2000),
        },
    },
    cost_options={
        "gauge": ["V7155035"],
    },
)

# model sans barage avec param calibré pour pouvoir récupérer les débits entrant dans le barrage 
model_event_no_dam_calib = model_event_no_dam.copy()
model_event_no_dam_calib.rr_parameters = model_event_calib.rr_parameters.copy()
model_event_no_dam_calib.forward_run()

fig, ax = sb.plot.plot.plot_hydrograph(
    model_event,
    plot_rainfall=True,
    plot_settings_sim={"color": "blue", "label": "Simulated with dam"},
)
fig, ax = sb.plot.plot.plot_hydrograph(
    model_event_no_dam,
    figure=(fig, *ax),
    plot_rainfall=False,
    plot_qobs=False,
    plot_settings_sim={"color": "green", "label": "Simulated with no dam"},
)
fig, ax = sb.plot.plot.plot_hydrograph(
    model_event_calib,
    figure=(fig, *ax),
    plot_rainfall=False,
    plot_qobs=False,
    plot_settings_sim={"color": "red", "label": "Calibrated with dam"},
)
fig, ax = sb.plot.plot.plot_hydrograph(
    model_event_no_dam_calib,
    figure=(fig, *ax),
    plot_rainfall=False,
    plot_qobs=False,
    plot_settings_sim={"color": "orange", "label": "Calibrated with no dam"},
)
sb.plot.plot.save_figure(
    fig, figname=os.path.join("/home/maxime/Images", "QSmash_Saint-cecile2.pdf")
)

Voici le graphique obtenu:

QSmash_Saint-cecile2.pdf

et les données du barrage

loi_H_V.txt
loi_Z_Q.txt

maximejay and others added 8 commits February 13, 2026 17:02
Nouveau type dérivé mwd_discontinuities.f90
Modificaion du mesh
Ajout de routine pour le laminage dans md_routing_operator.f90
…ntrée/prélèvement

Status: le code est fonctionnelle
Todo : nettoyage, amélioration de la syntaxe et de l'intégration dans le mesh, différentiation par tapenade
Simplify the mesh
rename hydraulics_discontinuities by hydraulic_structure
Deux ouvrages  sont actuellement codées :
- les barrages
- des injections d'eau

Dans le code fortran:
- Nouveau type dérivé : hydraulic_structureDT
- Nouveaux sous type dérivés : DamDT et InflowDT
Ces type dérivés permettent de stocker des data relative aux fonctionnements des ouvrages hydrauliques dans input_dataDT

Le calcul lié à ces ouvrages est réalisé dans la routine de routage, pour l'instant uniquement lr. Une nouvelle variable d'états permet de stocker la cote d'eau des barrages

Au niveau du maillage:
- chaque ouvrage est défini comme un exutoire
- une nouvelle liste dans generate_mesh permet de définir sir un exutoire est: un exutoire, une jauge, un barrage, une injection d'eau ...

A prévoir:
- ne pas lire les débits si l'exutoire est typé autre que 'jauge'
- standardiser les données relatives aux barrages et aux injection d'eau pour permettre leurs lectures automatiques dans smash
- Comment sortir le débit avant ouvrage si un exutoire est défini, et après l'ouvrage si une structure hydraulique est définie ? utilitée ?
- Gestion des doublons dans le tableau des qsim jauge/structure hydraulique si des jauge et des ouvrages sont identifiés au memes coordonnées ? est-ce problématique ?
- Documentation
- Tests
@maximejay
maximejay requested review from inoelloc and pag13 February 26, 2026 09:43
@maximejay maximejay self-assigned this Feb 26, 2026
@maximejay maximejay added the enhancement New feature or request label Feb 26, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant