19 cache look up table

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

include/troll.hpp

@@ -218,8 +218,6 @@ vector<float> 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

src/generate-lut-header.cpp

@@ -0,0 +1,134 @@
+#include <iostream>
+#include <fstream>
+#include <array>
+#include <map>
+#include <functional>
+#include <string>
+#include <cmath>
+#include <iomanip>
+#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<float(float)>& compute,
+    size_t n_bins
+) {
+    out << "    constexpr std::array<float, " << n_bins << "> " << 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 (<https://doi.org/10.1111/j.1365-3040.2010.02119.x>), 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 (<https://doi.org/10.1111/j.1365-3040.2010.02119.x>), 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 <array>\n\n";
+    out << "namespace lut {\n";
+
+    // Map LUT names to their computation functions
+    std::map<std::string, std::function<float(float)>> 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;
+}

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