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 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/.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/CMakeLists.txt b/CMakeLists.txt index 0323f6c93..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) @@ -94,6 +95,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/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/include/daisy/column.h b/include/daisy/column.h index 0a458e986..6a97db259 100644 --- a/include/daisy/column.h +++ b/include/daisy/column.h @@ -136,6 +136,20 @@ 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 + // 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..a11cba49e 100644 --- a/include/daisy/daisy.h +++ b/include/daisy/daisy.h @@ -37,7 +37,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 +65,23 @@ 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); + /** 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 +90,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/programs/bmi.h b/include/programs/bmi.h new file mode 100644 index 000000000..138d53beb --- /dev/null +++ b/include/programs/bmi.h @@ -0,0 +1,116 @@ +// 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 "programs/daisy_bmi.h" + +class Daisy; // forward declaration for protected daisy() accessor + +/** + * @class BMI + * @brief BMI 2.0 wrapper around DaisyBMI + * + * 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 BMI +{ +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); + + // ===== CONSTRUCTOR / DESTRUCTOR ===== + BMI(); + ~BMI(); + +private: + DaisyBMI 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; + +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 new file mode 100644 index 000000000..a87df6dbf --- /dev/null +++ b/include/programs/daisy_bmi.h @@ -0,0 +1,473 @@ +// 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 + +// Forward declarations +class Toplevel; +class Daisy; + +/** + * @class DaisyBMI + * @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 DaisyBMI +{ +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); + + /** + * 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 + */ + DaisyBMI(); + + /** + * Destructor + */ + ~DaisyBMI(); + + // Delete copy operations + DaisyBMI(const DaisyBMI&) = delete; + DaisyBMI& operator=(const DaisyBMI&) = delete; + + // Allow move operations + 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 DaisyBMI. */ + Daisy& daisy_ref (); +}; + +#endif // DAISY_BMI_H diff --git a/python/examples/daisy_bmi_example.py b/python/examples/daisy_bmi_example.py new file mode 100644 index 000000000..2ff3d1659 --- /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 BMI +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_bmi = BMI() +os.chdir(config_file.parent) +daisy_bmi.initialize(str(config_file)) + +# get static internals +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: {daisy_bmi.get_time_units()}") + +for itime in np.arange(end_time-1): + # Get current time + t = daisy_bmi.get_current_time() + + # Advance Daisy one day + daisy_bmi.update_until(t + 1.0) + + # get pressure heads and GW table depth + 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}') + +daisy_bmi.finalize() diff --git a/python/daisy_py_fun_example.py b/python/examples/daisy_py_fun_example.py similarity index 95% rename from python/daisy_py_fun_example.py rename to python/examples/daisy_py_fun_example.py index a6837d526..d69604731 100644 --- a/python/daisy_py_fun_example.py +++ b/python/examples/daisy_py_fun_example.py @@ -1,7 +1,8 @@ '''Module illustrating how to implement a concrete DaisyPyFun that can be called from Daisy''' import math -from daisy_py_fun_interface import DaisyPyFun +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. @@ -9,7 +10,6 @@ class Tmin(DaisyPyFun): range = 'None' domain = 'dg C' - def __init__(self, args): '''Initialize T_min. @@ -58,4 +58,4 @@ def _apply(self, t): 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 + 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..8380bbadd --- /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 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"] 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/daisy_py_fun_interface.py b/python/src/daisy/plugins/daisy_py_fun_interface.py similarity index 100% rename from python/daisy_py_fun_interface.py rename to python/src/daisy/plugins/daisy_py_fun_interface.py 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/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..07b6645e9 100644 --- a/src/daisy/column_std.C +++ b/src/daisy/column_std.C @@ -87,6 +87,14 @@ 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 + // Log variables. double yield_DM; double yield_N; @@ -184,6 +192,19 @@ 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; + // Simulation. void clear (); void tick_source (const Scope&, const Time&, Treelog&); @@ -848,6 +869,7 @@ 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. + chemistry->mass_balance (geometry, *soil_water); soil_water->tick_ice (geometry, *soil, dt, msg); movement->tick (*soil, *soil_water, *soil_heat, @@ -891,6 +913,24 @@ 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 (); + + 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 +1251,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 +1411,65 @@ 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; +} diff --git a/src/daisy/daisy.C b/src/daisy/daisy.C index 0fd0603a9..cf26a32ba 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,99 @@ 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); +} + +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/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..d89ad9bf1 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 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 +) + + +# ===== 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() diff --git a/src/programs/api_bindings.cpp b/src/programs/api_bindings.cpp new file mode 100644 index 000000000..bc526f327 --- /dev/null +++ b/src/programs/api_bindings.cpp @@ -0,0 +1,138 @@ +// api_bindings.cpp -- pybind11 bindings for BMI and DaisyAPI. +// +// 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 +// $(python3 -m pybind11 --includes) +// daisy_api_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.BMI() + 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, "BMI") + + .def(py::init<>(), "Create a new Daisy BMI instance") + + // --- Lifecycle --- + .def("initialize", &BMI::initialize, py::arg("config_file"), + "Initialize from a .dai config file") + .def("update", &BMI::update, + "Advance model by one timestep") + .def("update_until", &BMI::update_until, py::arg("time"), + "Advance model until time (days since start)") + .def("finalize", &BMI::finalize, + "Finalize and release all resources") + + // --- Information --- + .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", &BMI::get_start_time, + "Simulation start time [days]") + .def("get_end_time", &BMI::get_end_time, + "Simulation end time [days] (NaN = open-ended)") + .def("get_current_time", &BMI::get_current_time, + "Current simulation time [days since start]") + .def("get_time_step", &BMI::get_time_step, + "Current timestep size [days]") + .def("get_time_units", &BMI::get_time_units, + "Time unit string, e.g. 'd'") + + // --- Variable metadata --- + .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", &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 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 + // 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", + [](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 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()); + return arr; + }, + py::arg("name"), + "Get an array output value as a numpy array of doubles"); +} diff --git a/src/programs/bmi.C b/src/programs/bmi.C new file mode 100644 index 000000000..f8cd8000c --- /dev/null +++ b/src/programs/bmi.C @@ -0,0 +1,317 @@ +// daisy_bmi.C -- BMI 2.0 wrapper implementation for Daisy +// +// This file is part of Daisy. + +#include "programs/bmi.h" +#include "daisy/daisy.h" +#include +#include +#include + +// ===== STATIC VARIABLE LISTS ===== + +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 BMI::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 ===== + +BMI::BMI() + : start_time_days_(0.0) + , current_time_days_(0.0) + , dt_days_(1.0 / 24.0) // default: hourly timestep +{ +} + +BMI::~BMI() +{ + // DaisyBMI destructor handles cleanup +} + +// ===== LIFECYCLE ===== + +void BMI::initialize(const std::string& config_file) +{ + try + { + ctrl_.initialize(config_file); + } + catch (const std::exception& e) + { + throw std::runtime_error( + std::string("BMI::initialize failed for config: ") + config_file + + "\n Reason: " + e.what()); + } + catch (const int code) + { + throw std::runtime_error( + 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("BMI::initialize failed for config: ") + config_file + + "\n Reason: " + msg); + } + catch (...) + { + throw std::runtime_error( + std::string("BMI::initialize failed for config: ") + config_file + + "\n Reason: unknown exception"); + } + + if (!ctrl_.is_initialized()) + throw std::runtime_error( + 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 + // 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 BMI::update() +{ + if (!ctrl_.tick()) + 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(); + + // 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 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). + 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 BMI::finalize() +{ + ctrl_.finalize(); +} + +// ===== INFORMATION ===== + +std::string BMI::get_component_name() const +{ + return "Daisy"; +} + +std::vector BMI::get_input_var_names() const +{ + return INPUT_VARS; +} + +std::vector BMI::get_output_var_names() const +{ + return OUTPUT_VARS; +} + +// ===== TIME ===== + +double BMI::get_start_time() const +{ + return start_time_days_; +} + +double BMI::get_end_time() const +{ + return ctrl_.get_stop_days(); +} + +double BMI::get_current_time() const +{ + return current_time_days_; +} + +double BMI::get_time_step() const +{ + return dt_days_; +} + +std::string BMI::get_time_units() const +{ + return "d"; +} + +// ===== VARIABLE INFO ===== + +std::string BMI::get_var_type(const std::string& /*name*/) const +{ + return "double"; +} + +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"; + 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 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" + || name == "soil_water__theta" + || name == "soil__theta_sat") + return ctrl_.get_soil_layer_count() * sizeof(double); + return sizeof(double); +} + +int BMI::get_var_itemsize(const std::string& /*name*/) const +{ + return static_cast(sizeof(double)); +} + +std::string BMI::get_var_location(const std::string& /*name*/) const +{ + return "none"; // scalars have no grid location +} + +int BMI::get_var_grid(const std::string& /*name*/) const +{ + return 0; // all vars on the single-cell scalar grid +} + +// ===== GRID INFO ===== + +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 BMI::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 BMI::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 BMI::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 DaisyBMI + // ctrl_.set_irrigation_rate(src[0]); + throw std::runtime_error("irrigation__rate setter not yet implemented in DaisyBMI"); + } + throw std::invalid_argument("Unknown input variable: " + name); +} + diff --git a/python/call_python_function.cpp b/src/programs/call_python_function.cpp similarity index 100% rename from python/call_python_function.cpp rename to src/programs/call_python_function.cpp diff --git a/src/programs/daisy_bmi.C b/src/programs/daisy_bmi.C new file mode 100644 index 000000000..c823577e0 --- /dev/null +++ b/src/programs/daisy_bmi.C @@ -0,0 +1,360 @@ +// 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& DaisyBMI::daisy() const +{ + return dynamic_cast(toplevel_->program()); +} + +void DaisyBMI::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 DaisyBMI::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 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. + load_config_file(config_file); + start(); + return initialized; +} + +// ===== LIFECYCLE ===== + +void DaisyBMI::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 DaisyBMI::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 DaisyBMI::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 DaisyBMI::stop() +{ + running = false; +} + +bool DaisyBMI::is_running() const +{ + return running; +} + +bool DaisyBMI::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 DaisyBMI::get_days_since_start() const { return elapsed_days_; } + +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 DaisyBMI::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 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 DaisyBMI::get_groundwater_depth() const +{ + return daisy().get_groundwater_table(); +} + +bool DaisyBMI::set_groundwater_depth(double depth_cm) +{ + daisy().set_groundwater_table(depth_cm); + return true; +} + +Daisy& DaisyBMI::daisy_ref () +{ + return daisy (); +} + +double DaisyBMI::get_pressure_head_at_depth(double) const +{ + return 0.0; // stub +} + +// ===== SOIL WATER ===== + +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 DaisyBMI::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 DaisyBMI::get_cumulative_recharge() const { return 0.0; } + +std::vector DaisyBMI::get_pressure_head_array() const +{ + return daisy().get_h_array(); // [cm], nlayers +} + +std::vector DaisyBMI::get_theta_array() const +{ + return daisy().get_theta_array(); // [-], nlayers +} + +std::vector DaisyBMI::get_theta_sat_array() const +{ + return daisy().get_theta_sat_array(); // [-], nlayers, static soil param +} + + +// ===== WATER BALANCE STUBS ===== + +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 DaisyBMI::get_cumulative_runoff() const { return 0.0; } + +// ===== CROP STUBS ===== + +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 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 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 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 DaisyBMI::get_last_message() const { return ""; } +bool DaisyBMI::has_errors() const { return false; } +double DaisyBMI::get_simulation_time() const { return 0.0; } + +// ===== CONSTRUCTOR / DESTRUCTOR ===== + +DaisyBMI::DaisyBMI() + : initialized(false), running(false) +{ +} + +DaisyBMI::~DaisyBMI() +{ + if (initialized) + { + try { finalize(); } catch (...) {} + } +} + +DaisyBMI::DaisyBMI(DaisyBMI&& 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; +} + +DaisyBMI& DaisyBMI::operator=(DaisyBMI&& 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 DaisyBMI::get_column_area() const +{ return daisy().get_column_area(); } + +std::vector DaisyBMI::get_layer_tops() const +{ return daisy().get_layer_tops(); } + +std::vector DaisyBMI::get_layer_bottoms() const +{ return daisy().get_layer_bottoms(); } + +int DaisyBMI::get_soil_layer_count() const +{ return static_cast(daisy().get_layer_tops().size()); } 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/python/CMakeLists.txt b/tools/call_python/CMakeLists.txt similarity index 100% rename from python/CMakeLists.txt rename to tools/call_python/CMakeLists.txt 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; +}