diff --git a/CMakeLists.txt b/CMakeLists.txt index d237d62..53cf0cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,16 +16,35 @@ MESSAGE(STATUS "GSL found") MESSAGE(STATUS " - version ${GSL_VERSION}") MESSAGE(STATUS " - included files located in ${GSL_INCLUDE_DIRS}") -add_executable(${PROJECT_NAME} src/troll.cpp src/lookUpTables_cache.cpp) +add_executable(${PROJECT_NAME} src/troll.cpp) + +target_compile_options(${PROJECT_NAME} PRIVATE -O3 -Wall) +target_link_libraries(${PROJECT_NAME} PRIVATE GSL::gsl GSL::gslcblas) + +# Add the LUT generator executable +add_executable(generate_lut_header src/generate-lut-header.cpp) + +# Include the 'include' directory for the generator +target_include_directories(generate_lut_header PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include) + +# Add a custom lut target to generate lut.hpp in the include directory +add_custom_target(lut + COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_CURRENT_SOURCE_DIR}/include + COMMAND generate_lut_header + BYPRODUCTS ${CMAKE_CURRENT_SOURCE_DIR}/include/lut.hpp + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include + DEPENDS generate_lut_header + COMMENT "Generating LUT header (include/lut.hpp)..." +) + +# Make the main TROLL target depend on the lut target +add_dependencies(${PROJECT_NAME} lut) # Include directories target_include_directories(${PROJECT_NAME} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/include # Public headers ) -target_compile_options(${PROJECT_NAME} PRIVATE -O3 -Wall) -target_link_libraries(${PROJECT_NAME} PRIVATE GSL::gsl GSL::gslcblas) - # sanity check add_custom_target(sanity COMMAND ${CMAKE_SOURCE_DIR}/sanity/check.sh diff --git a/include/troll.hpp b/include/troll.hpp index 4981077..bf76805 100644 --- a/include/troll.hpp +++ b/include/troll.hpp @@ -218,8 +218,6 @@ vector DailyMeanWindSpeed; //!< Global vector: win // Complex temperature-dependent functions used in the Farquhar model are computed once at 'Taccuracy' resolution. Leaf temperature must be comprised between 0°C and 60°C, Values are stored every 0.5°C step in Tleaf, so 120 values in total int nbTbins; //!< Global variable: number of bins for the temperature lookup tables float iTaccuracy; //!< Global variable: inverse of accuracy of a temperature bin (e.g. if Taccuracy is 0.1 or 0.5 °C, then iTaccuracy is 10.0 or 2.0, respectively) -float *LookUp_KmT(0); //!< Global vector: lookup table for Km(T) in Farquhar model -float *LookUp_GammaT(0); //!< Global vector: lookup table for Gamma(T) in Farquhar model float *LookUp_VcmaxT(0); //!< Global vector: lookup table for Vcmax(T) in Farquhar model float *LookUp_JmaxT(0); //!< Global vector: lookup table for Jmax(T) in Farquhar model //float *LookUp_Rday(0); //!< Global vector: lookup table for Rday(T) in Farquhar model //new IM: no redundancy anymore between LookUp_Rday and LookUp_Rnight diff --git a/src/generate-lut-header.cpp b/src/generate-lut-header.cpp new file mode 100644 index 0000000..8477a10 --- /dev/null +++ b/src/generate-lut-header.cpp @@ -0,0 +1,134 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include "troll_defines.hpp" + +/** + * @file generate-lut-header.cpp + * @brief Generates lut.hpp with precomputed lookup tables. + * This is a build-time tool; not part of the runtime code. + */ + +/** + * @brief Number of bins for temperature lookup tables. + * Determines the resolution of the LUTs. + */ +constexpr int n_temperature_bins = 500; + +/** + * @brief Temperature accuracy for lookup tables. + * Defines the temperature step size (in °C) between bins. + */ +constexpr float temperature_accuracy = 0.1f; + +/** + * @brief Generates a lookup table and writes it to the output file. + * @tparam N Size of the lookup table. + * @param out Output file stream. + * @param name Name of the lookup table. + * @param n_bins Size of the lookup table. + * @param compute Function to compute the value for a given temperature. + */ +void generate_lut( + std::ofstream& out, + const std::string& name, + const std::function& compute, + size_t n_bins +) { + out << " constexpr std::array " << name << " = {\n"; + for (int i = 0; i < n_bins; ++i) { + float T = i * temperature_accuracy; + /// setprecision to 10 prints all significant digits of a 32-bits float without truncation + out << std::scientific << std::setprecision(10) << compute(T); + if (i < n_bins - 1) out << ", "; + if (i % 10 == 0) out << "\n"; + } + out << "\n };\n\n"; +} + +/** + * @brief Computes the effective Michaelis-Menten constant for CO₂, \( K_m(T) \). + * + * Taken from von Caemmerer 2000, as in Domingues et al. 2010 (), for consistency. + * + * \( K_m(T) = K_c(T) \cdot \left(1 + \frac{O_2}{K_o(T)}\right) \), + * where \( K_c(T) \) and \( K_o(T) \) are temperature-dependent constants. + * + * @param T Temperature in °C. + * @return \( K_m \) in Pa. + * + * + * Km <- Kc_Tleaf*(1+O2/Ko_Tleaf) + * Kc_Tleaf <- Kc * exp((Tleaf+273.15-298.15)*delta_Kc/(298.15*8.314*(Tleaf+273.15))) + * delta_Kc <- 59356 # activation energy (kJ mol−1) - Badger & Collatz (1977) + * Ko_Tleaf <- Ko * exp((Tleaf+273.15-298.15)*delta_Ko/(298.15*8.314*(Tleaf+273.15))) + * Kc <- 40.49 # Michaelis-Menten constant for CO2 (Pa) - von Caemmerer et al. (1994) + * Ko <- 24.8 # Michaelis-Menten constant for CO2 (kPa)- von Caemmerer et al. (1994) + * O2 <- 21 # Estimated Oxygen concentration at mesophyll - (kPa) + * delta_Ko <- 35948 # activation energy (kJ mol−1) - Badger & Collatz (1977) + * + * @param T Temperature in °C. + * @return \( K_m \) value at temperature T. + */ +float compute_KmT(float T) { + +#ifdef WATER + return 404.0 * exp(((T - 25.0) / (298 * 0.00831 * (273 + T))) * 59.36) * + (1 + 210 * 1.0 / 248.0 * exp(-(T - 25.0) / (298 * 0.00831 * (273 + T)) * 35.94)); +#else + return 404.0 * exp(((T - 25.0) / (298 * 0.00831 * (273 + T))) * 59.36) * + (1 + 210 * 1.0 / 248.0 * exp(-(T - 25.0) / (298 * 0.00831 * (273 + T)) * 35.94)) * iCair; +#endif +} + +/** + * @brief Computes the \( Gamma(T) \) lookup table. + * + * Taken from von Caemmerer 2000, as in Domingues et al. 2010 (), for consistency. + * + * @param T Temperature in °C. + * @return Gamma value at temperature T. + */ +float compute_GammaT(float T) { + +#ifdef WATER + return 37.0 * exp(((T - 25.0) / (298 * 0.00831 * (273 + T))) * 23.4); +#else + return 37.0 * exp(((T - 25.0) / (298 * 0.00831 * (273 + T))) * 23.4) * iCair; +#endif +} + +int main() { + std::ofstream out("lut.hpp"); + if (!out) { + std::cerr << "Failed to open lut.hpp for writing.\n"; + return 1; + } + + out << "/// Auto-generated LUT header. Do not edit! Refer to src/generate-lut-header.cpp for details. \n"; + out << "#ifndef LUT_HPP \n"; + out << "#define LUT_HPP \n"; + out << "#include \n\n"; + out << "namespace lut {\n"; + + // Map LUT names to their computation functions + std::map> lut_functions = { + {"LookUp_KmT", compute_KmT}, + {"LookUp_GammaT", compute_GammaT} + // Add other LUTs here + }; + + // Generate all LUTs + for (const auto& [name, compute] : lut_functions) { + generate_lut(out, name, compute, n_temperature_bins); + } + + out << "} /// namespace lut \n"; + out << "#endif /// LUT_HPP \n"; + return 0; +} diff --git a/src/troll.cpp b/src/troll.cpp index 45d83fd..7124276 100644 --- a/src/troll.cpp +++ b/src/troll.cpp @@ -1,6 +1,6 @@ #include "troll.hpp" #include "constants.hpp" -#include "lookUpTables_cache.hpp" +#include "lut.hpp" // ############################################# // Species constructor @@ -1634,8 +1634,8 @@ leafFluxes Tree::Photosyn(float PPFD, float TLEAF, float CS, float DS) // Parameters for the Farquhar model int convT = int(iTaccuracy * TLEAF); - float KmT = LookUp_KmT[convT]; // with temperature dependencies - float GammaT = LookUp_GammaT[convT]; // with temperature dependencies + float KmT = lut::LookUp_KmT[convT]; // with temperature dependencies + float GammaT = lut::LookUp_GammaT[convT]; // with temperature dependencies float Rday = t_Rdark * LookUp_Rleaf[convT] * DAYRESP; // leaf respiration with temperature dependencies and inhibition due to light (DAYRESP). ////////////// Model of stomatal conductance for CO2 ///////////// @@ -1720,8 +1720,8 @@ leafFluxes Tree::Photosyn(float PPFD, float TLEAF, float CS, float DS) // Parameters for the Farquhar model int convT = int(iTaccuracy * TLEAF); - float KmT = LookUp_KmT[convT]; // with temperature dependencies - float GammaT = LookUp_GammaT[convT]; // with temperature dependencies + float KmT = lut::LookUp_KmT[convT]; // with temperature dependencies + float GammaT = lut::LookUp_GammaT[convT]; // with temperature dependencies float Rday = t_Rdark * LookUp_Rleaf[convT] * DAYRESP; // leaf respiration with temperature dependencies and inhibition due to light (DAYRESP). // Model of stomatal conductance for CO2 @@ -1933,8 +1933,8 @@ float Tree::GPPleaf(float PPFD, float VPD, float T) int convT = int(iTaccuracy * T); // temperature data at a resolution of Taccuracy=0.1°C -- stored in lookup tables ranging from 0°C to 50°C --- // if(convT>500 || isnan(convT) || convT <0) cout << t_site << " | convT: " << convT << " | T: " << T << " | PPFD: " << PPFD << " | VPD: " << VPD << endl; - float KmT = LookUp_KmT[convT]; - float GammaT = LookUp_GammaT[convT]; + float KmT = lut::LookUp_KmT[convT]; + float GammaT = lut::LookUp_GammaT[convT]; // float g1 = -3.97 * t_wsg + 6.53 (Lin et al. 2015) @@ -5554,10 +5554,6 @@ void InitialiseLookUpTables() iTaccuracy = 1.0 / Taccuracy; cout << endl << "Built-in maximal temperature: " << float(nbTbins) * Taccuracy << endl; - if (NULL == (LookUp_KmT = new float[nbTbins])) - cerr << "!!! Mem_Alloc LookUp_KmT" << endl; - if (NULL == (LookUp_GammaT = new float[nbTbins])) - cerr << "!!! Mem_Alloc LookUp_GammaT" << endl; if (NULL == (LookUp_VcmaxT = new float[nbTbins])) cerr << "!!! Mem_Alloc LookUp_VcmaxT" << endl; if (NULL == (LookUp_JmaxT = new float[nbTbins])) @@ -5592,15 +5588,6 @@ void InitialiseLookUpTables() float temper = float(i) * Taccuracy; // !!!UPDATE provide references for these equations -#ifdef WATER - LookUp_KmT[i] = 404.0 * exp(((temper - 25.0) / (298 * 0.00831 * (273 + temper))) * 59.36) * - (1 + 210 * 1.0 / 248.0 * exp(-(temper - 25.0) / (298 * 0.00831 * (273 + temper)) * 35.94)); // taken from von Caemmerer 2000, as in Domingues et al. 2010, for consistency - LookUp_GammaT[i] = 37.0 * exp(((temper - 25.0) / (298 * 0.00831 * (273 + temper))) * 23.4); // taken from von Caemmerer 2000, as in Domingues et al. 2010, for consistency -#else - LookUp_KmT[i] = 404.0 * exp(((temper - 25.0) / (298 * 0.00831 * (273 + temper))) * 59.36) * - (1 + 210 * 1.0 / 248.0 * exp(-(temper - 25.0) / (298 * 0.00831 * (273 + temper)) * 35.94)) * iCair; // taken from von Caemmerer 2000, as in Domingues et al. 2010, for consistency - LookUp_GammaT[i] = 37.0 * exp(((temper - 25.0) / (298 * 0.00831 * (273 + temper))) * 23.4) * iCair; // taken from von Caemmerer 2000, as in Domingues et al. 2010, for consistency -#endif LookUp_VcmaxT[i] = exp(26.35 - 65.33 / (0.00831 * (temper + 273.15))); // taken from Bernacchi et al. 2003 PCE, as in Domingues et al. 2010 for consistency LookUp_JmaxT[i] = exp(17.57 - 43.54 / (0.00831 * (temper + 273.15))); // taken from Bernacchi et al. 2003 PCE, as in Domingues et al. 2010 for consistency // LookUp_Rday[i]=exp((temper-25.0)*0.1*log(3.09-0.0215*(25.0+temper))); //newIM: no redundancy anymore between LookUp_Rday and LookUP_Rnight @@ -10866,7 +10853,6 @@ void FreeMem() #endif delete[] LookUp_T; - delete[] LookUp_KmT; delete[] LookUp_VPD; delete[] LookUp_flux_absorption; delete[] LookUp_flux; @@ -10875,7 +10861,6 @@ void FreeMem() delete[] LookUp_JmaxT; delete[] LookUp_Rstem; delete[] LookUp_Rleaf; // newIM: no redundancy anymore between LookUp_Rday and LookUp_Rleaf - delete[] LookUp_GammaT; // delete [] LookUp_Rnight; //newIM: no redundancy anymore between LookUp_Rday and LookUp_Rleaf delete[] LookUp_VcmaxT; #ifdef WATER