From af40f96148ccea7d6107390746cb33b5fb647cf9 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 19 Jun 2026 15:58:13 +0200 Subject: [PATCH 01/11] Add BMI 2.0 C++ interface: DaisyBMI and DaisyPythonController --- CMakeLists.txt | 3 + include/daisy/column.h | 20 + include/daisy/daisy.h | 33 +- include/daisy/lower_boundary/groundwater.h | 1 + include/daisy/soil/soil_water.h | 2 + include/programs/bmi.h | 119 ++++++ include/programs/daisy_bmi.h | 476 +++++++++++++++++++++ src/daisy/column.C | 14 + src/daisy/column_std.C | 208 +++++++++ src/daisy/daisy.C | 103 ++++- src/daisy/soil/soil_water.C | 7 + src/object_model/toplevel.C | 11 +- src/programs/CMakeLists.txt | 65 +++ src/programs/bmi.C | 322 ++++++++++++++ src/programs/bmi_bindings.cpp | 135 ++++++ src/programs/call_python_function.cpp | 33 ++ src/programs/daisy_bmi.C | 361 ++++++++++++++++ 17 files changed, 1908 insertions(+), 5 deletions(-) create mode 100644 include/programs/bmi.h create mode 100644 include/programs/daisy_bmi.h create mode 100644 src/programs/bmi.C create mode 100644 src/programs/bmi_bindings.cpp create mode 100644 src/programs/call_python_function.cpp create mode 100644 src/programs/daisy_bmi.C diff --git a/CMakeLists.txt b/CMakeLists.txt index 0323f6c93..ff7668b4f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -94,6 +94,9 @@ endif() # Sources are added with target_sources in CMakeLists in the source tree add_subdirectory(src) +# Add tools for python interface +add_subdirectory(tools) + # Packaging # lib/ and sample/ contain .dai files that define functionality that should be installed add_subdirectory(lib) diff --git a/include/daisy/column.h b/include/daisy/column.h index 0a458e986..a54b80d51 100644 --- a/include/daisy/column.h +++ b/include/daisy/column.h @@ -25,6 +25,7 @@ #include "object_model/model_framed.h" #include "daisy/manager/irrigate.h" +#include #include class Frame; @@ -136,6 +137,25 @@ class Column : public ModelFramed virtual double soil_inorganic_nitrogen (double from, // [kg N/ha] double to) const = 0; virtual double second_year_utilization () const = 0; + + // Python/BMI coupling: groundwater table and soil state arrays. + virtual double get_groundwater_table () const { return 0.0; } // [cm], neg = below surface + virtual void set_groundwater_table (double) {} // [cm] + virtual double get_bottom_flux () const { return 0.0; } // [cm/h], + = downward + virtual std::vector get_flux_array () const { return {}; } // [cm/h], bottom edge of each cell + virtual std::vector get_h_array () const { return {}; } // [cm], pressure head per layer + virtual std::vector get_theta_array () const { return {}; } // [-], volumetric water content + virtual std::vector get_theta_sat_array () const { return {}; } // [-], saturated water content + virtual double get_runoff_rate () const { return 0.0; } // [mm/day] + virtual double get_column_area () const; // [cm²] + virtual std::vector get_layer_tops () const; // [cm], negative downward + virtual std::vector get_layer_bottoms () const; // [cm], negative downward + // Perturb GW table by dh_cm, re-run Richards, compute Sy = Σ(Δθ·Δz)/dh, + // then restore state. Default: no-op returning 0. + virtual std::tuple, std::vector> + estimate_sy_perturbation (double /*dh_cm*/) + { return {0.0, {}, {}}; } + // Current development stage for the crop named "crop", or // Crop::DSremove if no such crop is present. virtual double crop_ds (symbol crop) const = 0; diff --git a/include/daisy/daisy.h b/include/daisy/daisy.h index 34787ca0a..ef0d86d5b 100644 --- a/include/daisy/daisy.h +++ b/include/daisy/daisy.h @@ -24,6 +24,7 @@ #define DAISY_H #include "programs/program.h" +#include #include #include @@ -37,7 +38,15 @@ class Scope; class Frame; class FrameModel; -class Daisy : public Program +#ifdef __unix +#define DAISY_EXPORT /* nothing */ +#elif defined (BUILD_DLL) +#define DAISY_EXPORT __declspec(dllexport) +#else +#define DAISY_EXPORT __declspec(dllimport) +#endif + +class DAISY_EXPORT Daisy : public Program { public: static const char *const default_description; @@ -57,7 +66,25 @@ class Daisy : public Program void start (); bool is_running () const; void stop (); - + + // Python/BMI coupling: forwarding to first (or pos-th) column. + double get_groundwater_table (unsigned int pos = 0u) const; // [cm] + void set_groundwater_table (double cm, unsigned int pos = 0u); + std::tuple, std::vector> + estimate_sy_perturbation (double dh_cm, unsigned int pos = 0u); + /** Hours from simulation start to the configured stop time, or -1 if open-ended. */ + double stop_duration_hours() const; + + double get_bottom_flux (unsigned int pos = 0u) const; // [cm/h] + std::vector get_flux_array (unsigned int pos = 0u) const; // [cm/h] + std::vector get_h_array (unsigned int pos = 0u) const; // [cm] + std::vector get_theta_array (unsigned int pos = 0u) const; // [-] + std::vector get_theta_sat_array (unsigned int pos = 0u) const; // [-] + double get_runoff_rate (unsigned int pos = 0u) const; // [mm/day] + double get_column_area (unsigned int pos = 0u) const; // [cm²] + std::vector get_layer_tops (unsigned int pos = 0u) const; // [cm] + std::vector get_layer_bottoms (unsigned int pos = 0u) const; // [cm] + // UI. public: void attach_ui (Run* run, const std::vector& logs); @@ -66,6 +93,8 @@ class Daisy : public Program public: bool run (Treelog&); void tick (Treelog&); + void summarize (Treelog&) const; + void close_output (); void output (Log&) const; // Create and Destroy. diff --git a/include/daisy/lower_boundary/groundwater.h b/include/daisy/lower_boundary/groundwater.h index de0471469..47145e2d0 100644 --- a/include/daisy/lower_boundary/groundwater.h +++ b/include/daisy/lower_boundary/groundwater.h @@ -59,6 +59,7 @@ class Groundwater : public ModelDerived // Accessors. public: virtual double table () const = 0; + virtual void set_table (double) {} // [cm], default no-op (override in fixed-table models) // Create and Destroy. public: diff --git a/include/daisy/soil/soil_water.h b/include/daisy/soil/soil_water.h index 65ffd522f..540524e4d 100644 --- a/include/daisy/soil/soil_water.h +++ b/include/daisy/soil/soil_water.h @@ -171,6 +171,8 @@ class SoilWater void set_matrix (const std::vector& h, const std::vector& Theta, const std::vector& q); + // Restore S_sum_ from a saved snapshot (used by estimate_sy_perturbation). + void restore_S_sum (const std::vector& snapshot); void set_tertiary (const std::vector& Theta_p, const std::vector& q_p, const std::vector& S_B2M, diff --git a/include/programs/bmi.h b/include/programs/bmi.h new file mode 100644 index 000000000..b665adc47 --- /dev/null +++ b/include/programs/bmi.h @@ -0,0 +1,119 @@ +// daisy_bmi.h -- BMI (Basic Model Interface) wrapper for Daisy +// +// Implements the BMI 2.0 standard so that Daisy can be driven by +// any BMI-aware framework. +// +// Reference: https://bmi.readthedocs.io/en/stable/ + +#ifndef BMI_H +#define BMI_H + +#include +#include +#include +#include "programs/daisy_bmi.h" + +/** + * @class DaisyBMI + * @brief BMI 2.0 wrapper around DaisyPythonController + * + * Variable names follow a "_" convention. + * + * Input variables (set_value): + * "groundwater__depth" [cm] - groundwater table depth + * "irrigation__rate" [mm/day] + * "land_surface__air_temperature" [degC] (if exposed by Daisy) + * + * Output variables (get_value): + * "soil_water__recharge_rate" [mm/day] + * "land_surface__evapotranspiration_rate" [mm/day] + * "groundwater__depth" [cm] + * "soil_water__transpiration_rate" [mm/day] + * "soil_water__evaporation_rate" [mm/day] + * "soil_water__drainage_rate" [mm/day] + * "soil_water__runoff_rate" [mm/day] + * "vegetation__leaf_area_index" [-] + * "vegetation__root_depth" [cm] + * "vegetation__aboveground_biomass" [g/m2] + * "soil_water__content" [m3/m3] (first layer) + */ +class DaisyBMI +{ +public: + // ===== BMI LIFECYCLE ===== + + /** Initialize model from a .dai config file. */ + void initialize(const std::string& config_file); + + /** Advance model by one timestep (dt = get_time_step()). */ + void update(); + + /** Advance model to a specific time (in days since start). */ + void update_until(double time); + + /** Finalize and release all resources. */ + void finalize(); + + // ===== BMI INFORMATION ===== + + std::string get_component_name() const; + std::vector get_input_var_names() const; + std::vector get_output_var_names() const; + + // ===== BMI TIME ===== + + double get_start_time() const; // always 0.0 + double get_end_time() const; // NaN = open-ended + double get_current_time() const; // days since simulation start + double get_time_step() const; // in days (e.g. 1/24 = hourly) + std::string get_time_units() const; // "d" + + // ===== BMI VARIABLE INFO ===== + + std::string get_var_type(const std::string& name) const; // "double" + std::string get_var_units(const std::string& name) const; + int get_var_itemsize(const std::string& name) const; // sizeof(double) + int get_var_nbytes(const std::string& name) const; // 1 * sizeof(double) + std::string get_var_location(const std::string& name) const; // "none" (scalar) + int get_var_grid(const std::string& name) const; // 0 (single-cell) + + // ===== BMI GRID INFO (single-cell, scalar grid) ===== + + int get_grid_rank(int grid) const; // 0 + int get_grid_size(int grid) const; // 1 + std::string get_grid_type(int grid) const; // "scalar" + + // ===== BMI GET/SET ===== + + /** Get scalar output value by variable name into dest[0]. */ + void get_value(const std::string& name, double* dest) const; + + /** Set scalar input value by variable name from src[0]. */ + void set_value(const std::string& name, const double* src); + + /** + * Estimate specific yield for column `col` via head perturbation + Richards re-solve. + * @param dh_cm Perturbation magnitude in cm (default 1 cm) + * @param col Column index (default 0) + * @return Estimated Sy [-] + */ + std::tuple, std::vector> + estimate_sy_perturbation(double dh_cm = 1.0, unsigned int col = 0u); + + // ===== CONSTRUCTOR / DESTRUCTOR ===== + DaisyBMI(); + ~DaisyBMI(); + +private: + DaisyPythonController ctrl_; + double start_time_days_; // always 0 + double current_time_days_; // updated each update() + double dt_days_; // timestep size in days + + static const std::vector INPUT_VARS; + static const std::vector OUTPUT_VARS; + + double get_output_value(const std::string& name) const; +}; + +#endif // DAISY_BMI_H diff --git a/include/programs/daisy_bmi.h b/include/programs/daisy_bmi.h new file mode 100644 index 000000000..bc6bfe3c2 --- /dev/null +++ b/include/programs/daisy_bmi.h @@ -0,0 +1,476 @@ +// daisy_bmi.h -- Python-controllable Daisy interface +// +// This file is part of Daisy. +// +// Daisy is free software; you can redistribute it and/or modify +// it under the terms of the GNU Lesser Public License as published by +// the Free Software Foundation; either version 2.1 of the License, or +// (at your option) any later version. + +#ifndef DAISY_BMI_H +#define DAISY_BMI_H + +#include +#include +#include +#include + +// Forward declarations +class Toplevel; +class Daisy; + +/** + * @class DaisyPythonController + * @brief Python-friendly interface to control Daisy simulation externally + * + * Allows Python code to: + * - Initialize simulation from .dai config file + * - Advance simulation per timestep + * - Read simulation state (time, water, ET, recharge, etc.) + * - Set boundary conditions (GW depth, rainfall, irrigation, etc.) + * - Control simulation parameters + */ +class DaisyPythonController +{ +private: + // Toplevel owns Metalib, parser, Treelog, and the Daisy program instance. + std::unique_ptr toplevel_; + bool initialized; + bool running; + + // Elapsed simulation time in days, incremented each tick(). + // Matches MODFLOW's approach: a simple counter, independent of calendar dates. + double elapsed_days_ = 0.0; + + // Private helper — returns Daisy& from toplevel (after initialize) + Daisy& daisy() const; + + // Private helper methods + void setup_logging(); + bool load_config_file(const std::string& config_file); + +public: + // ===== INITIALIZATION ===== + + /** + * Initialize simulation from .dai configuration file + * @param config_file Path to .dai configuration file + * @return true if successful, false otherwise + */ + bool initialize(const std::string& config_file); + + /** + * Check if simulation is initialized + * @return true if initialized + */ + bool is_initialized() const { return initialized; } + + + // ===== SIMULATION CONTROL ===== + + /** + * Advance simulation by specified number of days + * @param days Number of days to advance + * @return true if successful, false if error + */ + bool advance(double days); + + /** + * Advance simulation one internal timestep + * @return true if successful, false if error + */ + bool tick(); + + /** + * Start simulation + */ + void start(); + + /** + * Stop simulation + */ + void stop(); + + /** + * Check if simulation is still running + * @return true if simulation is running + */ + bool is_running() const; + + /** + * Finalize simulation and cleanup + * @return true if successful + */ + bool finalize(); + + + // ===== TIME ACCESS ===== + + /** + * Get current simulation time as string + * @return Time formatted as string (YYYY-MM-DD HH:MM) + */ + std::string get_time_string() const; + + /** + * Get current year + * @return Year (e.g., 2020) + */ + int get_year() const; + + /** + * Get current month (1-12) + * @return Month + */ + int get_month() const; + + /** + * Get current day of month (1-31) + * @return Day + */ + int get_day() const; + + /** + * Get current hour (0-23) + * @return Hour + */ + int get_hour() const; + + /** + * Get number of days simulated so far + * @return Days from simulation start + */ + double get_days_since_start() const; + + /** + * Get the total duration of the simulation in days. + * @return Days from simulation start to configured stop time, + * or NaN if the simulation has no fixed end date (open-ended). + */ + double get_stop_days() const; + + /** + * Get current timestep size in hours + * @return Current dt in hours + */ + double get_current_dt_hours() const; + + /** + * Get current timestep size in days + * @return Current dt in days + */ + double get_current_dt_days() const; + //====== DISCRETISATION INFO ====== + double get_column_area () const; + std::vector get_layer_tops () const; + std::vector get_layer_bottoms () const; + + // ===== GROUNDWATER STATE ===== + + /** + * Get groundwater table depth + * @return Groundwater depth in cm below surface (negative = above surface) + */ + double get_groundwater_depth() const; + + /** + * Set groundwater table depth (boundary condition) + * @param depth_cm Depth in cm below surface (negative = above surface) + * @return true if successful + */ + bool set_groundwater_depth(double depth_cm); + + /** + * Estimate specific yield via GW head perturbation + Richards re-solve. + * @param dh_cm Head perturbation in cm (default 1 cm) + * @return Estimated Sy [-] + */ + std::tuple, std::vector> + estimate_sy_perturbation(double dh_cm = 1.0); + + /** + * Get pressure head at specific depth + * @param z_cm Depth in cm (positive downward) + * @return Pressure head in cm at that depth + */ + double get_pressure_head_at_depth(double z_cm) const; + + + // ===== SOIL WATER STATE ===== + + /** + * Get volumetric water content for all soil layers + * @return Vector of theta values (m3/m3) for each layer, top to bottom + */ + std::vector get_soil_water_content() const; + + /** + * Get volumetric water content at specific depth + * @param z_cm Depth in cm (positive downward) + * @return Theta (m3/m3) + */ + double get_soil_water_at_depth(double z_cm) const; + + /** + * Get water potential at specific depth + * @param z_cm Depth in cm (positive downward) + * @return Matric potential in cm + */ + double get_matric_potential_at_depth(double z_cm) const; + + /** + * Get cumulative water content from surface to depth + * @param z_cm Depth in cm (positive downward) + * @return Cumulative water content in mm (0 to depth) + */ + double get_cumulative_water_to_depth(double z_cm) const; + + /** + * Get total plant available water in profile + * @return Available water in mm + */ + double get_available_water() const; + + + // ===== WATER BALANCE ===== + + /** + * Get current actual evapotranspiration + * @return ET rate in mm/day + */ + double get_et_rate() const; + + /** + * Get cumulative evapotranspiration since simulation start + * @return Cumulative ET in mm + */ + double get_cumulative_et() const; + + /** + * Get current actual transpiration from plant + * @return Transpiration in mm/day + */ + double get_transpiration_rate() const; + + /** + * Get current actual evaporation from soil + * @return Evaporation in mm/day + */ + double get_evaporation_rate() const; + + /** + * Get recharge/percolation to groundwater per soil layer + * @return Flux at bottom edge of each layer [mm/day], nlayers long + */ + std::vector get_recharge_rate() const; + + /** + * Get soil pressure head per layer + * @return Pressure head [cm] per layer, nlayers long + */ + std::vector get_pressure_head_array() const; + + /** + * Get volumetric water content per layer + * @return Theta [-] per layer, nlayers long + */ + std::vector get_theta_array() const; + + /** + * Get saturated volumetric water content per layer (static soil parameter) + * @return Theta_sat [-] per layer, nlayers long + */ + std::vector get_theta_sat_array() const; + + /** + * Get cumulative recharge since simulation start + * @return Cumulative recharge in mm + */ + double get_cumulative_recharge() const; + + /** + * Get tile drainage flux + * @return Drainage in mm/day + */ + double get_drainage_rate() const; + + /** + * Get cumulative tile drainage since simulation start + * @return Cumulative drainage in mm + */ + double get_cumulative_drainage() const; + + /** + * Get surface runoff + * @return Runoff in mm/day + */ + double get_runoff_rate() const; + + /** + * Get cumulative runoff since simulation start + * @return Cumulative runoff in mm + */ + double get_cumulative_runoff() const; + + + // ===== CROP/VEGETATION STATE ===== + + /** + * Get leaf area index (LAI) + * @return LAI (m2 leaf / m2 ground) + */ + double get_leaf_area_index() const; + + /** + * Get root depth + * @return Root depth in cm + */ + double get_root_depth() const; + + /** + * Get aboveground biomass + * @return Biomass in g/m2 dry matter + */ + double get_aboveground_biomass() const; + + + // ===== SOIL STATE ===== + + /** + * Get soil temperature at specific depth + * @param z_cm Depth in cm (positive downward) + * @return Temperature in Celsius + */ + double get_soil_temperature_at_depth(double z_cm) const; + + /** + * Get number of soil layers + * @return Number of layers + */ + int get_soil_layer_count() const; + + /** + * Get depth of soil layer + * @param layer Layer index (0-based, top to bottom) + * @return Layer thickness in cm + */ + double get_soil_layer_thickness(int layer) const; + + /** + * Get soil layer top depth + * @param layer Layer index + * @return Depth to top of layer in cm + */ + double get_soil_layer_top(int layer) const; + + /** + * Get soil layer bottom depth + * @param layer Layer index + * @return Depth to bottom of layer in cm + */ + double get_soil_layer_bottom(int layer) const; + + + // ===== WEATHER INPUT ===== + + /** + * Get current rainfall rate + * @return Rainfall in mm/day + */ + double get_rainfall_rate() const; + + /** + * Get cumulative rainfall since simulation start + * @return Cumulative rainfall in mm + */ + double get_cumulative_rainfall() const; + + /** + * Get reference evapotranspiration (PET) + * @return Reference ET in mm/day + */ + double get_reference_et() const; + + /** + * Get air temperature + * @return Temperature in Celsius + */ + double get_air_temperature() const; + + /** + * Get air relative humidity + * @return Relative humidity (0-100) in percent + */ + double get_air_humidity() const; + + /** + * Get wind speed + * @return Wind speed in m/s + */ + double get_wind_speed() const; + + + // ===== NITROGEN STATE ===== + + /** + * Get mineral nitrogen in soil profile + * @return Nitrogen in kg/ha + */ + double get_mineral_nitrogen() const; + + /** + * Get cumulative N leaching + * @return N lost in mm + */ + double get_cumulative_n_leaching() const; + + /** + * Get cumulative N uptake by crop + * @return N uptake in kg/ha + */ + double get_cumulative_n_uptake() const; + + + // ===== DIAGNOSTICS ===== + + /** + * Get last error/warning message + * @return Error message string + */ + std::string get_last_message() const; + + /** + * Check if last operation had errors + * @return true if errors occurred + */ + bool has_errors() const; + + /** + * Get simulation duration for profiling + * @return CPU seconds used + */ + double get_simulation_time() const; + + + // ===== CONSTRUCTOR/DESTRUCTOR ===== + + /** + * Constructor + */ + DaisyPythonController(); + + /** + * Destructor + */ + ~DaisyPythonController(); + + // Delete copy operations + DaisyPythonController(const DaisyPythonController&) = delete; + DaisyPythonController& operator=(const DaisyPythonController&) = delete; + + // Allow move operations + DaisyPythonController(DaisyPythonController&&) noexcept; + DaisyPythonController& operator=(DaisyPythonController&&) noexcept; +}; + +#endif // DAISY_PYTHON_CONTROLLER_H diff --git a/src/daisy/column.C b/src/daisy/column.C index addcf7009..d1a017e3b 100644 --- a/src/daisy/column.C +++ b/src/daisy/column.C @@ -87,4 +87,18 @@ the other processes in Daisy as submodels.") { } } Column_init; +// Default BMI implementations on Column base class. +// Subclasses (e.g. ColumnStandard) override these. +double +Column::get_column_area () const +{ return area; } + +std::vector +Column::get_layer_tops () const +{ return {}; } + +std::vector +Column::get_layer_bottoms () const +{ return {}; } + // column.C ends here. diff --git a/src/daisy/column_std.C b/src/daisy/column_std.C index 2ffd981fd..de150aa74 100644 --- a/src/daisy/column_std.C +++ b/src/daisy/column_std.C @@ -87,6 +87,23 @@ struct ColumnStandard : public Column std::vector tillage_age; std::unique_ptr irrigation; + // Python/BMI coupling: cached arrays (filled in tick_move, valid after each tick). + double bottom_flux_cached_ = 0.0; // [cm/h] + std::vector flux_array_cached_; // [cm/h], bottom edge of each cell + std::vector h_array_cached_; // [cm], pressure head per layer + std::vector theta_array_cached_; // [-], volumetric water content + std::vector theta_sat_cached_; // [-], saturated θ (static soil param) + mutable double runoff_rate_cached_ = 0.0; // [mm], accumulated since last read + + // Snapshots taken just before Richards solve (needed by estimate_sy_perturbation). + Time time_last_; + const Weather* weather_last_ = nullptr; + const Scope* scope_last_ = nullptr; + std::vector h_pretick_; + std::vector theta_pretick_; + double table_pretick_ = 0.0; + std::vector S_sum_pretick_; + // Log variables. double yield_DM; double yield_N; @@ -184,6 +201,21 @@ public: std::string crop_names () const; double bottom () const; + // Python/BMI coupling. + double get_groundwater_table () const override; + void set_groundwater_table (double cm) override; + double get_bottom_flux () const override; + double get_column_area () const override; + std::vector get_layer_tops () const override; + std::vector get_layer_bottoms () const override; + std::vector get_flux_array () const override; + std::vector get_h_array () const override; + std::vector get_theta_array () const override; + std::vector get_theta_sat_array () const override; + double get_runoff_rate () const override; + std::tuple, std::vector> + estimate_sy_perturbation (double dh_cm) override; + // Simulation. void clear (); void tick_source (const Scope&, const Time&, Treelog&); @@ -848,6 +880,22 @@ ColumnStandard::tick_move (const Metalib& metalib, soil_heat->tick (geometry, *soil, *soil_water, T_bottom, *movement, surface->temperature (), dt, msg); soil_water->reset_old (); // Set Theta_old to Theta here. + + // BMI: snapshot h + theta + S_sum + GW table just before Richards solve. + { + const size_t _n = geometry.cell_size (); + h_pretick_.resize (_n); + theta_pretick_.resize (_n); + for (size_t _i = 0; _i < _n; ++_i) + { + h_pretick_[_i] = soil_water->h (_i); + theta_pretick_[_i] = soil_water->Theta (_i); + } + table_pretick_ = groundwater->table (); + S_sum_pretick_.resize (_n); + for (size_t _i = 0; _i < _n; ++_i) + S_sum_pretick_[_i] = soil_water->S_sum (_i); + } chemistry->mass_balance (geometry, *soil_water); soil_water->tick_ice (geometry, *soil, dt, msg); movement->tick (*soil, *soil_water, *soil_heat, @@ -891,6 +939,27 @@ ColumnStandard::tick_move (const Metalib& metalib, } chemistry->update_C (*soil, *soil_water, *soil_heat, *awi); chemistry->mass_balance (geometry, *soil_water); + + // BMI: cache post-Richards arrays and runoff for BMI getters. + { + const size_t n = geometry.cell_size (); + time_last_ = time_end; + weather_last_ = weather.get () ? weather.get () : global_weather; + scope_last_ = extern_scope; + + bottom_flux_cached_ = soil_water->q_matrix (n); + + flux_array_cached_.resize (n); + h_array_cached_.resize (n); + theta_array_cached_.resize (n); + for (size_t i = 0; i < n; ++i) + { + flux_array_cached_[i] = soil_water->q_matrix (i + 1); // [cm/h] + h_array_cached_[i] = soil_water->h (i); // [cm] + theta_array_cached_[i] = soil_water->Theta (i); // [-] + } + runoff_rate_cached_ += surface->runoff_rate () * surface->ponding_average () * dt; // [mm] + } } bool @@ -1211,6 +1280,14 @@ ColumnStandard::initialize (const Metalib& metalib, // Soil conductivity and capacity logs. soil_heat->tick_after (geometry.cell_size (), *soil, *soil_water, msg); + // BMI: cache static soil parameter theta_sat once after initialization. + { + const size_t n = geometry.cell_size (); + theta_sat_cached_.resize (n); + for (size_t i = 0; i < n; ++i) + theta_sat_cached_[i] = soil->Theta_sat (i); + } + // Litter layer. litter->tick (*bioclimate, geometry, *soil, *soil_water, *soil_heat, *organic_matter, *chemistry, 0.0, msg); @@ -1363,3 +1440,134 @@ Hansen et.al. 1990. with generic movement in soil.") } column_syntax; // column_std.C ends here. + +// ===== Python/BMI getter implementations ===== + +double +ColumnStandard::get_groundwater_table () const +{ return groundwater->table (); } + +void +ColumnStandard::set_groundwater_table (double cm) +{ groundwater->set_table (cm); } + +double +ColumnStandard::get_bottom_flux () const +{ return bottom_flux_cached_; } + +double +ColumnStandard::get_column_area () const +{ return area; } + +std::vector +ColumnStandard::get_layer_tops () const +{ + const size_t n = geometry.cell_size (); + std::vector tops (n); + for (size_t i = 0; i < n; ++i) + tops[i] = geometry.cell_top (i); + return tops; +} + +std::vector +ColumnStandard::get_layer_bottoms () const +{ + const size_t n = geometry.cell_size (); + std::vector bottoms (n); + for (size_t i = 0; i < n; ++i) + bottoms[i] = geometry.cell_bottom (i); + return bottoms; +} + +std::vector +ColumnStandard::get_flux_array () const +{ return flux_array_cached_; } + +std::vector +ColumnStandard::get_h_array () const +{ return h_array_cached_; } + +std::vector +ColumnStandard::get_theta_array () const +{ return theta_array_cached_; } + +std::vector +ColumnStandard::get_theta_sat_array () const +{ return theta_sat_cached_; } + +double +ColumnStandard::get_runoff_rate () const +{ + const double v = runoff_rate_cached_; + runoff_rate_cached_ = 0.0; + return v; +} + +auto ColumnStandard::estimate_sy_perturbation (double dh_cm) + -> std::tuple, std::vector> +{ + const std::tuple, std::vector> failed + = {0.0, {}, {}}; + + if (!weather_last_ || !scope_last_ || h_pretick_.empty () || S_sum_pretick_.empty ()) + return failed; + + const double dh_safe = std::min (dh_cm, -table_pretick_); + if (dh_safe <= 0.0) + return failed; + + const size_t n = geometry.cell_size (); + + auto restore = [&]() + { + for (size_t i = 0; i < n; ++i) + soil_water->set_content (i, h_pretick_[i], theta_pretick_[i]); + soil_water->restore_S_sum (S_sum_pretick_); + }; + + struct TickResult { std::vector theta, flux_mm_d, h_cm; }; + auto run_tick = [&](double gw_table) -> TickResult + { + groundwater->set_table (gw_table); + soil_water->reset_old (); + try + { + movement->tick (*soil, *soil_water, *soil_heat, + *surface, *groundwater, + time_last_, *scope_last_, *weather_last_, + 24.0, Treelog::null ()); + } + catch (...) {} + TickResult r; + r.theta.resize (n); r.flux_mm_d.resize (n); r.h_cm.resize (n); + for (size_t i = 0; i < n; ++i) + { + r.theta[i] = soil_water->Theta (i); + r.h_cm[i] = soil_water->h (i); + r.flux_mm_d[i] = soil_water->q_matrix (i + 1) * 10.0 * 24.0; + } + return r; + }; + + restore (); + const TickResult A = run_tick (table_pretick_); + + restore (); + const TickResult B = run_tick (table_pretick_ + dh_safe); + + double numerator = 0.0; + for (size_t i = 0; i < n; ++i) + { + const double dtheta = B.theta[i] - A.theta[i]; + const double dz_m = (geometry.cell_top (i) - geometry.cell_bottom (i)) * 1e-2; + numerator += dtheta * dz_m; + } + const double sy = std::abs (numerator) / std::abs (dh_safe * 1e-2); + + groundwater->set_table (table_pretick_); + for (size_t i = 0; i < n; ++i) + soil_water->set_content (i, h_array_cached_[i], theta_array_cached_[i]); + soil_water->restore_S_sum (std::vector (n, 0.0)); + + return {sy, B.flux_mm_d, B.h_cm}; +} diff --git a/src/daisy/daisy.C b/src/daisy/daisy.C index 0fd0603a9..14a046c56 100644 --- a/src/daisy/daisy.C +++ b/src/daisy/daisy.C @@ -61,7 +61,7 @@ public: const std::unique_ptr scopesel; const Scope* extern_scope; const std::unique_ptr print_time; - const std::unique_ptr output_log; + std::unique_ptr output_log; // non-const: close_output() resets it const bool message_timestep; const Timestep timestep; const double max_dt; @@ -568,6 +568,107 @@ the simulation. Can be overwritten by column specific weather."); Daisy::~Daisy () { } +// ===== Python/BMI forwarding methods ===== + +double +Daisy::get_groundwater_table (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_groundwater_table () : 0.0; +} + +void +Daisy::set_groundwater_table (double cm, unsigned int pos) +{ + Column* col = impl->field->find (pos); + if (col) col->set_groundwater_table (cm); +} + +auto Daisy::estimate_sy_perturbation (double dh_cm, unsigned int pos) + -> std::tuple, std::vector> +{ + Column* col = impl->field->find (pos); + return col ? col->estimate_sy_perturbation (dh_cm) + : std::make_tuple (0.0, std::vector{}, std::vector{}); +} + +double +Daisy::stop_duration_hours() const +{ + return impl->duration; +} + +double +Daisy::get_bottom_flux (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_bottom_flux () : 0.0; +} + +std::vector +Daisy::get_flux_array (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_flux_array () : std::vector{}; +} + +std::vector +Daisy::get_h_array (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_h_array () : std::vector{}; +} + +std::vector +Daisy::get_theta_array (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_theta_array () : std::vector{}; +} + +std::vector +Daisy::get_theta_sat_array (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_theta_sat_array () : std::vector{}; +} + +double +Daisy::get_runoff_rate (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_runoff_rate () : 0.0; +} + +double +Daisy::get_column_area (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_column_area () : 0.0; +} + +std::vector +Daisy::get_layer_tops (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_layer_tops () : std::vector{}; +} + +std::vector +Daisy::get_layer_bottoms (unsigned int pos) const +{ + Column* col = impl->field->find (pos); + return col ? col->get_layer_bottoms () : std::vector{}; +} + +void +Daisy::summarize (Treelog& msg) const +{ impl->summarize (msg); } + +void +Daisy::close_output () +{ impl->output_log.reset (); } + static struct ProgramDaisySyntax : public DeclareModel { Model* make (const BlockModel& al) const diff --git a/src/daisy/soil/soil_water.C b/src/daisy/soil/soil_water.C index ad0f6c5d0..64736bd2b 100644 --- a/src/daisy/soil/soil_water.C +++ b/src/daisy/soil/soil_water.C @@ -239,6 +239,13 @@ SoilWater::set_matrix (const std::vector& h, q_matrix_ = q; } +void +SoilWater::restore_S_sum (const std::vector& snapshot) +{ + daisy_assert (S_sum_.size () == snapshot.size ()); + S_sum_ = snapshot; +} + void SoilWater::set_tertiary (const std::vector& Theta_p, const std::vector& q_p, diff --git a/src/object_model/toplevel.C b/src/object_model/toplevel.C index 666f1aeaa..d7ba24e81 100644 --- a/src/object_model/toplevel.C +++ b/src/object_model/toplevel.C @@ -46,13 +46,16 @@ #ifdef BUILD_PYTHON #include +#include #endif struct Toplevel::Implementation : boost::noncopyable { #ifdef BUILD_PYTHON - // Start python interpreter - pybind11::scoped_interpreter guard; + // When running as a standalone executable, we own the Python interpreter. + // When loaded as a Python extension (BMI), Python is already running — + // skip initialisation to avoid "interpreter is already running" error. + std::optional guard; #endif const symbol preferred_ui; const std::string program_name; @@ -195,6 +198,10 @@ Toplevel::Implementation::Implementation (Metalib::load_frame_t load_syntax, has_daisy_log (false) { (void) setlocale (LC_ALL, "C"); +#ifdef BUILD_PYTHON + if (!Py_IsInitialized ()) + guard.emplace (); +#endif } Toplevel::Implementation::~Implementation () diff --git a/src/programs/CMakeLists.txt b/src/programs/CMakeLists.txt index 37361a525..76d0934e2 100644 --- a/src/programs/CMakeLists.txt +++ b/src/programs/CMakeLists.txt @@ -18,3 +18,68 @@ target_sources(${DAISY_CORE_NAME} PRIVATE program_spawn.C program_weather.C ) + +# ===== daisy_bmi pybind11 module ===== +pybind11_add_module(daisy_bmi bmi_bindings.cpp bmi.C daisy_bmi.C) +target_link_libraries(daisy_bmi PRIVATE ${DAISY_CORE_NAME}) +target_include_directories(daisy_bmi PRIVATE + ${CMAKE_SOURCE_DIR}/include +) + + +# ===== Copy Python module and runtime dependencies ===== +if (WIN32) + # Where the Python package lives (src layout) + set(PYTHON_PACKAGE_DIR "${CMAKE_SOURCE_DIR}/python/src/daisy") + + # Derive the MinGW/MSYS2 bin directory from the compiler path. + # This avoids any hardcoded path: wherever CMake found the compiler is + # where the runtime DLLs live (e.g. C:/msys64/ucrt64/bin). + get_filename_component(COMPILER_DIR "${CMAKE_CXX_COMPILER}" DIRECTORY) + + message(STATUS "daisy_bmi: Python package dir : ${PYTHON_PACKAGE_DIR}") + message(STATUS "daisy_bmi: MinGW runtime dir : ${COMPILER_DIR}") + + add_custom_command(TARGET daisy_bmi POST_BUILD + + # Ensure target directory exists + COMMAND ${CMAKE_COMMAND} -E make_directory "${PYTHON_PACKAGE_DIR}" + + # --- Copy the Python module (.pyd) --- + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "$" + "${PYTHON_PACKAGE_DIR}" + + # --- Copy Daisy core DLL --- + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "$" + "${PYTHON_PACKAGE_DIR}" + + # --- Copy MinGW/MSYS2 runtime DLLs --- + # Direct runtime dependencies of libcore.dll: + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libstdc++-6.dll" + "${PYTHON_PACKAGE_DIR}" + + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libgcc_s_seh-1.dll" + "${PYTHON_PACKAGE_DIR}" + + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libwinpthread-1.dll" + "${PYTHON_PACKAGE_DIR}" + + # Transitive dependencies (libcore -> libcxsparse -> libsuitesparseconfig -> libgomp-1): + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libcxsparse.dll" + "${PYTHON_PACKAGE_DIR}" + + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libsuitesparseconfig.dll" + "${PYTHON_PACKAGE_DIR}" + + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "${COMPILER_DIR}/libgomp-1.dll" + "${PYTHON_PACKAGE_DIR}" + ) +endif() \ No newline at end of file diff --git a/src/programs/bmi.C b/src/programs/bmi.C new file mode 100644 index 000000000..8e97cb929 --- /dev/null +++ b/src/programs/bmi.C @@ -0,0 +1,322 @@ +// daisy_bmi.C -- BMI 2.0 wrapper implementation for Daisy +// +// This file is part of Daisy. + +#include "programs/bmi.h" +#include +#include +#include + +// ===== STATIC VARIABLE LISTS ===== + +const std::vector DaisyBMI::INPUT_VARS = { + "groundwater__depth", // [cm] below surface + "irrigation__rate", // [mm/day] +}; + +const std::vector DaisyBMI::OUTPUT_VARS = { + "soil_water__recharge_rate", // [mm/day] + "land_surface__evapotranspiration_rate", // [mm/day] + "soil_water__transpiration_rate", // [mm/day] + "soil_water__evaporation_rate", // [mm/day] + "soil_water__drainage_rate", // [mm/day] + "soil_water__runoff_rate", // [mm/day] + "groundwater__depth", // [cm] + "vegetation__leaf_area_index", // [-] + "vegetation__root_depth", // [cm] + "vegetation__aboveground_biomass", // [g/m2] + "soil_water__content", // [m3/m3] top layer + "column__area", // [cm2] scalar + "soil_layer__top_depth", // [cm] array, N layers + "soil_layer__bottom_depth", // [cm] array, N layers + "soil_water__pressure_head", // [cm] array, N layers + "soil_water__theta", // [-] array, N layers + "soil__theta_sat", // [-] array, N layers (static soil param) +}; + +// ===== CONSTRUCTOR / DESTRUCTOR ===== + +DaisyBMI::DaisyBMI() + : start_time_days_(0.0) + , current_time_days_(0.0) + , dt_days_(1.0 / 24.0) // default: hourly timestep +{ +} + +DaisyBMI::~DaisyBMI() +{ + // DaisyPythonController destructor handles cleanup +} + +// ===== LIFECYCLE ===== + +void DaisyBMI::initialize(const std::string& config_file) +{ + try + { + ctrl_.initialize(config_file); + } + catch (const std::exception& e) + { + throw std::runtime_error( + std::string("DaisyBMI::initialize failed for config: ") + config_file + + "\n Reason: " + e.what()); + } + catch (const int code) + { + throw std::runtime_error( + std::string("DaisyBMI::initialize failed for config: ") + config_file + + "\n Reason: Daisy reported an error (check daisy.log for details). Exit code: " + + std::to_string(code)); + } + catch (const char* msg) + { + throw std::runtime_error( + std::string("DaisyBMI::initialize failed for config: ") + config_file + + "\n Reason: " + msg); + } + catch (...) + { + throw std::runtime_error( + std::string("DaisyBMI::initialize failed for config: ") + config_file + + "\n Reason: unknown exception"); + } + + if (!ctrl_.is_initialized()) + throw std::runtime_error( + std::string("DaisyBMI::initialize failed for config: ") + config_file); + + // Use Daisy's actual internal time as the reference point. + // get_days_since_start() returns 0 right after init, but reading it through + // the controller ensures we are in sync with Daisy's clock from the start. + start_time_days_ = ctrl_.get_days_since_start(); // = 0.0 at init + current_time_days_ = start_time_days_; + + dt_days_ = ctrl_.get_current_dt_days(); + if (dt_days_ <= 0.0) + dt_days_ = 1.0 / 24.0; // fallback: hourly +} + +void DaisyBMI::update() +{ + if (!ctrl_.tick()) + throw std::runtime_error("DaisyBMI::update (tick) failed"); + + // Sync current time directly from Daisy's internal clock. + current_time_days_ = ctrl_.get_days_since_start(); + + // Keep dt in sync with Daisy's actual dt. + double actual_dt = ctrl_.get_current_dt_days(); + if (actual_dt > 0.0) + dt_days_ = actual_dt; +} + +void DaisyBMI::update_until(double time) +{ + // Use half-timestep epsilon to avoid running one tick too many due to + // floating-point accumulation (e.g. 24 * 1/24 is not exactly 1.0). + const double eps = dt_days_ * 0.5; + while (current_time_days_ < time - eps && ctrl_.is_running()) + update(); + + // Snap to the exact requested time so callers always see a clean value + // and drift does not accumulate across successive update_until calls. + if (ctrl_.is_running()) + current_time_days_ = time; +} + +void DaisyBMI::finalize() +{ + ctrl_.finalize(); +} + +// ===== INFORMATION ===== + +std::string DaisyBMI::get_component_name() const +{ + return "Daisy"; +} + +std::vector DaisyBMI::get_input_var_names() const +{ + return INPUT_VARS; +} + +std::vector DaisyBMI::get_output_var_names() const +{ + return OUTPUT_VARS; +} + +// ===== TIME ===== + +double DaisyBMI::get_start_time() const +{ + return start_time_days_; +} + +double DaisyBMI::get_end_time() const +{ + return ctrl_.get_stop_days(); +} + +double DaisyBMI::get_current_time() const +{ + return current_time_days_; +} + +double DaisyBMI::get_time_step() const +{ + return dt_days_; +} + +std::string DaisyBMI::get_time_units() const +{ + return "d"; +} + +// ===== VARIABLE INFO ===== + +std::string DaisyBMI::get_var_type(const std::string& /*name*/) const +{ + return "double"; +} + +std::string DaisyBMI::get_var_units(const std::string& name) const +{ + if (name == "groundwater__depth") return "cm"; + if (name == "irrigation__rate") return "mm d-1"; + if (name == "soil_water__recharge_rate") return "mm d-1"; + if (name == "land_surface__evapotranspiration_rate") return "mm d-1"; + if (name == "soil_water__transpiration_rate") return "mm d-1"; + if (name == "soil_water__evaporation_rate") return "mm d-1"; + if (name == "soil_water__drainage_rate") return "mm d-1"; + if (name == "soil_water__runoff_rate") return "mm d-1"; + if (name == "vegetation__leaf_area_index") return "1"; + if (name == "vegetation__root_depth") return "cm"; + if (name == "vegetation__aboveground_biomass") return "g m-2"; + if (name == "soil_water__content") return "m3 m-3"; + if (name == "column__area") return "cm2"; + if (name == "soil_layer__top_depth") return "cm"; + if (name == "soil_layer__bottom_depth") return "cm"; + if (name == "soil_water__pressure_head") return "cm"; + if (name == "soil_water__theta") return "1"; + if (name == "soil__theta_sat") return "1"; + throw std::invalid_argument("Unknown variable: " + name); +} + +int DaisyBMI::get_var_nbytes(const std::string& name) const { + if (name == "soil_layer__top_depth" || name == "soil_layer__bottom_depth" + || name == "soil_water__recharge_rate" + || name == "soil_water__pressure_head" + || name == "soil_water__theta" + || name == "soil__theta_sat") + return ctrl_.get_soil_layer_count() * sizeof(double); + return sizeof(double); +} + +int DaisyBMI::get_var_itemsize(const std::string& /*name*/) const +{ + return static_cast(sizeof(double)); +} + +std::string DaisyBMI::get_var_location(const std::string& /*name*/) const +{ + return "none"; // scalars have no grid location +} + +int DaisyBMI::get_var_grid(const std::string& /*name*/) const +{ + return 0; // all vars on the single-cell scalar grid +} + +// ===== GRID INFO ===== + +int DaisyBMI::get_grid_rank(int /*grid*/) const { return 0; } +int DaisyBMI::get_grid_size(int /*grid*/) const { return 1; } +std::string DaisyBMI::get_grid_type(int /*grid*/) const { return "scalar"; } + +// ===== GET / SET ===== + +double DaisyBMI::get_output_value(const std::string& name) const +{ + if (name == "soil_water__recharge_rate") { + auto v = ctrl_.get_recharge_rate(); + return v.empty() ? 0.0 : v.back(); // scalar fallback = bottom-layer value + } + if (name == "land_surface__evapotranspiration_rate") return ctrl_.get_et_rate(); + if (name == "soil_water__transpiration_rate") return ctrl_.get_transpiration_rate(); + if (name == "soil_water__evaporation_rate") return ctrl_.get_evaporation_rate(); + if (name == "soil_water__drainage_rate") return ctrl_.get_drainage_rate(); + if (name == "soil_water__runoff_rate") return ctrl_.get_runoff_rate(); + if (name == "groundwater__depth") return ctrl_.get_groundwater_depth(); + if (name == "vegetation__leaf_area_index") return ctrl_.get_leaf_area_index(); + if (name == "vegetation__root_depth") return ctrl_.get_root_depth(); + if (name == "vegetation__aboveground_biomass") return ctrl_.get_aboveground_biomass(); + if (name == "soil_water__content") + { + auto layers = ctrl_.get_soil_water_content(); + return layers.empty() ? std::numeric_limits::quiet_NaN() : layers[0]; + } + if (name == "column__area") return ctrl_.get_column_area(); + throw std::invalid_argument("Unknown output variable: " + name); +} + +void DaisyBMI::get_value(const std::string& name, double* dest) const +{ + if (name == "soil_layer__top_depth") { + auto v = ctrl_.get_layer_tops(); + std::copy(v.begin(), v.end(), dest); + return; + } + if (name == "soil_layer__bottom_depth") { + auto v = ctrl_.get_layer_bottoms(); + std::copy(v.begin(), v.end(), dest); + return; + } + if (name == "soil_water__recharge_rate") { + auto v = ctrl_.get_recharge_rate(); + std::copy(v.begin(), v.end(), dest); + return; + } + if (name == "soil_water__pressure_head") { + auto v = ctrl_.get_pressure_head_array(); + std::copy(v.begin(), v.end(), dest); + return; + } + if (name == "soil_water__theta") { + auto v = ctrl_.get_theta_array(); + std::copy(v.begin(), v.end(), dest); + return; + } + if (name == "soil__theta_sat") { + auto v = ctrl_.get_theta_sat_array(); + std::copy(v.begin(), v.end(), dest); + return; + } + dest[0] = get_output_value(name); // bestaand pad voor scalars +} + +void DaisyBMI::set_value(const std::string& name, const double* src) +{ + if (name == "groundwater__depth") + { + ctrl_.set_groundwater_depth(src[0]); + return; + } + if (name == "irrigation__rate") + { + // TODO: expose irrigation setter in DaisyPythonController + // ctrl_.set_irrigation_rate(src[0]); + throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyPythonController"); + } + throw std::invalid_argument("Unknown input variable: " + name); +} + +auto DaisyBMI::estimate_sy_perturbation(double dh_cm, unsigned int col) + -> std::tuple, std::vector> +{ + if (col != 0u) + throw std::invalid_argument("estimate_sy_perturbation: col > 0 not yet supported via DaisyBMI"); + return ctrl_.estimate_sy_perturbation(dh_cm); +} + diff --git a/src/programs/bmi_bindings.cpp b/src/programs/bmi_bindings.cpp new file mode 100644 index 000000000..94001f0ab --- /dev/null +++ b/src/programs/bmi_bindings.cpp @@ -0,0 +1,135 @@ +// bmi_bindings.cpp -- pybind11 bindings for DaisyBMI +// +// Exposes DaisyBMI to Python as the module "daisy_bmi". +// +// Compile (example, adjust paths): +// g++ -O3 -Wall -shared -std=c++17 -fPIC +// $(python3 -m pybind11 --includes) +// daisy_bmi_bindings.cpp daisy_bmi.C daisy_python_controller.C +// -o daisy_bmi$(python3-config --extension-suffix) +// -L. -ldaisy_core +// +// Usage from Python: +// import daisy_bmi +// m = daisy_bmi.DaisyBMI() +// m.initialize("myconfig.dai") +// while m.get_current_time() < 365.0: +// recharge = m.get_value("soil_water__recharge_rate") +// m.set_value("groundwater__depth", 150.0) +// m.update() +// m.finalize() + +#include +#include +#include +#include "programs/bmi.h" + +namespace py = pybind11; + +PYBIND11_MODULE(daisy_bmi, m) +{ + m.doc() = R"pbdoc( + daisy_bmi - BMI 2.0 interface to the Daisy soil-crop-atmosphere model + ======================================================================= + + Example + ------- + import daisy_bmi + + bmi = daisy_bmi.DaisyBMI() + bmi.initialize("myconfig.dai") + + print("Component:", bmi.get_component_name()) + print("Time step:", bmi.get_time_step(), bmi.get_time_units()) + print("Inputs:", bmi.get_input_var_names()) + print("Outputs:", bmi.get_output_var_names()) + + while bmi.get_current_time() < 365.0: + # Pull outputs + recharge = bmi.get_value("soil_water__recharge_rate") + et = bmi.get_value("land_surface__evapotranspiration_rate") + + # Push inputs + bmi.set_value("groundwater__depth", 150.0) # 150 cm below surface + + bmi.update() + + bmi.finalize() + )pbdoc"; + + py::class_(m, "DaisyBMI") + + .def(py::init<>(), "Create a new Daisy BMI instance") + + // --- Lifecycle --- + .def("initialize", &DaisyBMI::initialize, py::arg("config_file"), + "Initialize from a .dai config file") + .def("update", &DaisyBMI::update, + "Advance model by one timestep") + .def("update_until", &DaisyBMI::update_until, py::arg("time"), + "Advance model until time (days since start)") + .def("finalize", &DaisyBMI::finalize, + "Finalize and release all resources") + + // --- Information --- + .def("get_component_name", &DaisyBMI::get_component_name) + .def("get_input_var_names", &DaisyBMI::get_input_var_names) + .def("get_output_var_names", &DaisyBMI::get_output_var_names) + + // --- Time --- + .def("get_start_time", &DaisyBMI::get_start_time, + "Simulation start time [days]") + .def("get_end_time", &DaisyBMI::get_end_time, + "Simulation end time [days] (NaN = open-ended)") + .def("get_current_time", &DaisyBMI::get_current_time, + "Current simulation time [days since start]") + .def("get_time_step", &DaisyBMI::get_time_step, + "Current timestep size [days]") + .def("get_time_units", &DaisyBMI::get_time_units, + "Time unit string, e.g. 'd'") + + // --- Variable metadata --- + .def("get_var_type", &DaisyBMI::get_var_type, py::arg("name")) + .def("get_var_units", &DaisyBMI::get_var_units, py::arg("name")) + .def("get_var_itemsize", &DaisyBMI::get_var_itemsize, py::arg("name")) + .def("get_var_nbytes", &DaisyBMI::get_var_nbytes, py::arg("name")) + .def("get_var_location", &DaisyBMI::get_var_location, py::arg("name")) + .def("get_var_grid", &DaisyBMI::get_var_grid, py::arg("name")) + + // --- Grid metadata --- + .def("get_grid_rank", &DaisyBMI::get_grid_rank, py::arg("grid")) + .def("get_grid_size", &DaisyBMI::get_grid_size, py::arg("grid")) + .def("get_grid_type", &DaisyBMI::get_grid_type, py::arg("grid")) + + // --- Get / Set (scalar convenience wrappers) --- + .def("get_value", + [](const DaisyBMI& self, const std::string& name) -> double { + // Allocate the correct-sized buffer; array vars (e.g. + // soil_water__recharge_rate) hold one value per soil layer. + // Writing N elements into a 1-element stack variable is UB and + // caused WinError 10054 crashes via stack corruption. + int n = self.get_var_nbytes(name) / static_cast(sizeof(double)); + std::vector buf(n); + self.get_value(name, buf.data()); + return buf.back(); // scalar = bottom-layer value, matches get_output_value + }, + py::arg("name"), + "Get a scalar output value by BMI variable name") + + .def("set_value", + [](DaisyBMI& self, const std::string& name, double value) { + self.set_value(name, &value); + }, + py::arg("name"), py::arg("value"), + "Set a scalar input value by BMI variable name") + + .def("get_value_array", + [](const DaisyBMI& self, const std::string& name) { + int n = self.get_var_nbytes(name) / static_cast(sizeof(double)); + py::array_t arr(n); + self.get_value(name, arr.mutable_data()); + return arr; + }, + py::arg("name"), + "Get an array output value as a numpy array of doubles"); +} diff --git a/src/programs/call_python_function.cpp b/src/programs/call_python_function.cpp new file mode 100644 index 000000000..2db1eaf73 --- /dev/null +++ b/src/programs/call_python_function.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +namespace py = pybind11; +using namespace py::literals; + +int main(int argc, char* argv[]) { + if (argc < 3) { + std::cout << "Usage: call_python_function " \ + "[]*" << std::endl; + return 1; + } + py::scoped_interpreter guard{}; + + auto py_module = py::module_::import(argv[1]); + + std::vector parameters; + for (std::size_t i = 3; i < argc; i++) { + parameters.push_back(argv[i]); + } + py::object py_object = py_module.attr(argv[2])(parameters); + + std::vector inputs {1, 10, 20, 30, 40}; + std::string domain = py_object.attr("domain").cast(); + std::string range = py_object.attr("range").cast(); + for (auto input : inputs) { + double py_result = py_object(input).cast(); + std::cout << "Input: " << input << " " << domain << std::endl + << "Result: " << py_result << " " << range << std::endl; + } + return 0; +} diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C new file mode 100644 index 000000000..93a50f438 --- /dev/null +++ b/src/programs/daisy_bmi.C @@ -0,0 +1,361 @@ +// daisy_bmi.C -- Python-controllable Daisy interface implementation +// +// This file is part of Daisy. + +#define BUILD_DLL + +#include "programs/daisy_bmi.h" +#include "object_model/toplevel.h" +#include "object_model/treelog.h" +#include "daisy/daisy.h" +#include "daisy/daisy_time.h" +#include "daisy/field.h" +#include +#include +#include +#include +#include +#include + +// ===== PRIVATE HELPERS ===== + +Daisy& DaisyPythonController::daisy() const +{ + return dynamic_cast(toplevel_->program()); +} + +void DaisyPythonController::setup_logging() +{ + if (toplevel_) toplevel_->set_ui_none(); +} + +// Set DAISYHOME to the parent of the directory containing libcore.dll +// (i.e. the daisy source/install root), but only if DAISYHOME is not +// already set by the user. This ensures that when Daisy runs as a Python +// extension, it finds lib/ and sample/ relative to itself — not relative +// to python.exe. +static void ensure_daisy_home() +{ + if (getenv("DAISYHOME")) + return; // user override takes priority + + try + { + // boost::dll::this_line_location() returns the path of the DLL that + // contains the calling code — i.e. libcore.dll (or daisy_bmi.pyd). + std::filesystem::path dll_path = boost::dll::this_line_location().native(); + // The DLLs sit in /python/src/daisy/ (installed) or + // / (build tree). Walk up to find a directory that contains + // both lib/ and sample/ — that is the Daisy home. + std::filesystem::path candidate = dll_path.parent_path(); + for (int i = 0; i < 5; ++i) + { + if (std::filesystem::is_directory(candidate / "lib") && + std::filesystem::is_directory(candidate / "sample")) + { +#ifdef _WIN32 + _putenv_s("DAISYHOME", candidate.string().c_str()); +#else + setenv("DAISYHOME", candidate.string().c_str(), 0); +#endif + return; + } + candidate = candidate.parent_path(); + } + } + catch (...) {} // non-fatal: Daisy will fall back to its own detection +} + +bool DaisyPythonController::load_config_file(const std::string& config_file) +{ + ensure_daisy_home(); + toplevel_ = std::make_unique("daisy"); + toplevel_->set_ui_none(); + + std::filesystem::path cfg_path = std::filesystem::absolute(config_file); + std::filesystem::path cfg_dir = cfg_path.parent_path(); + std::filesystem::path old_dir = std::filesystem::current_path(); + try + { + std::filesystem::current_path(cfg_dir); + toplevel_->parse_file(cfg_path.filename().string().c_str()); + std::filesystem::current_path(old_dir); + } + catch (...) + { + std::filesystem::current_path(old_dir); + throw; + } + return true; +} + +// ===== INITIALIZATION ===== + +bool DaisyPythonController::initialize(const std::string& config_file) +{ + // Let exceptions propagate — the caller (DaisyBMI) will wrap them + // with a meaningful message. A silent 'return false' hides the root cause. + load_config_file(config_file); + start(); + return initialized; +} + +// ===== LIFECYCLE ===== + +void DaisyPythonController::start() +{ + if (!toplevel_) return; + toplevel_->initialize(); // moves state from is_uninitialized -> is_ready + toplevel_->msg().open("Running"); + try + { + daisy().start(); + initialized = true; + running = daisy().is_running(); + } + catch (...) + { + toplevel_->msg().close(); + throw; + } +} + +bool DaisyPythonController::advance(double days) +{ + if (!initialized) return false; + try + { + int steps = static_cast(days * 24); + for (int i = 0; i < steps && running; ++i) + { + daisy().tick(toplevel_->msg()); + running = daisy().is_running(); + } + elapsed_days_ += days; + return running; + } + catch (...) { return false; } +} + +bool DaisyPythonController::tick() +{ + if (!initialized) return false; + try + { + daisy().tick(toplevel_->msg()); + running = daisy().is_running(); + if (running) elapsed_days_ += 1.0 / 24.0; + return running; + } + catch (...) { return false; } +} + +void DaisyPythonController::stop() +{ + running = false; +} + +bool DaisyPythonController::is_running() const +{ + return running; +} + +bool DaisyPythonController::finalize() +{ + if (!initialized) return true; + daisy().stop(); + daisy().summarize (toplevel_->msg()); + daisy().close_output (); // flush + close all DLF file handles + toplevel_->msg().close(); + initialized = false; + running = false; + return true; +} + +// ===== TIME ===== + +double DaisyPythonController::get_days_since_start() const { return elapsed_days_; } + +double DaisyPythonController::get_stop_days() const +{ + if (!initialized) return std::numeric_limits::quiet_NaN(); + const double h = daisy().stop_duration_hours(); + return (h >= 0.0) ? h / 24.0 : std::numeric_limits::quiet_NaN(); +} + +std::string DaisyPythonController::get_time_string() const +{ + if (!initialized) return ""; + const Time& t = daisy().time(); + char buf[32]; + std::snprintf(buf, sizeof(buf), "%04d-%02d-%02d %02d:00", + t.year(), t.month(), t.mday(), t.hour()); + return buf; +} + +int DaisyPythonController::get_year() const { return initialized ? daisy().time().year() : 0; } +int DaisyPythonController::get_month() const { return initialized ? daisy().time().month() : 0; } +int DaisyPythonController::get_day() const { return initialized ? daisy().time().mday() : 0; } +int DaisyPythonController::get_hour() const { return initialized ? daisy().time().hour() : 0; } +double DaisyPythonController::get_current_dt_hours() const { return 1.0; } +double DaisyPythonController::get_current_dt_days() const { return 1.0 / 24.0; } + +// ===== GROUNDWATER ===== + +double DaisyPythonController::get_groundwater_depth() const +{ + return daisy().get_groundwater_table(); +} + +bool DaisyPythonController::set_groundwater_depth(double depth_cm) +{ + daisy().set_groundwater_table(depth_cm); + return true; +} + +auto DaisyPythonController::estimate_sy_perturbation(double dh_cm) + -> std::tuple, std::vector> +{ + return daisy().estimate_sy_perturbation(dh_cm); +} + +double DaisyPythonController::get_pressure_head_at_depth(double) const +{ + return 0.0; // stub +} + +// ===== SOIL WATER ===== + +std::vector DaisyPythonController::get_soil_water_content() const { return {}; } +double DaisyPythonController::get_soil_water_at_depth(double) const { return 0.0; } +double DaisyPythonController::get_matric_potential_at_depth(double) const { return 0.0; } +double DaisyPythonController::get_cumulative_water_to_depth(double) const { return 0.0; } +double DaisyPythonController::get_available_water() const { return 0.0; } + +// ===== RECHARGE ===== + +std::vector DaisyPythonController::get_recharge_rate() const +{ + std::vector flux = daisy().get_flux_array(); + for (double& v : flux) + v *= 10.0 * 24.0; // cm/h -> mm/day + return flux; +} + +double DaisyPythonController::get_cumulative_recharge() const { return 0.0; } + +std::vector DaisyPythonController::get_pressure_head_array() const +{ + return daisy().get_h_array(); // [cm], nlayers +} + +std::vector DaisyPythonController::get_theta_array() const +{ + return daisy().get_theta_array(); // [-], nlayers +} + +std::vector DaisyPythonController::get_theta_sat_array() const +{ + return daisy().get_theta_sat_array(); // [-], nlayers, static soil param +} + + +// ===== WATER BALANCE STUBS ===== + +double DaisyPythonController::get_et_rate() const { return 0.0; } +double DaisyPythonController::get_cumulative_et() const { return 0.0; } +double DaisyPythonController::get_transpiration_rate() const { return 0.0; } +double DaisyPythonController::get_evaporation_rate() const { return 0.0; } +double DaisyPythonController::get_drainage_rate() const { return 0.0; } +double DaisyPythonController::get_cumulative_drainage() const { return 0.0; } +double DaisyPythonController::get_runoff_rate() const +{ return daisy().get_runoff_rate(); } // [mm/day] +double DaisyPythonController::get_cumulative_runoff() const { return 0.0; } + +// ===== CROP STUBS ===== + +double DaisyPythonController::get_leaf_area_index() const { return 0.0; } +double DaisyPythonController::get_root_depth() const { return 0.0; } +double DaisyPythonController::get_aboveground_biomass() const { return 0.0; } + +// ===== SOIL STUBS ===== + +double DaisyPythonController::get_soil_temperature_at_depth(double) const { return 0.0; } +double DaisyPythonController::get_soil_layer_thickness(int) const { return 0.0; } +double DaisyPythonController::get_soil_layer_top(int) const { return 0.0; } +double DaisyPythonController::get_soil_layer_bottom(int) const { return 0.0; } + +// ===== WEATHER STUBS ===== + +double DaisyPythonController::get_rainfall_rate() const { return 0.0; } +double DaisyPythonController::get_cumulative_rainfall() const { return 0.0; } +double DaisyPythonController::get_reference_et() const { return 0.0; } +double DaisyPythonController::get_air_temperature() const { return 0.0; } +double DaisyPythonController::get_air_humidity() const { return 0.0; } +double DaisyPythonController::get_wind_speed() const { return 0.0; } + +// ===== NITROGEN STUBS ===== + +double DaisyPythonController::get_mineral_nitrogen() const { return 0.0; } +double DaisyPythonController::get_cumulative_n_leaching() const { return 0.0; } +double DaisyPythonController::get_cumulative_n_uptake() const { return 0.0; } + +// ===== DIAGNOSTICS STUBS ===== + +std::string DaisyPythonController::get_last_message() const { return ""; } +bool DaisyPythonController::has_errors() const { return false; } +double DaisyPythonController::get_simulation_time() const { return 0.0; } + +// ===== CONSTRUCTOR / DESTRUCTOR ===== + +DaisyPythonController::DaisyPythonController() + : initialized(false), running(false) +{ +} + +DaisyPythonController::~DaisyPythonController() +{ + if (initialized) + { + try { finalize(); } catch (...) {} + } +} + +DaisyPythonController::DaisyPythonController(DaisyPythonController&& other) noexcept + : toplevel_(std::move(other.toplevel_)), + initialized(other.initialized), + running(other.running), + elapsed_days_(other.elapsed_days_) +{ + other.initialized = false; + other.running = false; + other.elapsed_days_ = 0.0; +} + +DaisyPythonController& DaisyPythonController::operator=(DaisyPythonController&& other) noexcept +{ + if (this != &other) + { + if (initialized) try { finalize(); } catch (...) {} + toplevel_ = std::move(other.toplevel_); + initialized = other.initialized; + running = other.running; + elapsed_days_ = other.elapsed_days_; + other.initialized = false; + other.running = false; + other.elapsed_days_ = 0.0; + } + return *this; +} + +double DaisyPythonController::get_column_area() const +{ return daisy().get_column_area(); } + +std::vector DaisyPythonController::get_layer_tops() const +{ return daisy().get_layer_tops(); } + +std::vector DaisyPythonController::get_layer_bottoms() const +{ return daisy().get_layer_bottoms(); } + +int DaisyPythonController::get_soil_layer_count() const +{ return static_cast(daisy().get_layer_tops().size()); } From 61d0c65da9d08e67a67b6b82f55079c632d57647 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 19 Jun 2026 15:58:34 +0200 Subject: [PATCH 02/11] Add daisy Python package with BMI bindings and DaisyPyFun plugin interface --- python/examples/daisy_bmi_example.py | 71 +++++++++++++++++++ python/examples/daisy_py_fun_example.py | 61 ++++++++++++++++ python/pyproject.toml | 22 ++++++ python/src/daisy/__init__.py | 19 +++++ python/src/daisy/plugins/__init__.py | 10 +++ .../daisy/plugins/daisy_py_fun_interface.py | 29 ++++++++ 6 files changed, 212 insertions(+) create mode 100644 python/examples/daisy_bmi_example.py create mode 100644 python/examples/daisy_py_fun_example.py create mode 100644 python/pyproject.toml create mode 100644 python/src/daisy/__init__.py create mode 100644 python/src/daisy/plugins/__init__.py create mode 100644 python/src/daisy/plugins/daisy_py_fun_interface.py diff --git a/python/examples/daisy_bmi_example.py b/python/examples/daisy_bmi_example.py new file mode 100644 index 000000000..f9df0dfb6 --- /dev/null +++ b/python/examples/daisy_bmi_example.py @@ -0,0 +1,71 @@ +""" +Example: Driving Daisy via the BMI interface (daisy_bmi Python module). + +This shows the standard BMI coupling pattern, and is ready to be plugged +into pymt, OpenEarth, or a custom MODFLOW6 coupling loop. + +Requirements: + - install daisy_bmi Python module via 'pip install -e python' from the root dir of this project +""" + +import os +import numpy as np +from daisy import DaisyBMI +from pathlib import Path + +# daisy config file +config_file = Path(r"c:\src\daisy\sample\sample_bmi.dai") + +# definitions +def compute_zh0(tops: np.ndarray, h_daisy: np.ndarray) -> float | None: + """Depth of h=0 [m, +down] by linear interpolation. None if WT outside column.""" + h = np.array(h_daisy) + tops_m = np.array(tops) * 1e-2 + sat_mask = h > 0.0 + if not sat_mask.any(): + return None + sat_node = int(np.argmax(sat_mask)) + if sat_node == 0: + return 0.0 + i_below = sat_node - 1 + i_above = sat_node + h_below = float(h[i_below]) + h_above = float(h[i_above]) + z_below = float(tops_m[i_below]) + z_above = float(tops_m[i_above]) + return z_below + (0.0 - h_below) / (h_above - h_below) * (z_above - z_below) + +# -- Daisy simulation -- + +# initialise Daisy +daisy = DaisyBMI() +os.chdir(config_file.parent) +daisy.initialize(str(config_file)) + +# get static internals +area = daisy.get_value("column__area") # scalar cm² +tops = daisy.get_value_array("soil_layer__top_depth") # np.NDarray[float] cm +bots = daisy.get_value_array("soil_layer__bottom_depth") # np.NDarray[float] cm +n_lay = len(tops) + +# runs simulations per daily time step +end_time = 365 + +print(f"BMI time units are: {daisy.get_time_units()}") + +for itime in np.arange(end_time-1): + # Get current time + t = daisy.get_current_time() + + # Advance Daisy one day + daisy.update_until(t + 1.0) + + # get presure heads of column + h = daisy.get_value_array("soil_water__pressure_head") + gwl = compute_zh0(tops, h_daisy=h) + if gwl is None: + print(f'gwl at day {t:.1f} is below Daisy column') + else: + print(f'gwl at day {t:.1f} is {gwl}') + +daisy.finalize() diff --git a/python/examples/daisy_py_fun_example.py b/python/examples/daisy_py_fun_example.py new file mode 100644 index 000000000..d69604731 --- /dev/null +++ b/python/examples/daisy_py_fun_example.py @@ -0,0 +1,61 @@ +'''Module illustrating how to implement a concrete DaisyPyFun that can be called from Daisy''' +import math +from daisy.plugins import DaisyPyFun + +# define Tmin plugin +class Tmin(DaisyPyFun): + '''A class that calculates the T_min function built into Daisy. See the reference manual for + details. + ''' + range = 'None' + domain = 'dg C' + + def __init__(self, args): + '''Initialize T_min. + + Parameters + ---------- + args : list + If args is not empty, the first parameter is assumed to be a reference value at which the + function value should be 1. + If args is empty, the reference value is set to 10. + ''' + if len(args) == 0: + self.ref = 10 + else: + self.ref = float(args[0]) + + + def __call__(self, t): + '''Apply Tmin to t. This is the interface exposed to Daisy + + Parameters + ---------- + t : float + Temperature in dg C + + Returns + ------- + float + ''' + value = self._apply(t) + if self.ref != 10: + scale = self._apply(self.ref) + value = math.nan if scale == 0 else value / scale + return value + + def _apply(self, t): + '''The actual application of Tmin to t. You can have as many properties and methods as you + like, as long as you implement the interface defined by DaisyPyFun''' + if t < 0: + return 0 + if t < 20: + return 1/self.ref * t + if t < 37: + return math.exp(0.47 - 0.027 * t + 0.00193 * t * t) + if t < 60: + # J.A. van Veen and M.J.Frissel. + t_max = 37 + max_val = math.exp(0.47 - 0.027 * t_max + 0.00193 * math.sqrt(t_max)) + return max_val * (1.0 - (t - 37.0) / (60.0 - 37.0)) + return 0 \ No newline at end of file diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 000000000..fde8fef89 --- /dev/null +++ b/python/pyproject.toml @@ -0,0 +1,22 @@ +[build-system] +requires = ["setuptools>=64"] +build-backend = "setuptools.build_meta" + +[project] +name = "daisy" +version = "0.1.0" +description = "Daisy BMI Python bindings" +readme = "README.md" +requires-python = ">=3.9" +authors = [ + { name = "Hendrik Kok" } +] + +# This ensures setuptools finds your package +[tool.setuptools.packages.find] +where = ["src"] +include = ["daisy*"] + +# Include compiled binaries in the package +[tool.setuptools.package-data] +"daisy" = ["*.pyd", "*.dll"] diff --git a/python/src/daisy/__init__.py b/python/src/daisy/__init__.py new file mode 100644 index 000000000..ef2a2290c --- /dev/null +++ b/python/src/daisy/__init__.py @@ -0,0 +1,19 @@ +import os as _os +import pathlib as _pathlib + +# On Windows, DLLs bundled next to the .pyd are not found automatically. +# Register the package directory so Windows finds libcore.dll and the +# other MinGW runtime DLLs that the build step copies here. +_pkg_dir = str(_pathlib.Path(__file__).parent) +if hasattr(_os, "add_dll_directory"): + _os.add_dll_directory(_pkg_dir) + +try: + from .daisy_bmi import DaisyBMI +except ImportError as e: + raise ImportError( + "Failed to import daisy_bmi. Make sure the .pyd and required DLLs are present.\n" + f" .pyd directory: {_pkg_dir}" + ) from e + +__all__ = ["DaisyBMI"] diff --git a/python/src/daisy/plugins/__init__.py b/python/src/daisy/plugins/__init__.py new file mode 100644 index 000000000..f716d7b80 --- /dev/null +++ b/python/src/daisy/plugins/__init__.py @@ -0,0 +1,10 @@ +""" +Plugin system for Daisy. + +This module exposes the base interface for defining Python plugins +that can be called from Daisy (C++ side). +""" + +from .daisy_py_fun_interface import DaisyPyFun + +__all__ = ["DaisyPyFun"] diff --git a/python/src/daisy/plugins/daisy_py_fun_interface.py b/python/src/daisy/plugins/daisy_py_fun_interface.py new file mode 100644 index 000000000..a992b2dab --- /dev/null +++ b/python/src/daisy/plugins/daisy_py_fun_interface.py @@ -0,0 +1,29 @@ +'''Abstract base class for python functions callable from Daisy. +See daisy_py_fun_example.py for a concrete usage example. +''' +from abc import ABCMeta, abstractmethod + +class DaisyPyFun(metaclass=ABCMeta): + '''A class that can be instantiated and called from Daisy''' + + @abstractmethod + def __init__(self, args: list): + '''The constructor can be used to parameterize the function. It must accept a single list as + parameter. + ''' + + @abstractmethod + def __call__(self, x: float) -> float: + '''Apply the function on x.''' + + @property + @abstractmethod + def range(self) -> str: + '''String representation of the units of the range of the function. Note that this should be + overriden by a property, not by a method.''' + + @property + @abstractmethod + def domain(self) -> str: + '''String representation of the units of the domain of the function. Note that this should be + overriden by a property, not by a method.''' From 5b0c1ff228c8c3d6d3724a1261970e6a0a8abdd6 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 19 Jun 2026 15:58:52 +0200 Subject: [PATCH 03/11] Remove old python/ files moved into src/programs/ and python/src/ --- python/CMakeLists.txt | 8 ----- python/call_python_function.cpp | 33 ----------------- python/daisy_py_fun_example.py | 61 -------------------------------- python/daisy_py_fun_interface.py | 29 --------------- 4 files changed, 131 deletions(-) delete mode 100644 python/CMakeLists.txt delete mode 100644 python/call_python_function.cpp delete mode 100644 python/daisy_py_fun_example.py delete mode 100644 python/daisy_py_fun_interface.py diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt deleted file mode 100644 index 0bbf65f9b..000000000 --- a/python/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.15) -project(call_python_from_daisy LANGUAGES CXX) - -set(PYBIND11_FINDPYTHON ON) -find_package(pybind11 REQUIRED) - -add_executable(call_python_function call_python_function.cpp) -target_link_libraries(call_python_function PRIVATE pybind11::embed) diff --git a/python/call_python_function.cpp b/python/call_python_function.cpp deleted file mode 100644 index 2db1eaf73..000000000 --- a/python/call_python_function.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include -#include -#include - -namespace py = pybind11; -using namespace py::literals; - -int main(int argc, char* argv[]) { - if (argc < 3) { - std::cout << "Usage: call_python_function " \ - "[]*" << std::endl; - return 1; - } - py::scoped_interpreter guard{}; - - auto py_module = py::module_::import(argv[1]); - - std::vector parameters; - for (std::size_t i = 3; i < argc; i++) { - parameters.push_back(argv[i]); - } - py::object py_object = py_module.attr(argv[2])(parameters); - - std::vector inputs {1, 10, 20, 30, 40}; - std::string domain = py_object.attr("domain").cast(); - std::string range = py_object.attr("range").cast(); - for (auto input : inputs) { - double py_result = py_object(input).cast(); - std::cout << "Input: " << input << " " << domain << std::endl - << "Result: " << py_result << " " << range << std::endl; - } - return 0; -} diff --git a/python/daisy_py_fun_example.py b/python/daisy_py_fun_example.py deleted file mode 100644 index a6837d526..000000000 --- a/python/daisy_py_fun_example.py +++ /dev/null @@ -1,61 +0,0 @@ -'''Module illustrating how to implement a concrete DaisyPyFun that can be called from Daisy''' -import math -from daisy_py_fun_interface import DaisyPyFun - -class Tmin(DaisyPyFun): - '''A class that calculates the T_min function built into Daisy. See the reference manual for - details. - ''' - range = 'None' - domain = 'dg C' - - - def __init__(self, args): - '''Initialize T_min. - - Parameters - ---------- - args : list - If args is not empty, the first parameter is assumed to be a reference value at which the - function value should be 1. - If args is empty, the reference value is set to 10. - ''' - if len(args) == 0: - self.ref = 10 - else: - self.ref = float(args[0]) - - - def __call__(self, t): - '''Apply Tmin to t. This is the interface exposed to Daisy - - Parameters - ---------- - t : float - Temperature in dg C - - Returns - ------- - float - ''' - value = self._apply(t) - if self.ref != 10: - scale = self._apply(self.ref) - value = math.nan if scale == 0 else value / scale - return value - - def _apply(self, t): - '''The actual application of Tmin to t. You can have as many properties and methods as you - like, as long as you implement the interface defined by DaisyPyFun''' - if t < 0: - return 0 - if t < 20: - return 1/self.ref * t - if t < 37: - return math.exp(0.47 - 0.027 * t + 0.00193 * t * t) - if t < 60: - # J.A. van Veen and M.J.Frissel. - t_max = 37 - max_val = math.exp(0.47 - 0.027 * t_max + 0.00193 * math.sqrt(t_max)) - return max_val * (1.0 - (t - 37.0) / (60.0 - 37.0)) - return 0 diff --git a/python/daisy_py_fun_interface.py b/python/daisy_py_fun_interface.py deleted file mode 100644 index a992b2dab..000000000 --- a/python/daisy_py_fun_interface.py +++ /dev/null @@ -1,29 +0,0 @@ -'''Abstract base class for python functions callable from Daisy. -See daisy_py_fun_example.py for a concrete usage example. -''' -from abc import ABCMeta, abstractmethod - -class DaisyPyFun(metaclass=ABCMeta): - '''A class that can be instantiated and called from Daisy''' - - @abstractmethod - def __init__(self, args: list): - '''The constructor can be used to parameterize the function. It must accept a single list as - parameter. - ''' - - @abstractmethod - def __call__(self, x: float) -> float: - '''Apply the function on x.''' - - @property - @abstractmethod - def range(self) -> str: - '''String representation of the units of the range of the function. Note that this should be - overriden by a property, not by a method.''' - - @property - @abstractmethod - def domain(self) -> str: - '''String representation of the units of the domain of the function. Note that this should be - overriden by a property, not by a method.''' From 8967961deb6f86c2c8583fe2d3fdf312c2b985c8 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 19 Jun 2026 15:59:22 +0200 Subject: [PATCH 04/11] Add sample_bmi.dai and build scripts --- build.bat | 18 ++++++++++ build_exe.bat | 39 +++++++++++++++++++++ sample/sample_bmi.dai | 40 ++++++++++++++++++++++ tools/CMakeLists.txt | 1 + tools/call_python/CMakeLists.txt | 8 +++++ tools/call_python/call_python_function.cpp | 33 ++++++++++++++++++ 6 files changed, 139 insertions(+) create mode 100644 build.bat create mode 100644 build_exe.bat create mode 100644 sample/sample_bmi.dai create mode 100644 tools/CMakeLists.txt create mode 100644 tools/call_python/CMakeLists.txt create mode 100644 tools/call_python/call_python_function.cpp diff --git a/build.bat b/build.bat new file mode 100644 index 000000000..77d579767 --- /dev/null +++ b/build.bat @@ -0,0 +1,18 @@ +rem set PATH=C:\msys64\ucrt64\bin;C:\msys64\usr\bin;%PATH% +rem cd C:\werkmap\PHISHIS\coupling_daisy\src\daisy-7.0.7\daisy-7.0.7 +rem cmake --preset mingw-gcc-native -B build/release +rem cmake --build build/release + +set PATH=C:\msys64\ucrt64\bin;C:\msys64\usr\bin;%PATH% +cd c:\src\daisy + +rem configure (paths defined in CMakeUserPresets.json) +cmake --preset mingw-gcc-native-local + +rem build the core +cmake --build build/release --target core + +rem rebuild only the .pyd target +cmake --build build/release --target daisy_bmi + +pause \ No newline at end of file diff --git a/build_exe.bat b/build_exe.bat new file mode 100644 index 000000000..6cbae9bd7 --- /dev/null +++ b/build_exe.bat @@ -0,0 +1,39 @@ +@echo off +:: build_exe.bat -- Configure (if needed) and build only the daisy-bin executable. +:: +:: Requires MSYS2 ucrt64 installed at C:\msys64. +:: Uses the CMake preset "mingw-gcc-native-local" defined in CMakeUserPresets.json. + +set CMAKE=C:\msys64\ucrt64\bin\cmake.exe +set SOURCE_DIR=%~dp0 +set BUILD_DIR=%~dp0build\release + +if not exist "%CMAKE%" ( + echo ERROR: cmake not found at %CMAKE% + echo Please install MSYS2 and run: pacman -S mingw-w64-ucrt-x86_64-cmake + exit /b 1 +) + +:: MSYS2 ucrt64 must be on PATH so cmake can find ninja.exe and g++.exe +set PATH=C:\msys64\ucrt64\bin;%PATH% + +:: Configure (or reconfigure) using the local preset. +:: The preset defines both source and build dirs, so just cd to the source root. +echo Configuring... +cd /d "%~dp0" +"%CMAKE%" --preset mingw-gcc-native-local +if %ERRORLEVEL% neq 0 ( + echo CONFIGURE FAILED + exit /b %ERRORLEVEL% +) + +echo Building daisy-bin... +"%CMAKE%" --build "%BUILD_DIR%" --target daisy-bin -j%NUMBER_OF_PROCESSORS% + +if %ERRORLEVEL% neq 0 ( + echo BUILD FAILED + exit /b %ERRORLEVEL% +) + +echo. +echo Build succeeded: %BUILD_DIR%\daisy-bin.exe diff --git a/sample/sample_bmi.dai b/sample/sample_bmi.dai new file mode 100644 index 000000000..3da309423 --- /dev/null +++ b/sample/sample_bmi.dai @@ -0,0 +1,40 @@ +;;; sample_bmi.dai --- Daisy setup for BMI coupling with fixed groundwater. + +;; Same inputs as sample.dai +(input file "fertilizer.dai") +(input file "tillage.dai") +(input file "crop.dai") +(input file "dk-management.dai") +(input file "irrigation.dai") +(input file "init-soil.dai") + +;; Inherit JB1_init_Cosby (30 cm Ap + 120 cm C, total 150 cm). +;; Closed bottom (zero flux) with initial water table at -130 cm. +;; Below the WT, h > 0 (saturated). Above, capillary equilibrium. +(defcolumn JB1_bmi JB1_init_Cosby + "JB1 column with closed lower boundary and initial water table at -130 cm." + (Groundwater flux 0) + (SoilWater (initial_h (-50 [cm] -80 [cm]) ; unsaturated near surface + (-100 [cm] -30 [cm]) ; approaching saturation + (-130 [cm] 0 [cm]) ; water table at -130 cm + (-150 [cm] 20 [cm])))) ; saturated at bottom + +(weather default "dk-taastrup.dwf") + +(time 1990 3 1 1) +(stop 1994 4 1 1) + +(column "JB1_bmi") + +(manager activity + (while "SBarley w. MF" + (repeat irrigate_30_tensiometer_overhead)) + "WBarley w. OF" + (while "WRape w. MF" + (activity irrigate_30_content_overhead irrigate_30_content_overhead)) + "SBarley & Pea") + +;; No output (driven externally via BMI) +(output) + +;;; sample_bmi.dai ends here diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt new file mode 100644 index 000000000..342af7232 --- /dev/null +++ b/tools/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(call_python) \ No newline at end of file diff --git a/tools/call_python/CMakeLists.txt b/tools/call_python/CMakeLists.txt new file mode 100644 index 000000000..0bbf65f9b --- /dev/null +++ b/tools/call_python/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.15) +project(call_python_from_daisy LANGUAGES CXX) + +set(PYBIND11_FINDPYTHON ON) +find_package(pybind11 REQUIRED) + +add_executable(call_python_function call_python_function.cpp) +target_link_libraries(call_python_function PRIVATE pybind11::embed) diff --git a/tools/call_python/call_python_function.cpp b/tools/call_python/call_python_function.cpp new file mode 100644 index 000000000..2db1eaf73 --- /dev/null +++ b/tools/call_python/call_python_function.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +namespace py = pybind11; +using namespace py::literals; + +int main(int argc, char* argv[]) { + if (argc < 3) { + std::cout << "Usage: call_python_function " \ + "[]*" << std::endl; + return 1; + } + py::scoped_interpreter guard{}; + + auto py_module = py::module_::import(argv[1]); + + std::vector parameters; + for (std::size_t i = 3; i < argc; i++) { + parameters.push_back(argv[i]); + } + py::object py_object = py_module.attr(argv[2])(parameters); + + std::vector inputs {1, 10, 20, 30, 40}; + std::string domain = py_object.attr("domain").cast(); + std::string range = py_object.attr("range").cast(); + for (auto input : inputs) { + double py_result = py_object(input).cast(); + std::cout << "Input: " << input << " " << domain << std::endl + << "Result: " << py_result << " " << range << std::endl; + } + return 0; +} From 637c9da1b5bfd5ba536c158fe7e285ecc978f2a4 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 19 Jun 2026 15:59:42 +0200 Subject: [PATCH 05/11] Update .gitignore: exclude build artifacts, .pyd, .dll and egg-info --- .gitignore | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index ca91f1aa7..cdfd5bb4a 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,12 @@ build TAGS # python stuff -__pycache__ \ No newline at end of file +__pycache__ +python/src/*.egg-info/ + +# build artifacts copied into the Python package by CMake +python/src/daisy/*.pyd +python/src/daisy/*.dll + +# local machine-specific CMake overrides (paths, Python envs, etc.) +CMakeUserPresets.json \ No newline at end of file From fd8e0a81d108493e37d85c7d1ead458de6be2b13 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Thu, 25 Jun 2026 14:23:40 +0200 Subject: [PATCH 06/11] Rename: DaisyPythonController -> DaisyController, bmi_bindings.cpp -> api_bindings.cpp --- include/programs/bmi.h | 4 +- include/programs/daisy_bmi.h | 16 +-- src/programs/CMakeLists.txt | 4 +- .../{bmi_bindings.cpp => api_bindings.cpp} | 4 +- src/programs/bmi.C | 6 +- src/programs/daisy_bmi.C | 136 +++++++++--------- 6 files changed, 85 insertions(+), 85 deletions(-) rename src/programs/{bmi_bindings.cpp => api_bindings.cpp} (97%) diff --git a/include/programs/bmi.h b/include/programs/bmi.h index b665adc47..16cfe94d0 100644 --- a/include/programs/bmi.h +++ b/include/programs/bmi.h @@ -15,7 +15,7 @@ /** * @class DaisyBMI - * @brief BMI 2.0 wrapper around DaisyPythonController + * @brief BMI 2.0 wrapper around DaisyController * * Variable names follow a "_" convention. * @@ -105,7 +105,7 @@ class DaisyBMI ~DaisyBMI(); private: - DaisyPythonController ctrl_; + DaisyController ctrl_; double start_time_days_; // always 0 double current_time_days_; // updated each update() double dt_days_; // timestep size in days diff --git a/include/programs/daisy_bmi.h b/include/programs/daisy_bmi.h index bc6bfe3c2..b31b213b0 100644 --- a/include/programs/daisy_bmi.h +++ b/include/programs/daisy_bmi.h @@ -20,7 +20,7 @@ class Toplevel; class Daisy; /** - * @class DaisyPythonController + * @class DaisyController * @brief Python-friendly interface to control Daisy simulation externally * * Allows Python code to: @@ -30,7 +30,7 @@ class Daisy; * - Set boundary conditions (GW depth, rainfall, irrigation, etc.) * - Control simulation parameters */ -class DaisyPythonController +class DaisyController { private: // Toplevel owns Metalib, parser, Treelog, and the Daisy program instance. @@ -457,20 +457,20 @@ class DaisyPythonController /** * Constructor */ - DaisyPythonController(); + DaisyController(); /** * Destructor */ - ~DaisyPythonController(); + ~DaisyController(); // Delete copy operations - DaisyPythonController(const DaisyPythonController&) = delete; - DaisyPythonController& operator=(const DaisyPythonController&) = delete; + DaisyController(const DaisyController&) = delete; + DaisyController& operator=(const DaisyController&) = delete; // Allow move operations - DaisyPythonController(DaisyPythonController&&) noexcept; - DaisyPythonController& operator=(DaisyPythonController&&) noexcept; + DaisyController(DaisyController&&) noexcept; + DaisyController& operator=(DaisyController&&) noexcept; }; #endif // DAISY_PYTHON_CONTROLLER_H diff --git a/src/programs/CMakeLists.txt b/src/programs/CMakeLists.txt index 76d0934e2..d89ad9bf1 100644 --- a/src/programs/CMakeLists.txt +++ b/src/programs/CMakeLists.txt @@ -20,7 +20,7 @@ target_sources(${DAISY_CORE_NAME} PRIVATE ) # ===== daisy_bmi pybind11 module ===== -pybind11_add_module(daisy_bmi bmi_bindings.cpp bmi.C daisy_bmi.C) +pybind11_add_module(daisy_bmi api_bindings.cpp bmi.C daisy_bmi.C) target_link_libraries(daisy_bmi PRIVATE ${DAISY_CORE_NAME}) target_include_directories(daisy_bmi PRIVATE ${CMAKE_SOURCE_DIR}/include @@ -82,4 +82,4 @@ if (WIN32) "${COMPILER_DIR}/libgomp-1.dll" "${PYTHON_PACKAGE_DIR}" ) -endif() \ No newline at end of file +endif() diff --git a/src/programs/bmi_bindings.cpp b/src/programs/api_bindings.cpp similarity index 97% rename from src/programs/bmi_bindings.cpp rename to src/programs/api_bindings.cpp index 94001f0ab..5ce49bd04 100644 --- a/src/programs/bmi_bindings.cpp +++ b/src/programs/api_bindings.cpp @@ -1,11 +1,11 @@ -// bmi_bindings.cpp -- pybind11 bindings for DaisyBMI +// api_bindings.cpp -- pybind11 bindings for DaisyBMI // // Exposes DaisyBMI to Python as the module "daisy_bmi". // // Compile (example, adjust paths): // g++ -O3 -Wall -shared -std=c++17 -fPIC // $(python3 -m pybind11 --includes) -// daisy_bmi_bindings.cpp daisy_bmi.C daisy_python_controller.C +// daisy_api_bindings.cpp daisy_bmi.C daisy_python_controller.C // -o daisy_bmi$(python3-config --extension-suffix) // -L. -ldaisy_core // diff --git a/src/programs/bmi.C b/src/programs/bmi.C index 8e97cb929..4614c6922 100644 --- a/src/programs/bmi.C +++ b/src/programs/bmi.C @@ -45,7 +45,7 @@ DaisyBMI::DaisyBMI() DaisyBMI::~DaisyBMI() { - // DaisyPythonController destructor handles cleanup + // DaisyController destructor handles cleanup } // ===== LIFECYCLE ===== @@ -305,9 +305,9 @@ void DaisyBMI::set_value(const std::string& name, const double* src) } if (name == "irrigation__rate") { - // TODO: expose irrigation setter in DaisyPythonController + // TODO: expose irrigation setter in DaisyController // ctrl_.set_irrigation_rate(src[0]); - throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyPythonController"); + throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyController"); } throw std::invalid_argument("Unknown input variable: " + name); } diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C index 93a50f438..f365180c8 100644 --- a/src/programs/daisy_bmi.C +++ b/src/programs/daisy_bmi.C @@ -19,12 +19,12 @@ // ===== PRIVATE HELPERS ===== -Daisy& DaisyPythonController::daisy() const +Daisy& DaisyController::daisy() const { return dynamic_cast(toplevel_->program()); } -void DaisyPythonController::setup_logging() +void DaisyController::setup_logging() { if (toplevel_) toplevel_->set_ui_none(); } @@ -66,7 +66,7 @@ static void ensure_daisy_home() catch (...) {} // non-fatal: Daisy will fall back to its own detection } -bool DaisyPythonController::load_config_file(const std::string& config_file) +bool DaisyController::load_config_file(const std::string& config_file) { ensure_daisy_home(); toplevel_ = std::make_unique("daisy"); @@ -91,7 +91,7 @@ bool DaisyPythonController::load_config_file(const std::string& config_file) // ===== INITIALIZATION ===== -bool DaisyPythonController::initialize(const std::string& config_file) +bool DaisyController::initialize(const std::string& config_file) { // Let exceptions propagate — the caller (DaisyBMI) will wrap them // with a meaningful message. A silent 'return false' hides the root cause. @@ -102,7 +102,7 @@ bool DaisyPythonController::initialize(const std::string& config_file) // ===== LIFECYCLE ===== -void DaisyPythonController::start() +void DaisyController::start() { if (!toplevel_) return; toplevel_->initialize(); // moves state from is_uninitialized -> is_ready @@ -120,7 +120,7 @@ void DaisyPythonController::start() } } -bool DaisyPythonController::advance(double days) +bool DaisyController::advance(double days) { if (!initialized) return false; try @@ -137,7 +137,7 @@ bool DaisyPythonController::advance(double days) catch (...) { return false; } } -bool DaisyPythonController::tick() +bool DaisyController::tick() { if (!initialized) return false; try @@ -150,17 +150,17 @@ bool DaisyPythonController::tick() catch (...) { return false; } } -void DaisyPythonController::stop() +void DaisyController::stop() { running = false; } -bool DaisyPythonController::is_running() const +bool DaisyController::is_running() const { return running; } -bool DaisyPythonController::finalize() +bool DaisyController::finalize() { if (!initialized) return true; daisy().stop(); @@ -174,16 +174,16 @@ bool DaisyPythonController::finalize() // ===== TIME ===== -double DaisyPythonController::get_days_since_start() const { return elapsed_days_; } +double DaisyController::get_days_since_start() const { return elapsed_days_; } -double DaisyPythonController::get_stop_days() const +double DaisyController::get_stop_days() const { if (!initialized) return std::numeric_limits::quiet_NaN(); const double h = daisy().stop_duration_hours(); return (h >= 0.0) ? h / 24.0 : std::numeric_limits::quiet_NaN(); } -std::string DaisyPythonController::get_time_string() const +std::string DaisyController::get_time_string() const { if (!initialized) return ""; const Time& t = daisy().time(); @@ -193,48 +193,48 @@ std::string DaisyPythonController::get_time_string() const return buf; } -int DaisyPythonController::get_year() const { return initialized ? daisy().time().year() : 0; } -int DaisyPythonController::get_month() const { return initialized ? daisy().time().month() : 0; } -int DaisyPythonController::get_day() const { return initialized ? daisy().time().mday() : 0; } -int DaisyPythonController::get_hour() const { return initialized ? daisy().time().hour() : 0; } -double DaisyPythonController::get_current_dt_hours() const { return 1.0; } -double DaisyPythonController::get_current_dt_days() const { return 1.0 / 24.0; } +int DaisyController::get_year() const { return initialized ? daisy().time().year() : 0; } +int DaisyController::get_month() const { return initialized ? daisy().time().month() : 0; } +int DaisyController::get_day() const { return initialized ? daisy().time().mday() : 0; } +int DaisyController::get_hour() const { return initialized ? daisy().time().hour() : 0; } +double DaisyController::get_current_dt_hours() const { return 1.0; } +double DaisyController::get_current_dt_days() const { return 1.0 / 24.0; } // ===== GROUNDWATER ===== -double DaisyPythonController::get_groundwater_depth() const +double DaisyController::get_groundwater_depth() const { return daisy().get_groundwater_table(); } -bool DaisyPythonController::set_groundwater_depth(double depth_cm) +bool DaisyController::set_groundwater_depth(double depth_cm) { daisy().set_groundwater_table(depth_cm); return true; } -auto DaisyPythonController::estimate_sy_perturbation(double dh_cm) +auto DaisyController::estimate_sy_perturbation(double dh_cm) -> std::tuple, std::vector> { return daisy().estimate_sy_perturbation(dh_cm); } -double DaisyPythonController::get_pressure_head_at_depth(double) const +double DaisyController::get_pressure_head_at_depth(double) const { return 0.0; // stub } // ===== SOIL WATER ===== -std::vector DaisyPythonController::get_soil_water_content() const { return {}; } -double DaisyPythonController::get_soil_water_at_depth(double) const { return 0.0; } -double DaisyPythonController::get_matric_potential_at_depth(double) const { return 0.0; } -double DaisyPythonController::get_cumulative_water_to_depth(double) const { return 0.0; } -double DaisyPythonController::get_available_water() const { return 0.0; } +std::vector DaisyController::get_soil_water_content() const { return {}; } +double DaisyController::get_soil_water_at_depth(double) const { return 0.0; } +double DaisyController::get_matric_potential_at_depth(double) const { return 0.0; } +double DaisyController::get_cumulative_water_to_depth(double) const { return 0.0; } +double DaisyController::get_available_water() const { return 0.0; } // ===== RECHARGE ===== -std::vector DaisyPythonController::get_recharge_rate() const +std::vector DaisyController::get_recharge_rate() const { std::vector flux = daisy().get_flux_array(); for (double& v : flux) @@ -242,19 +242,19 @@ std::vector DaisyPythonController::get_recharge_rate() const return flux; } -double DaisyPythonController::get_cumulative_recharge() const { return 0.0; } +double DaisyController::get_cumulative_recharge() const { return 0.0; } -std::vector DaisyPythonController::get_pressure_head_array() const +std::vector DaisyController::get_pressure_head_array() const { return daisy().get_h_array(); // [cm], nlayers } -std::vector DaisyPythonController::get_theta_array() const +std::vector DaisyController::get_theta_array() const { return daisy().get_theta_array(); // [-], nlayers } -std::vector DaisyPythonController::get_theta_sat_array() const +std::vector DaisyController::get_theta_sat_array() const { return daisy().get_theta_sat_array(); // [-], nlayers, static soil param } @@ -262,58 +262,58 @@ std::vector DaisyPythonController::get_theta_sat_array() const // ===== WATER BALANCE STUBS ===== -double DaisyPythonController::get_et_rate() const { return 0.0; } -double DaisyPythonController::get_cumulative_et() const { return 0.0; } -double DaisyPythonController::get_transpiration_rate() const { return 0.0; } -double DaisyPythonController::get_evaporation_rate() const { return 0.0; } -double DaisyPythonController::get_drainage_rate() const { return 0.0; } -double DaisyPythonController::get_cumulative_drainage() const { return 0.0; } -double DaisyPythonController::get_runoff_rate() const +double DaisyController::get_et_rate() const { return 0.0; } +double DaisyController::get_cumulative_et() const { return 0.0; } +double DaisyController::get_transpiration_rate() const { return 0.0; } +double DaisyController::get_evaporation_rate() const { return 0.0; } +double DaisyController::get_drainage_rate() const { return 0.0; } +double DaisyController::get_cumulative_drainage() const { return 0.0; } +double DaisyController::get_runoff_rate() const { return daisy().get_runoff_rate(); } // [mm/day] -double DaisyPythonController::get_cumulative_runoff() const { return 0.0; } +double DaisyController::get_cumulative_runoff() const { return 0.0; } // ===== CROP STUBS ===== -double DaisyPythonController::get_leaf_area_index() const { return 0.0; } -double DaisyPythonController::get_root_depth() const { return 0.0; } -double DaisyPythonController::get_aboveground_biomass() const { return 0.0; } +double DaisyController::get_leaf_area_index() const { return 0.0; } +double DaisyController::get_root_depth() const { return 0.0; } +double DaisyController::get_aboveground_biomass() const { return 0.0; } // ===== SOIL STUBS ===== -double DaisyPythonController::get_soil_temperature_at_depth(double) const { return 0.0; } -double DaisyPythonController::get_soil_layer_thickness(int) const { return 0.0; } -double DaisyPythonController::get_soil_layer_top(int) const { return 0.0; } -double DaisyPythonController::get_soil_layer_bottom(int) const { return 0.0; } +double DaisyController::get_soil_temperature_at_depth(double) const { return 0.0; } +double DaisyController::get_soil_layer_thickness(int) const { return 0.0; } +double DaisyController::get_soil_layer_top(int) const { return 0.0; } +double DaisyController::get_soil_layer_bottom(int) const { return 0.0; } // ===== WEATHER STUBS ===== -double DaisyPythonController::get_rainfall_rate() const { return 0.0; } -double DaisyPythonController::get_cumulative_rainfall() const { return 0.0; } -double DaisyPythonController::get_reference_et() const { return 0.0; } -double DaisyPythonController::get_air_temperature() const { return 0.0; } -double DaisyPythonController::get_air_humidity() const { return 0.0; } -double DaisyPythonController::get_wind_speed() const { return 0.0; } +double DaisyController::get_rainfall_rate() const { return 0.0; } +double DaisyController::get_cumulative_rainfall() const { return 0.0; } +double DaisyController::get_reference_et() const { return 0.0; } +double DaisyController::get_air_temperature() const { return 0.0; } +double DaisyController::get_air_humidity() const { return 0.0; } +double DaisyController::get_wind_speed() const { return 0.0; } // ===== NITROGEN STUBS ===== -double DaisyPythonController::get_mineral_nitrogen() const { return 0.0; } -double DaisyPythonController::get_cumulative_n_leaching() const { return 0.0; } -double DaisyPythonController::get_cumulative_n_uptake() const { return 0.0; } +double DaisyController::get_mineral_nitrogen() const { return 0.0; } +double DaisyController::get_cumulative_n_leaching() const { return 0.0; } +double DaisyController::get_cumulative_n_uptake() const { return 0.0; } // ===== DIAGNOSTICS STUBS ===== -std::string DaisyPythonController::get_last_message() const { return ""; } -bool DaisyPythonController::has_errors() const { return false; } -double DaisyPythonController::get_simulation_time() const { return 0.0; } +std::string DaisyController::get_last_message() const { return ""; } +bool DaisyController::has_errors() const { return false; } +double DaisyController::get_simulation_time() const { return 0.0; } // ===== CONSTRUCTOR / DESTRUCTOR ===== -DaisyPythonController::DaisyPythonController() +DaisyController::DaisyController() : initialized(false), running(false) { } -DaisyPythonController::~DaisyPythonController() +DaisyController::~DaisyController() { if (initialized) { @@ -321,7 +321,7 @@ DaisyPythonController::~DaisyPythonController() } } -DaisyPythonController::DaisyPythonController(DaisyPythonController&& other) noexcept +DaisyController::DaisyController(DaisyController&& other) noexcept : toplevel_(std::move(other.toplevel_)), initialized(other.initialized), running(other.running), @@ -332,7 +332,7 @@ DaisyPythonController::DaisyPythonController(DaisyPythonController&& other) noex other.elapsed_days_ = 0.0; } -DaisyPythonController& DaisyPythonController::operator=(DaisyPythonController&& other) noexcept +DaisyController& DaisyController::operator=(DaisyController&& other) noexcept { if (this != &other) { @@ -348,14 +348,14 @@ DaisyPythonController& DaisyPythonController::operator=(DaisyPythonController&& return *this; } -double DaisyPythonController::get_column_area() const +double DaisyController::get_column_area() const { return daisy().get_column_area(); } -std::vector DaisyPythonController::get_layer_tops() const +std::vector DaisyController::get_layer_tops() const { return daisy().get_layer_tops(); } -std::vector DaisyPythonController::get_layer_bottoms() const +std::vector DaisyController::get_layer_bottoms() const { return daisy().get_layer_bottoms(); } -int DaisyPythonController::get_soil_layer_count() const +int DaisyController::get_soil_layer_count() const { return static_cast(daisy().get_layer_tops().size()); } From 9fb5c3a04dcbf6e88fdb1a0a38bcd65892c787fd Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Thu, 25 Jun 2026 15:05:05 +0200 Subject: [PATCH 07/11] update python package on new name --- .vscode/settings.json | 3 +++ python/examples/daisy_bmi_example.py | 4 ++-- python/src/daisy/__init__.py | 4 ++-- src/programs/api_bindings.cpp | 2 +- 4 files changed, 8 insertions(+), 5 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..c9ebf2d27 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:system" +} \ No newline at end of file diff --git a/python/examples/daisy_bmi_example.py b/python/examples/daisy_bmi_example.py index f9df0dfb6..e0b3b4370 100644 --- a/python/examples/daisy_bmi_example.py +++ b/python/examples/daisy_bmi_example.py @@ -10,7 +10,7 @@ import os import numpy as np -from daisy import DaisyBMI +from daisy import DaisyAPI from pathlib import Path # daisy config file @@ -38,7 +38,7 @@ def compute_zh0(tops: np.ndarray, h_daisy: np.ndarray) -> float | None: # -- Daisy simulation -- # initialise Daisy -daisy = DaisyBMI() +daisy = DaisyAPI() os.chdir(config_file.parent) daisy.initialize(str(config_file)) diff --git a/python/src/daisy/__init__.py b/python/src/daisy/__init__.py index ef2a2290c..baf39099c 100644 --- a/python/src/daisy/__init__.py +++ b/python/src/daisy/__init__.py @@ -9,11 +9,11 @@ _os.add_dll_directory(_pkg_dir) try: - from .daisy_bmi import DaisyBMI + from .daisy_bmi import DaisyAPI except ImportError as e: raise ImportError( "Failed to import daisy_bmi. Make sure the .pyd and required DLLs are present.\n" f" .pyd directory: {_pkg_dir}" ) from e -__all__ = ["DaisyBMI"] +__all__ = ["DaisyAPI"] diff --git a/src/programs/api_bindings.cpp b/src/programs/api_bindings.cpp index 5ce49bd04..514ba565f 100644 --- a/src/programs/api_bindings.cpp +++ b/src/programs/api_bindings.cpp @@ -57,7 +57,7 @@ PYBIND11_MODULE(daisy_bmi, m) bmi.finalize() )pbdoc"; - py::class_(m, "DaisyBMI") + py::class_(m, "DaisyAPI") .def(py::init<>(), "Create a new Daisy BMI instance") From fc1b93f5edd256aa6443321ec3200a8dfc795759 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Thu, 25 Jun 2026 17:03:02 +0200 Subject: [PATCH 08/11] cleaning up --- include/programs/daisy_bmi.h | 8 +++++++- src/programs/daisy_bmi.C | 5 +++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/include/programs/daisy_bmi.h b/include/programs/daisy_bmi.h index b31b213b0..eaad2df1d 100644 --- a/include/programs/daisy_bmi.h +++ b/include/programs/daisy_bmi.h @@ -471,6 +471,12 @@ class DaisyController // Allow move operations DaisyController(DaisyController&&) noexcept; DaisyController& operator=(DaisyController&&) noexcept; + +public: + /** Direct access to the Daisy simulation object. + * For use by DaisyAPI extension methods — not intended for BMI callers. + * Public because DaisyAPI does not inherit DaisyController. */ + Daisy& daisy_ref (); }; -#endif // DAISY_PYTHON_CONTROLLER_H +#endif // DAISY_BMI_H diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C index f365180c8..1d7812a81 100644 --- a/src/programs/daisy_bmi.C +++ b/src/programs/daisy_bmi.C @@ -219,6 +219,11 @@ auto DaisyController::estimate_sy_perturbation(double dh_cm) return daisy().estimate_sy_perturbation(dh_cm); } +Daisy& DaisyController::daisy_ref () +{ + return daisy (); +} + double DaisyController::get_pressure_head_at_depth(double) const { return 0.0; // stub From a9c18293e9bf4fec2850aac64bad624c4ebb7aa4 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Fri, 26 Jun 2026 13:17:21 +0200 Subject: [PATCH 09/11] cleaning up --- .vscode/c_cpp_properties.json | 10 ++ CMakeLists.txt | 1 + include/programs/bmi.h | 19 ++-- include/programs/daisy_bmi.h | 18 ++-- python/examples/daisy_bmi_example.py | 24 ++--- python/src/daisy/__init__.py | 4 +- src/programs/api_bindings.cpp | 59 ++++++------ src/programs/bmi.C | 77 ++++++++------- src/programs/daisy_bmi.C | 138 +++++++++++++-------------- 9 files changed, 187 insertions(+), 163 deletions(-) create mode 100644 .vscode/c_cpp_properties.json diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 000000000..2134e338c --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,10 @@ +{ + "configurations": [ + { + "name": "MinGW", + "compileCommands": "${workspaceFolder}/build/release/compile_commands.json", + "intelliSenseMode": "gcc-x64" + } + ], + "version": 4 +} diff --git a/CMakeLists.txt b/CMakeLists.txt index ff7668b4f..a7cfe039f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,7 @@ project( HOMEPAGE_URL https://daisy.ku.dk/ LANGUAGES CXX C ) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) include(cmake/AddCoverageBuildType.cmake) if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) diff --git a/include/programs/bmi.h b/include/programs/bmi.h index 16cfe94d0..095b79b37 100644 --- a/include/programs/bmi.h +++ b/include/programs/bmi.h @@ -13,9 +13,11 @@ #include #include "programs/daisy_bmi.h" +class Daisy; // forward declaration for protected daisy() accessor + /** - * @class DaisyBMI - * @brief BMI 2.0 wrapper around DaisyController + * @class BMI + * @brief BMI 2.0 wrapper around DaisyBMI * * Variable names follow a "_" convention. * @@ -37,7 +39,7 @@ * "vegetation__aboveground_biomass" [g/m2] * "soil_water__content" [m3/m3] (first layer) */ -class DaisyBMI +class BMI { public: // ===== BMI LIFECYCLE ===== @@ -101,11 +103,11 @@ class DaisyBMI estimate_sy_perturbation(double dh_cm = 1.0, unsigned int col = 0u); // ===== CONSTRUCTOR / DESTRUCTOR ===== - DaisyBMI(); - ~DaisyBMI(); + BMI(); + ~BMI(); private: - DaisyController ctrl_; + DaisyBMI ctrl_; double start_time_days_; // always 0 double current_time_days_; // updated each update() double dt_days_; // timestep size in days @@ -114,6 +116,11 @@ class DaisyBMI static const std::vector OUTPUT_VARS; double get_output_value(const std::string& name) const; + +protected: + /** Direct access to the Daisy simulation object for use by DaisyAPI + * extension methods. Same Daisy& that DaisyBMI uses internally. */ + Daisy& daisy (); }; #endif // DAISY_BMI_H diff --git a/include/programs/daisy_bmi.h b/include/programs/daisy_bmi.h index eaad2df1d..46a8b9d65 100644 --- a/include/programs/daisy_bmi.h +++ b/include/programs/daisy_bmi.h @@ -20,7 +20,7 @@ class Toplevel; class Daisy; /** - * @class DaisyController + * @class DaisyBMI * @brief Python-friendly interface to control Daisy simulation externally * * Allows Python code to: @@ -30,7 +30,7 @@ class Daisy; * - Set boundary conditions (GW depth, rainfall, irrigation, etc.) * - Control simulation parameters */ -class DaisyController +class DaisyBMI { private: // Toplevel owns Metalib, parser, Treelog, and the Daisy program instance. @@ -457,25 +457,25 @@ class DaisyController /** * Constructor */ - DaisyController(); + DaisyBMI(); /** * Destructor */ - ~DaisyController(); + ~DaisyBMI(); // Delete copy operations - DaisyController(const DaisyController&) = delete; - DaisyController& operator=(const DaisyController&) = delete; + DaisyBMI(const DaisyBMI&) = delete; + DaisyBMI& operator=(const DaisyBMI&) = delete; // Allow move operations - DaisyController(DaisyController&&) noexcept; - DaisyController& operator=(DaisyController&&) noexcept; + DaisyBMI(DaisyBMI&&) noexcept; + DaisyBMI& operator=(DaisyBMI&&) noexcept; public: /** Direct access to the Daisy simulation object. * For use by DaisyAPI extension methods — not intended for BMI callers. - * Public because DaisyAPI does not inherit DaisyController. */ + * Public because DaisyAPI does not inherit DaisyBMI. */ Daisy& daisy_ref (); }; diff --git a/python/examples/daisy_bmi_example.py b/python/examples/daisy_bmi_example.py index e0b3b4370..912b92a30 100644 --- a/python/examples/daisy_bmi_example.py +++ b/python/examples/daisy_bmi_example.py @@ -10,7 +10,7 @@ import os import numpy as np -from daisy import DaisyAPI +from daisy import API from pathlib import Path # daisy config file @@ -38,34 +38,34 @@ def compute_zh0(tops: np.ndarray, h_daisy: np.ndarray) -> float | None: # -- Daisy simulation -- # initialise Daisy -daisy = DaisyAPI() +api = API() os.chdir(config_file.parent) -daisy.initialize(str(config_file)) +api.initialize(str(config_file)) # get static internals -area = daisy.get_value("column__area") # scalar cm² -tops = daisy.get_value_array("soil_layer__top_depth") # np.NDarray[float] cm -bots = daisy.get_value_array("soil_layer__bottom_depth") # np.NDarray[float] cm +area = api.get_value("column__area") # scalar cm² +tops = api.get_value_array("soil_layer__top_depth") # np.NDarray[float] cm +bots = api.get_value_array("soil_layer__bottom_depth") # np.NDarray[float] cm n_lay = len(tops) # runs simulations per daily time step end_time = 365 -print(f"BMI time units are: {daisy.get_time_units()}") +print(f"BMI time units are: {api.get_time_units()}") for itime in np.arange(end_time-1): # Get current time - t = daisy.get_current_time() + t = api.get_current_time() # Advance Daisy one day - daisy.update_until(t + 1.0) + api.update_until(t + 1.0) - # get presure heads of column - h = daisy.get_value_array("soil_water__pressure_head") + # get pressure heads and GW table depth + h = api.get_value_array("soil_water__pressure_head") gwl = compute_zh0(tops, h_daisy=h) if gwl is None: print(f'gwl at day {t:.1f} is below Daisy column') else: print(f'gwl at day {t:.1f} is {gwl}') -daisy.finalize() +api.finalize() diff --git a/python/src/daisy/__init__.py b/python/src/daisy/__init__.py index baf39099c..a1740f3ed 100644 --- a/python/src/daisy/__init__.py +++ b/python/src/daisy/__init__.py @@ -9,11 +9,11 @@ _os.add_dll_directory(_pkg_dir) try: - from .daisy_bmi import DaisyAPI + from .daisy_bmi import BMI, API except ImportError as e: raise ImportError( "Failed to import daisy_bmi. Make sure the .pyd and required DLLs are present.\n" f" .pyd directory: {_pkg_dir}" ) from e -__all__ = ["DaisyAPI"] +__all__ = ["BMI", "API"] diff --git a/src/programs/api_bindings.cpp b/src/programs/api_bindings.cpp index 514ba565f..0ba77e5fe 100644 --- a/src/programs/api_bindings.cpp +++ b/src/programs/api_bindings.cpp @@ -1,6 +1,9 @@ -// api_bindings.cpp -- pybind11 bindings for DaisyBMI +// api_bindings.cpp -- pybind11 bindings for BMI and DaisyAPI. // -// Exposes DaisyBMI to Python as the module "daisy_bmi". +// Exposes both classes to Python as the module "daisy_bmi". +// BMI — pure BMI 2.0 standard interface. +// DaisyAPI — inherits BMI, adds Daisy-specific extensions +// (perturbation_tick, etc.). b3897eb3 (cleaning up) // // Compile (example, adjust paths): // g++ -O3 -Wall -shared -std=c++17 -fPIC @@ -36,7 +39,7 @@ PYBIND11_MODULE(daisy_bmi, m) ------- import daisy_bmi - bmi = daisy_bmi.DaisyBMI() + bmi = daisy_bmi.BMI() bmi.initialize("myconfig.dai") print("Component:", bmi.get_component_name()) @@ -57,53 +60,53 @@ PYBIND11_MODULE(daisy_bmi, m) bmi.finalize() )pbdoc"; - py::class_(m, "DaisyAPI") + py::class_(m, "DaisyBMI") .def(py::init<>(), "Create a new Daisy BMI instance") // --- Lifecycle --- - .def("initialize", &DaisyBMI::initialize, py::arg("config_file"), + .def("initialize", &BMI::initialize, py::arg("config_file"), "Initialize from a .dai config file") - .def("update", &DaisyBMI::update, + .def("update", &BMI::update, "Advance model by one timestep") - .def("update_until", &DaisyBMI::update_until, py::arg("time"), + .def("update_until", &BMI::update_until, py::arg("time"), "Advance model until time (days since start)") - .def("finalize", &DaisyBMI::finalize, + .def("finalize", &BMI::finalize, "Finalize and release all resources") // --- Information --- - .def("get_component_name", &DaisyBMI::get_component_name) - .def("get_input_var_names", &DaisyBMI::get_input_var_names) - .def("get_output_var_names", &DaisyBMI::get_output_var_names) + .def("get_component_name", &BMI::get_component_name) + .def("get_input_var_names", &BMI::get_input_var_names) + .def("get_output_var_names", &BMI::get_output_var_names) // --- Time --- - .def("get_start_time", &DaisyBMI::get_start_time, + .def("get_start_time", &BMI::get_start_time, "Simulation start time [days]") - .def("get_end_time", &DaisyBMI::get_end_time, + .def("get_end_time", &BMI::get_end_time, "Simulation end time [days] (NaN = open-ended)") - .def("get_current_time", &DaisyBMI::get_current_time, + .def("get_current_time", &BMI::get_current_time, "Current simulation time [days since start]") - .def("get_time_step", &DaisyBMI::get_time_step, + .def("get_time_step", &BMI::get_time_step, "Current timestep size [days]") - .def("get_time_units", &DaisyBMI::get_time_units, + .def("get_time_units", &BMI::get_time_units, "Time unit string, e.g. 'd'") // --- Variable metadata --- - .def("get_var_type", &DaisyBMI::get_var_type, py::arg("name")) - .def("get_var_units", &DaisyBMI::get_var_units, py::arg("name")) - .def("get_var_itemsize", &DaisyBMI::get_var_itemsize, py::arg("name")) - .def("get_var_nbytes", &DaisyBMI::get_var_nbytes, py::arg("name")) - .def("get_var_location", &DaisyBMI::get_var_location, py::arg("name")) - .def("get_var_grid", &DaisyBMI::get_var_grid, py::arg("name")) + .def("get_var_type", &BMI::get_var_type, py::arg("name")) + .def("get_var_units", &BMI::get_var_units, py::arg("name")) + .def("get_var_itemsize", &BMI::get_var_itemsize, py::arg("name")) + .def("get_var_nbytes", &BMI::get_var_nbytes, py::arg("name")) + .def("get_var_location", &BMI::get_var_location, py::arg("name")) + .def("get_var_grid", &BMI::get_var_grid, py::arg("name")) // --- Grid metadata --- - .def("get_grid_rank", &DaisyBMI::get_grid_rank, py::arg("grid")) - .def("get_grid_size", &DaisyBMI::get_grid_size, py::arg("grid")) - .def("get_grid_type", &DaisyBMI::get_grid_type, py::arg("grid")) + .def("get_grid_rank", &BMI::get_grid_rank, py::arg("grid")) + .def("get_grid_size", &BMI::get_grid_size, py::arg("grid")) + .def("get_grid_type", &BMI::get_grid_type, py::arg("grid")) // --- Get / Set (scalar convenience wrappers) --- .def("get_value", - [](const DaisyBMI& self, const std::string& name) -> double { + [](const BMI& self, const std::string& name) -> double { // Allocate the correct-sized buffer; array vars (e.g. // soil_water__recharge_rate) hold one value per soil layer. // Writing N elements into a 1-element stack variable is UB and @@ -117,14 +120,14 @@ PYBIND11_MODULE(daisy_bmi, m) "Get a scalar output value by BMI variable name") .def("set_value", - [](DaisyBMI& self, const std::string& name, double value) { + [](BMI& self, const std::string& name, double value) { self.set_value(name, &value); }, py::arg("name"), py::arg("value"), "Set a scalar input value by BMI variable name") .def("get_value_array", - [](const DaisyBMI& self, const std::string& name) { + [](const BMI& self, const std::string& name) { int n = self.get_var_nbytes(name) / static_cast(sizeof(double)); py::array_t arr(n); self.get_value(name, arr.mutable_data()); diff --git a/src/programs/bmi.C b/src/programs/bmi.C index 4614c6922..15039452f 100644 --- a/src/programs/bmi.C +++ b/src/programs/bmi.C @@ -3,18 +3,21 @@ // This file is part of Daisy. #include "programs/bmi.h" +#include "daisy/daisy.h" #include #include #include // ===== STATIC VARIABLE LISTS ===== -const std::vector DaisyBMI::INPUT_VARS = { +Daisy& BMI::daisy () { return ctrl_.daisy_ref (); } + +const std::vector BMI::INPUT_VARS = { "groundwater__depth", // [cm] below surface "irrigation__rate", // [mm/day] }; -const std::vector DaisyBMI::OUTPUT_VARS = { +const std::vector BMI::OUTPUT_VARS = { "soil_water__recharge_rate", // [mm/day] "land_surface__evapotranspiration_rate", // [mm/day] "soil_water__transpiration_rate", // [mm/day] @@ -36,21 +39,21 @@ const std::vector DaisyBMI::OUTPUT_VARS = { // ===== CONSTRUCTOR / DESTRUCTOR ===== -DaisyBMI::DaisyBMI() +BMI::BMI() : start_time_days_(0.0) , current_time_days_(0.0) , dt_days_(1.0 / 24.0) // default: hourly timestep { } -DaisyBMI::~DaisyBMI() +BMI::~BMI() { - // DaisyController destructor handles cleanup + // DaisyBMI destructor handles cleanup } // ===== LIFECYCLE ===== -void DaisyBMI::initialize(const std::string& config_file) +void BMI::initialize(const std::string& config_file) { try { @@ -59,32 +62,32 @@ void DaisyBMI::initialize(const std::string& config_file) catch (const std::exception& e) { throw std::runtime_error( - std::string("DaisyBMI::initialize failed for config: ") + config_file + std::string("BMI::initialize failed for config: ") + config_file + "\n Reason: " + e.what()); } catch (const int code) { throw std::runtime_error( - std::string("DaisyBMI::initialize failed for config: ") + config_file + std::string("BMI::initialize failed for config: ") + config_file + "\n Reason: Daisy reported an error (check daisy.log for details). Exit code: " + std::to_string(code)); } catch (const char* msg) { throw std::runtime_error( - std::string("DaisyBMI::initialize failed for config: ") + config_file + std::string("BMI::initialize failed for config: ") + config_file + "\n Reason: " + msg); } catch (...) { throw std::runtime_error( - std::string("DaisyBMI::initialize failed for config: ") + config_file + std::string("BMI::initialize failed for config: ") + config_file + "\n Reason: unknown exception"); } if (!ctrl_.is_initialized()) throw std::runtime_error( - std::string("DaisyBMI::initialize failed for config: ") + config_file); + std::string("BMI::initialize failed for config: ") + config_file); // Use Daisy's actual internal time as the reference point. // get_days_since_start() returns 0 right after init, but reading it through @@ -97,10 +100,10 @@ void DaisyBMI::initialize(const std::string& config_file) dt_days_ = 1.0 / 24.0; // fallback: hourly } -void DaisyBMI::update() +void BMI::update() { if (!ctrl_.tick()) - throw std::runtime_error("DaisyBMI::update (tick) failed"); + throw std::runtime_error("BMI::update (tick) failed"); // Sync current time directly from Daisy's internal clock. current_time_days_ = ctrl_.get_days_since_start(); @@ -111,7 +114,7 @@ void DaisyBMI::update() dt_days_ = actual_dt; } -void DaisyBMI::update_until(double time) +void BMI::update_until(double time) { // Use half-timestep epsilon to avoid running one tick too many due to // floating-point accumulation (e.g. 24 * 1/24 is not exactly 1.0). @@ -125,63 +128,63 @@ void DaisyBMI::update_until(double time) current_time_days_ = time; } -void DaisyBMI::finalize() +void BMI::finalize() { ctrl_.finalize(); } // ===== INFORMATION ===== -std::string DaisyBMI::get_component_name() const +std::string BMI::get_component_name() const { return "Daisy"; } -std::vector DaisyBMI::get_input_var_names() const +std::vector BMI::get_input_var_names() const { return INPUT_VARS; } -std::vector DaisyBMI::get_output_var_names() const +std::vector BMI::get_output_var_names() const { return OUTPUT_VARS; } // ===== TIME ===== -double DaisyBMI::get_start_time() const +double BMI::get_start_time() const { return start_time_days_; } -double DaisyBMI::get_end_time() const +double BMI::get_end_time() const { return ctrl_.get_stop_days(); } -double DaisyBMI::get_current_time() const +double BMI::get_current_time() const { return current_time_days_; } -double DaisyBMI::get_time_step() const +double BMI::get_time_step() const { return dt_days_; } -std::string DaisyBMI::get_time_units() const +std::string BMI::get_time_units() const { return "d"; } // ===== VARIABLE INFO ===== -std::string DaisyBMI::get_var_type(const std::string& /*name*/) const +std::string BMI::get_var_type(const std::string& /*name*/) const { return "double"; } -std::string DaisyBMI::get_var_units(const std::string& name) const +std::string BMI::get_var_units(const std::string& name) const { if (name == "groundwater__depth") return "cm"; if (name == "irrigation__rate") return "mm d-1"; @@ -204,7 +207,7 @@ std::string DaisyBMI::get_var_units(const std::string& name) const throw std::invalid_argument("Unknown variable: " + name); } -int DaisyBMI::get_var_nbytes(const std::string& name) const { +int BMI::get_var_nbytes(const std::string& name) const { if (name == "soil_layer__top_depth" || name == "soil_layer__bottom_depth" || name == "soil_water__recharge_rate" || name == "soil_water__pressure_head" @@ -214,30 +217,30 @@ int DaisyBMI::get_var_nbytes(const std::string& name) const { return sizeof(double); } -int DaisyBMI::get_var_itemsize(const std::string& /*name*/) const +int BMI::get_var_itemsize(const std::string& /*name*/) const { return static_cast(sizeof(double)); } -std::string DaisyBMI::get_var_location(const std::string& /*name*/) const +std::string BMI::get_var_location(const std::string& /*name*/) const { return "none"; // scalars have no grid location } -int DaisyBMI::get_var_grid(const std::string& /*name*/) const +int BMI::get_var_grid(const std::string& /*name*/) const { return 0; // all vars on the single-cell scalar grid } // ===== GRID INFO ===== -int DaisyBMI::get_grid_rank(int /*grid*/) const { return 0; } -int DaisyBMI::get_grid_size(int /*grid*/) const { return 1; } -std::string DaisyBMI::get_grid_type(int /*grid*/) const { return "scalar"; } +int BMI::get_grid_rank(int /*grid*/) const { return 0; } +int BMI::get_grid_size(int /*grid*/) const { return 1; } +std::string BMI::get_grid_type(int /*grid*/) const { return "scalar"; } // ===== GET / SET ===== -double DaisyBMI::get_output_value(const std::string& name) const +double BMI::get_output_value(const std::string& name) const { if (name == "soil_water__recharge_rate") { auto v = ctrl_.get_recharge_rate(); @@ -261,7 +264,7 @@ double DaisyBMI::get_output_value(const std::string& name) const throw std::invalid_argument("Unknown output variable: " + name); } -void DaisyBMI::get_value(const std::string& name, double* dest) const +void BMI::get_value(const std::string& name, double* dest) const { if (name == "soil_layer__top_depth") { auto v = ctrl_.get_layer_tops(); @@ -296,7 +299,7 @@ void DaisyBMI::get_value(const std::string& name, double* dest) const dest[0] = get_output_value(name); // bestaand pad voor scalars } -void DaisyBMI::set_value(const std::string& name, const double* src) +void BMI::set_value(const std::string& name, const double* src) { if (name == "groundwater__depth") { @@ -305,9 +308,9 @@ void DaisyBMI::set_value(const std::string& name, const double* src) } if (name == "irrigation__rate") { - // TODO: expose irrigation setter in DaisyController + // TODO: expose irrigation setter in DaisyBMI // ctrl_.set_irrigation_rate(src[0]); - throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyController"); + throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyBMI"); } throw std::invalid_argument("Unknown input variable: " + name); } diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C index 1d7812a81..6693497fd 100644 --- a/src/programs/daisy_bmi.C +++ b/src/programs/daisy_bmi.C @@ -19,12 +19,12 @@ // ===== PRIVATE HELPERS ===== -Daisy& DaisyController::daisy() const +Daisy& DaisyBMI::daisy() const { return dynamic_cast(toplevel_->program()); } -void DaisyController::setup_logging() +void DaisyBMI::setup_logging() { if (toplevel_) toplevel_->set_ui_none(); } @@ -66,7 +66,7 @@ static void ensure_daisy_home() catch (...) {} // non-fatal: Daisy will fall back to its own detection } -bool DaisyController::load_config_file(const std::string& config_file) +bool DaisyBMI::load_config_file(const std::string& config_file) { ensure_daisy_home(); toplevel_ = std::make_unique("daisy"); @@ -91,7 +91,7 @@ bool DaisyController::load_config_file(const std::string& config_file) // ===== INITIALIZATION ===== -bool DaisyController::initialize(const std::string& config_file) +bool DaisyBMI::initialize(const std::string& config_file) { // Let exceptions propagate — the caller (DaisyBMI) will wrap them // with a meaningful message. A silent 'return false' hides the root cause. @@ -102,7 +102,7 @@ bool DaisyController::initialize(const std::string& config_file) // ===== LIFECYCLE ===== -void DaisyController::start() +void DaisyBMI::start() { if (!toplevel_) return; toplevel_->initialize(); // moves state from is_uninitialized -> is_ready @@ -120,7 +120,7 @@ void DaisyController::start() } } -bool DaisyController::advance(double days) +bool DaisyBMI::advance(double days) { if (!initialized) return false; try @@ -137,7 +137,7 @@ bool DaisyController::advance(double days) catch (...) { return false; } } -bool DaisyController::tick() +bool DaisyBMI::tick() { if (!initialized) return false; try @@ -150,17 +150,17 @@ bool DaisyController::tick() catch (...) { return false; } } -void DaisyController::stop() +void DaisyBMI::stop() { running = false; } -bool DaisyController::is_running() const +bool DaisyBMI::is_running() const { return running; } -bool DaisyController::finalize() +bool DaisyBMI::finalize() { if (!initialized) return true; daisy().stop(); @@ -174,16 +174,16 @@ bool DaisyController::finalize() // ===== TIME ===== -double DaisyController::get_days_since_start() const { return elapsed_days_; } +double DaisyBMI::get_days_since_start() const { return elapsed_days_; } -double DaisyController::get_stop_days() const +double DaisyBMI::get_stop_days() const { if (!initialized) return std::numeric_limits::quiet_NaN(); const double h = daisy().stop_duration_hours(); return (h >= 0.0) ? h / 24.0 : std::numeric_limits::quiet_NaN(); } -std::string DaisyController::get_time_string() const +std::string DaisyBMI::get_time_string() const { if (!initialized) return ""; const Time& t = daisy().time(); @@ -193,53 +193,53 @@ std::string DaisyController::get_time_string() const return buf; } -int DaisyController::get_year() const { return initialized ? daisy().time().year() : 0; } -int DaisyController::get_month() const { return initialized ? daisy().time().month() : 0; } -int DaisyController::get_day() const { return initialized ? daisy().time().mday() : 0; } -int DaisyController::get_hour() const { return initialized ? daisy().time().hour() : 0; } -double DaisyController::get_current_dt_hours() const { return 1.0; } -double DaisyController::get_current_dt_days() const { return 1.0 / 24.0; } +int DaisyBMI::get_year() const { return initialized ? daisy().time().year() : 0; } +int DaisyBMI::get_month() const { return initialized ? daisy().time().month() : 0; } +int DaisyBMI::get_day() const { return initialized ? daisy().time().mday() : 0; } +int DaisyBMI::get_hour() const { return initialized ? daisy().time().hour() : 0; } +double DaisyBMI::get_current_dt_hours() const { return 1.0; } +double DaisyBMI::get_current_dt_days() const { return 1.0 / 24.0; } // ===== GROUNDWATER ===== -double DaisyController::get_groundwater_depth() const +double DaisyBMI::get_groundwater_depth() const { return daisy().get_groundwater_table(); } -bool DaisyController::set_groundwater_depth(double depth_cm) +bool DaisyBMI::set_groundwater_depth(double depth_cm) { daisy().set_groundwater_table(depth_cm); return true; } -auto DaisyController::estimate_sy_perturbation(double dh_cm) +auto DaisyBMI::estimate_sy_perturbation(double dh_cm) -> std::tuple, std::vector> { return daisy().estimate_sy_perturbation(dh_cm); } -Daisy& DaisyController::daisy_ref () +Daisy& DaisyBMI::daisy_ref () { return daisy (); } -double DaisyController::get_pressure_head_at_depth(double) const +double DaisyBMI::get_pressure_head_at_depth(double) const { return 0.0; // stub } // ===== SOIL WATER ===== -std::vector DaisyController::get_soil_water_content() const { return {}; } -double DaisyController::get_soil_water_at_depth(double) const { return 0.0; } -double DaisyController::get_matric_potential_at_depth(double) const { return 0.0; } -double DaisyController::get_cumulative_water_to_depth(double) const { return 0.0; } -double DaisyController::get_available_water() const { return 0.0; } +std::vector DaisyBMI::get_soil_water_content() const { return {}; } +double DaisyBMI::get_soil_water_at_depth(double) const { return 0.0; } +double DaisyBMI::get_matric_potential_at_depth(double) const { return 0.0; } +double DaisyBMI::get_cumulative_water_to_depth(double) const { return 0.0; } +double DaisyBMI::get_available_water() const { return 0.0; } // ===== RECHARGE ===== -std::vector DaisyController::get_recharge_rate() const +std::vector DaisyBMI::get_recharge_rate() const { std::vector flux = daisy().get_flux_array(); for (double& v : flux) @@ -247,19 +247,19 @@ std::vector DaisyController::get_recharge_rate() const return flux; } -double DaisyController::get_cumulative_recharge() const { return 0.0; } +double DaisyBMI::get_cumulative_recharge() const { return 0.0; } -std::vector DaisyController::get_pressure_head_array() const +std::vector DaisyBMI::get_pressure_head_array() const { return daisy().get_h_array(); // [cm], nlayers } -std::vector DaisyController::get_theta_array() const +std::vector DaisyBMI::get_theta_array() const { return daisy().get_theta_array(); // [-], nlayers } -std::vector DaisyController::get_theta_sat_array() const +std::vector DaisyBMI::get_theta_sat_array() const { return daisy().get_theta_sat_array(); // [-], nlayers, static soil param } @@ -267,58 +267,58 @@ std::vector DaisyController::get_theta_sat_array() const // ===== WATER BALANCE STUBS ===== -double DaisyController::get_et_rate() const { return 0.0; } -double DaisyController::get_cumulative_et() const { return 0.0; } -double DaisyController::get_transpiration_rate() const { return 0.0; } -double DaisyController::get_evaporation_rate() const { return 0.0; } -double DaisyController::get_drainage_rate() const { return 0.0; } -double DaisyController::get_cumulative_drainage() const { return 0.0; } -double DaisyController::get_runoff_rate() const +double DaisyBMI::get_et_rate() const { return 0.0; } +double DaisyBMI::get_cumulative_et() const { return 0.0; } +double DaisyBMI::get_transpiration_rate() const { return 0.0; } +double DaisyBMI::get_evaporation_rate() const { return 0.0; } +double DaisyBMI::get_drainage_rate() const { return 0.0; } +double DaisyBMI::get_cumulative_drainage() const { return 0.0; } +double DaisyBMI::get_runoff_rate() const { return daisy().get_runoff_rate(); } // [mm/day] -double DaisyController::get_cumulative_runoff() const { return 0.0; } +double DaisyBMI::get_cumulative_runoff() const { return 0.0; } // ===== CROP STUBS ===== -double DaisyController::get_leaf_area_index() const { return 0.0; } -double DaisyController::get_root_depth() const { return 0.0; } -double DaisyController::get_aboveground_biomass() const { return 0.0; } +double DaisyBMI::get_leaf_area_index() const { return 0.0; } +double DaisyBMI::get_root_depth() const { return 0.0; } +double DaisyBMI::get_aboveground_biomass() const { return 0.0; } // ===== SOIL STUBS ===== -double DaisyController::get_soil_temperature_at_depth(double) const { return 0.0; } -double DaisyController::get_soil_layer_thickness(int) const { return 0.0; } -double DaisyController::get_soil_layer_top(int) const { return 0.0; } -double DaisyController::get_soil_layer_bottom(int) const { return 0.0; } +double DaisyBMI::get_soil_temperature_at_depth(double) const { return 0.0; } +double DaisyBMI::get_soil_layer_thickness(int) const { return 0.0; } +double DaisyBMI::get_soil_layer_top(int) const { return 0.0; } +double DaisyBMI::get_soil_layer_bottom(int) const { return 0.0; } // ===== WEATHER STUBS ===== -double DaisyController::get_rainfall_rate() const { return 0.0; } -double DaisyController::get_cumulative_rainfall() const { return 0.0; } -double DaisyController::get_reference_et() const { return 0.0; } -double DaisyController::get_air_temperature() const { return 0.0; } -double DaisyController::get_air_humidity() const { return 0.0; } -double DaisyController::get_wind_speed() const { return 0.0; } +double DaisyBMI::get_rainfall_rate() const { return 0.0; } +double DaisyBMI::get_cumulative_rainfall() const { return 0.0; } +double DaisyBMI::get_reference_et() const { return 0.0; } +double DaisyBMI::get_air_temperature() const { return 0.0; } +double DaisyBMI::get_air_humidity() const { return 0.0; } +double DaisyBMI::get_wind_speed() const { return 0.0; } // ===== NITROGEN STUBS ===== -double DaisyController::get_mineral_nitrogen() const { return 0.0; } -double DaisyController::get_cumulative_n_leaching() const { return 0.0; } -double DaisyController::get_cumulative_n_uptake() const { return 0.0; } +double DaisyBMI::get_mineral_nitrogen() const { return 0.0; } +double DaisyBMI::get_cumulative_n_leaching() const { return 0.0; } +double DaisyBMI::get_cumulative_n_uptake() const { return 0.0; } // ===== DIAGNOSTICS STUBS ===== -std::string DaisyController::get_last_message() const { return ""; } -bool DaisyController::has_errors() const { return false; } -double DaisyController::get_simulation_time() const { return 0.0; } +std::string DaisyBMI::get_last_message() const { return ""; } +bool DaisyBMI::has_errors() const { return false; } +double DaisyBMI::get_simulation_time() const { return 0.0; } // ===== CONSTRUCTOR / DESTRUCTOR ===== -DaisyController::DaisyController() +DaisyBMI::DaisyBMI() : initialized(false), running(false) { } -DaisyController::~DaisyController() +DaisyBMI::~DaisyBMI() { if (initialized) { @@ -326,7 +326,7 @@ DaisyController::~DaisyController() } } -DaisyController::DaisyController(DaisyController&& other) noexcept +DaisyBMI::DaisyBMI(DaisyBMI&& other) noexcept : toplevel_(std::move(other.toplevel_)), initialized(other.initialized), running(other.running), @@ -337,7 +337,7 @@ DaisyController::DaisyController(DaisyController&& other) noexcept other.elapsed_days_ = 0.0; } -DaisyController& DaisyController::operator=(DaisyController&& other) noexcept +DaisyBMI& DaisyBMI::operator=(DaisyBMI&& other) noexcept { if (this != &other) { @@ -353,14 +353,14 @@ DaisyController& DaisyController::operator=(DaisyController&& other) noexcept return *this; } -double DaisyController::get_column_area() const +double DaisyBMI::get_column_area() const { return daisy().get_column_area(); } -std::vector DaisyController::get_layer_tops() const +std::vector DaisyBMI::get_layer_tops() const { return daisy().get_layer_tops(); } -std::vector DaisyController::get_layer_bottoms() const +std::vector DaisyBMI::get_layer_bottoms() const { return daisy().get_layer_bottoms(); } -int DaisyController::get_soil_layer_count() const +int DaisyBMI::get_soil_layer_count() const { return static_cast(daisy().get_layer_tops().size()); } From bba3fd14969d0fe5ec8dbfedc4995d56f4d8bf52 Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Mon, 29 Jun 2026 10:03:29 +0200 Subject: [PATCH 10/11] Remove accidental Sy perturbation from BMI branch --- include/daisy/column.h | 6 -- include/daisy/daisy.h | 3 - include/daisy/soil/soil_water.h | 2 - include/programs/bmi.h | 10 ---- include/programs/daisy_bmi.h | 9 --- src/daisy/column_std.C | 98 --------------------------------- src/daisy/daisy.C | 8 --- src/daisy/soil/soil_water.C | 7 --- src/programs/bmi.C | 8 --- src/programs/daisy_bmi.C | 6 -- 10 files changed, 157 deletions(-) diff --git a/include/daisy/column.h b/include/daisy/column.h index a54b80d51..6a97db259 100644 --- a/include/daisy/column.h +++ b/include/daisy/column.h @@ -25,7 +25,6 @@ #include "object_model/model_framed.h" #include "daisy/manager/irrigate.h" -#include #include class Frame; @@ -150,11 +149,6 @@ class Column : public ModelFramed virtual double get_column_area () const; // [cm²] virtual std::vector get_layer_tops () const; // [cm], negative downward virtual std::vector get_layer_bottoms () const; // [cm], negative downward - // Perturb GW table by dh_cm, re-run Richards, compute Sy = Σ(Δθ·Δz)/dh, - // then restore state. Default: no-op returning 0. - virtual std::tuple, std::vector> - estimate_sy_perturbation (double /*dh_cm*/) - { return {0.0, {}, {}}; } // Current development stage for the crop named "crop", or // Crop::DSremove if no such crop is present. diff --git a/include/daisy/daisy.h b/include/daisy/daisy.h index ef0d86d5b..a11cba49e 100644 --- a/include/daisy/daisy.h +++ b/include/daisy/daisy.h @@ -24,7 +24,6 @@ #define DAISY_H #include "programs/program.h" -#include #include #include @@ -70,8 +69,6 @@ class DAISY_EXPORT Daisy : public Program // Python/BMI coupling: forwarding to first (or pos-th) column. double get_groundwater_table (unsigned int pos = 0u) const; // [cm] void set_groundwater_table (double cm, unsigned int pos = 0u); - std::tuple, std::vector> - estimate_sy_perturbation (double dh_cm, unsigned int pos = 0u); /** Hours from simulation start to the configured stop time, or -1 if open-ended. */ double stop_duration_hours() const; diff --git a/include/daisy/soil/soil_water.h b/include/daisy/soil/soil_water.h index 540524e4d..65ffd522f 100644 --- a/include/daisy/soil/soil_water.h +++ b/include/daisy/soil/soil_water.h @@ -171,8 +171,6 @@ class SoilWater void set_matrix (const std::vector& h, const std::vector& Theta, const std::vector& q); - // Restore S_sum_ from a saved snapshot (used by estimate_sy_perturbation). - void restore_S_sum (const std::vector& snapshot); void set_tertiary (const std::vector& Theta_p, const std::vector& q_p, const std::vector& S_B2M, diff --git a/include/programs/bmi.h b/include/programs/bmi.h index 095b79b37..138d53beb 100644 --- a/include/programs/bmi.h +++ b/include/programs/bmi.h @@ -9,7 +9,6 @@ #define BMI_H #include -#include #include #include "programs/daisy_bmi.h" @@ -93,15 +92,6 @@ class BMI /** Set scalar input value by variable name from src[0]. */ void set_value(const std::string& name, const double* src); - /** - * Estimate specific yield for column `col` via head perturbation + Richards re-solve. - * @param dh_cm Perturbation magnitude in cm (default 1 cm) - * @param col Column index (default 0) - * @return Estimated Sy [-] - */ - std::tuple, std::vector> - estimate_sy_perturbation(double dh_cm = 1.0, unsigned int col = 0u); - // ===== CONSTRUCTOR / DESTRUCTOR ===== BMI(); ~BMI(); diff --git a/include/programs/daisy_bmi.h b/include/programs/daisy_bmi.h index 46a8b9d65..a87df6dbf 100644 --- a/include/programs/daisy_bmi.h +++ b/include/programs/daisy_bmi.h @@ -12,7 +12,6 @@ #include #include -#include #include // Forward declarations @@ -180,14 +179,6 @@ class DaisyBMI */ bool set_groundwater_depth(double depth_cm); - /** - * Estimate specific yield via GW head perturbation + Richards re-solve. - * @param dh_cm Head perturbation in cm (default 1 cm) - * @return Estimated Sy [-] - */ - std::tuple, std::vector> - estimate_sy_perturbation(double dh_cm = 1.0); - /** * Get pressure head at specific depth * @param z_cm Depth in cm (positive downward) diff --git a/src/daisy/column_std.C b/src/daisy/column_std.C index de150aa74..07b6645e9 100644 --- a/src/daisy/column_std.C +++ b/src/daisy/column_std.C @@ -95,15 +95,6 @@ struct ColumnStandard : public Column std::vector theta_sat_cached_; // [-], saturated θ (static soil param) mutable double runoff_rate_cached_ = 0.0; // [mm], accumulated since last read - // Snapshots taken just before Richards solve (needed by estimate_sy_perturbation). - Time time_last_; - const Weather* weather_last_ = nullptr; - const Scope* scope_last_ = nullptr; - std::vector h_pretick_; - std::vector theta_pretick_; - double table_pretick_ = 0.0; - std::vector S_sum_pretick_; - // Log variables. double yield_DM; double yield_N; @@ -213,8 +204,6 @@ public: std::vector get_theta_array () const override; std::vector get_theta_sat_array () const override; double get_runoff_rate () const override; - std::tuple, std::vector> - estimate_sy_perturbation (double dh_cm) override; // Simulation. void clear (); @@ -881,21 +870,6 @@ ColumnStandard::tick_move (const Metalib& metalib, surface->temperature (), dt, msg); soil_water->reset_old (); // Set Theta_old to Theta here. - // BMI: snapshot h + theta + S_sum + GW table just before Richards solve. - { - const size_t _n = geometry.cell_size (); - h_pretick_.resize (_n); - theta_pretick_.resize (_n); - for (size_t _i = 0; _i < _n; ++_i) - { - h_pretick_[_i] = soil_water->h (_i); - theta_pretick_[_i] = soil_water->Theta (_i); - } - table_pretick_ = groundwater->table (); - S_sum_pretick_.resize (_n); - for (size_t _i = 0; _i < _n; ++_i) - S_sum_pretick_[_i] = soil_water->S_sum (_i); - } chemistry->mass_balance (geometry, *soil_water); soil_water->tick_ice (geometry, *soil, dt, msg); movement->tick (*soil, *soil_water, *soil_heat, @@ -943,9 +917,6 @@ ColumnStandard::tick_move (const Metalib& metalib, // BMI: cache post-Richards arrays and runoff for BMI getters. { const size_t n = geometry.cell_size (); - time_last_ = time_end; - weather_last_ = weather.get () ? weather.get () : global_weather; - scope_last_ = extern_scope; bottom_flux_cached_ = soil_water->q_matrix (n); @@ -1502,72 +1473,3 @@ ColumnStandard::get_runoff_rate () const runoff_rate_cached_ = 0.0; return v; } - -auto ColumnStandard::estimate_sy_perturbation (double dh_cm) - -> std::tuple, std::vector> -{ - const std::tuple, std::vector> failed - = {0.0, {}, {}}; - - if (!weather_last_ || !scope_last_ || h_pretick_.empty () || S_sum_pretick_.empty ()) - return failed; - - const double dh_safe = std::min (dh_cm, -table_pretick_); - if (dh_safe <= 0.0) - return failed; - - const size_t n = geometry.cell_size (); - - auto restore = [&]() - { - for (size_t i = 0; i < n; ++i) - soil_water->set_content (i, h_pretick_[i], theta_pretick_[i]); - soil_water->restore_S_sum (S_sum_pretick_); - }; - - struct TickResult { std::vector theta, flux_mm_d, h_cm; }; - auto run_tick = [&](double gw_table) -> TickResult - { - groundwater->set_table (gw_table); - soil_water->reset_old (); - try - { - movement->tick (*soil, *soil_water, *soil_heat, - *surface, *groundwater, - time_last_, *scope_last_, *weather_last_, - 24.0, Treelog::null ()); - } - catch (...) {} - TickResult r; - r.theta.resize (n); r.flux_mm_d.resize (n); r.h_cm.resize (n); - for (size_t i = 0; i < n; ++i) - { - r.theta[i] = soil_water->Theta (i); - r.h_cm[i] = soil_water->h (i); - r.flux_mm_d[i] = soil_water->q_matrix (i + 1) * 10.0 * 24.0; - } - return r; - }; - - restore (); - const TickResult A = run_tick (table_pretick_); - - restore (); - const TickResult B = run_tick (table_pretick_ + dh_safe); - - double numerator = 0.0; - for (size_t i = 0; i < n; ++i) - { - const double dtheta = B.theta[i] - A.theta[i]; - const double dz_m = (geometry.cell_top (i) - geometry.cell_bottom (i)) * 1e-2; - numerator += dtheta * dz_m; - } - const double sy = std::abs (numerator) / std::abs (dh_safe * 1e-2); - - groundwater->set_table (table_pretick_); - for (size_t i = 0; i < n; ++i) - soil_water->set_content (i, h_array_cached_[i], theta_array_cached_[i]); - soil_water->restore_S_sum (std::vector (n, 0.0)); - - return {sy, B.flux_mm_d, B.h_cm}; -} diff --git a/src/daisy/daisy.C b/src/daisy/daisy.C index 14a046c56..cf26a32ba 100644 --- a/src/daisy/daisy.C +++ b/src/daisy/daisy.C @@ -584,14 +584,6 @@ Daisy::set_groundwater_table (double cm, unsigned int pos) if (col) col->set_groundwater_table (cm); } -auto Daisy::estimate_sy_perturbation (double dh_cm, unsigned int pos) - -> std::tuple, std::vector> -{ - Column* col = impl->field->find (pos); - return col ? col->estimate_sy_perturbation (dh_cm) - : std::make_tuple (0.0, std::vector{}, std::vector{}); -} - double Daisy::stop_duration_hours() const { diff --git a/src/daisy/soil/soil_water.C b/src/daisy/soil/soil_water.C index 64736bd2b..ad0f6c5d0 100644 --- a/src/daisy/soil/soil_water.C +++ b/src/daisy/soil/soil_water.C @@ -239,13 +239,6 @@ SoilWater::set_matrix (const std::vector& h, q_matrix_ = q; } -void -SoilWater::restore_S_sum (const std::vector& snapshot) -{ - daisy_assert (S_sum_.size () == snapshot.size ()); - S_sum_ = snapshot; -} - void SoilWater::set_tertiary (const std::vector& Theta_p, const std::vector& q_p, diff --git a/src/programs/bmi.C b/src/programs/bmi.C index 15039452f..f8cd8000c 100644 --- a/src/programs/bmi.C +++ b/src/programs/bmi.C @@ -315,11 +315,3 @@ void BMI::set_value(const std::string& name, const double* src) throw std::invalid_argument("Unknown input variable: " + name); } -auto DaisyBMI::estimate_sy_perturbation(double dh_cm, unsigned int col) - -> std::tuple, std::vector> -{ - if (col != 0u) - throw std::invalid_argument("estimate_sy_perturbation: col > 0 not yet supported via DaisyBMI"); - return ctrl_.estimate_sy_perturbation(dh_cm); -} - diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C index 6693497fd..c823577e0 100644 --- a/src/programs/daisy_bmi.C +++ b/src/programs/daisy_bmi.C @@ -213,12 +213,6 @@ bool DaisyBMI::set_groundwater_depth(double depth_cm) return true; } -auto DaisyBMI::estimate_sy_perturbation(double dh_cm) - -> std::tuple, std::vector> -{ - return daisy().estimate_sy_perturbation(dh_cm); -} - Daisy& DaisyBMI::daisy_ref () { return daisy (); From 32edb7023fd03a3942b87c08c724047aa395172f Mon Sep 17 00:00:00 2001 From: HendrikKok Date: Mon, 29 Jun 2026 10:51:19 +0200 Subject: [PATCH 11/11] cleanup --- python/examples/daisy_bmi_example.py | 22 +++++++++++----------- python/src/daisy/__init__.py | 4 ++-- src/programs/api_bindings.cpp | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/python/examples/daisy_bmi_example.py b/python/examples/daisy_bmi_example.py index 912b92a30..2ff3d1659 100644 --- a/python/examples/daisy_bmi_example.py +++ b/python/examples/daisy_bmi_example.py @@ -10,7 +10,7 @@ import os import numpy as np -from daisy import API +from daisy import BMI from pathlib import Path # daisy config file @@ -38,34 +38,34 @@ def compute_zh0(tops: np.ndarray, h_daisy: np.ndarray) -> float | None: # -- Daisy simulation -- # initialise Daisy -api = API() +daisy_bmi = BMI() os.chdir(config_file.parent) -api.initialize(str(config_file)) +daisy_bmi.initialize(str(config_file)) # get static internals -area = api.get_value("column__area") # scalar cm² -tops = api.get_value_array("soil_layer__top_depth") # np.NDarray[float] cm -bots = api.get_value_array("soil_layer__bottom_depth") # np.NDarray[float] cm +area = daisy_bmi.get_value("column__area") # scalar cm² +tops = daisy_bmi.get_value_array("soil_layer__top_depth") # np.NDarray[float] cm +bots = daisy_bmi.get_value_array("soil_layer__bottom_depth") # np.NDarray[float] cm n_lay = len(tops) # runs simulations per daily time step end_time = 365 -print(f"BMI time units are: {api.get_time_units()}") +print(f"BMI time units are: {daisy_bmi.get_time_units()}") for itime in np.arange(end_time-1): # Get current time - t = api.get_current_time() + t = daisy_bmi.get_current_time() # Advance Daisy one day - api.update_until(t + 1.0) + daisy_bmi.update_until(t + 1.0) # get pressure heads and GW table depth - h = api.get_value_array("soil_water__pressure_head") + h = daisy_bmi.get_value_array("soil_water__pressure_head") gwl = compute_zh0(tops, h_daisy=h) if gwl is None: print(f'gwl at day {t:.1f} is below Daisy column') else: print(f'gwl at day {t:.1f} is {gwl}') -api.finalize() +daisy_bmi.finalize() diff --git a/python/src/daisy/__init__.py b/python/src/daisy/__init__.py index a1740f3ed..8380bbadd 100644 --- a/python/src/daisy/__init__.py +++ b/python/src/daisy/__init__.py @@ -9,11 +9,11 @@ _os.add_dll_directory(_pkg_dir) try: - from .daisy_bmi import BMI, API + from .daisy_bmi import BMI except ImportError as e: raise ImportError( "Failed to import daisy_bmi. Make sure the .pyd and required DLLs are present.\n" f" .pyd directory: {_pkg_dir}" ) from e -__all__ = ["BMI", "API"] +__all__ = ["BMI"] diff --git a/src/programs/api_bindings.cpp b/src/programs/api_bindings.cpp index 0ba77e5fe..bc526f327 100644 --- a/src/programs/api_bindings.cpp +++ b/src/programs/api_bindings.cpp @@ -60,7 +60,7 @@ PYBIND11_MODULE(daisy_bmi, m) bmi.finalize() )pbdoc"; - py::class_(m, "DaisyBMI") + py::class_(m, "BMI") .def(py::init<>(), "Create a new Daisy BMI instance")